%% Options
% Global
opts.Order = 2;
opts.UseParams = 1:48;
opts.UseICs = [];
opts.UseModelICs = false;
nVP = length(opts.UseParams);
nVX = length(opts.UseICs);
nV = nVP + nVX;
opts.RelTol = 1e-6;
opts.AbsTol = 1e-9;
opts.Verbose = 1;

% Fitting
opts.TolOptim = 1e-6;
opts.Normalized = true;
opts.UseAdjoint = false;
opts.AdaptAbsTol = true;
opts.MaxIter = 300;
opts.MaxFunEvals = 5000;
opts.Restart = 4;
opts.RestartJump = 0.001;

% Experiment selection
opts.ReturnCount = 30;
opts.MaxGreedySize = 1;
opts.TerminalGoal = -48;
opts.AllowRepeats = false;

goal = @goalParameterSpaceBelowTenPercent;

%% Load model
m = loadModelSbml('EGFNGF_Apgar.xml', [], 1:32, [], opts);

ic = m.ic;
ic(1) = 1000; %EGF
ic(2) = 4560; %NGF
m = m.update(m.p, ic);

clear ic

%% Create nominal experiment
tF = 120;
u = [];
con = constructExperiment(m, tF, u, [], [], [], 'nominal');

%% Generate data
outputs = (1:m.nY).';
ntimes = 20; % 100, 80, 40, 20, 10, 5
lintimes = linspace(tF/ntimes, tF, ntimes).'; % ntimes evenly spaced time points
sd = @(t,yInd,yVal)(max(0.1*yVal,1)); % Standard deviation of measurement error is 10% floored at 1
seed = 1337;
[obj data] = generateFakeData(m, con, outputs, lintimes, sd, true, true, seed, opts);
com = compareData(data.OutputList, data.TimesList); % Comparison points

clear ntimes data

%% Initial Information
F = ObjectiveInformation(m, con, obj, opts);

%% Create possible experiments
% Possible doses
EGFdoses =        [1e7, 1e5, 1e3, 10, 0];
NGFdoses = 4.56 * [1e7, 1e5, 1e3, 10, 0];

% Species to be knocked-up or down
xchange = [3, 5, 7, 9, 11, 13, 14, 16, 18, 20, 22, 24, 26, 28, 30, 31, 32];
knocks = [100 0.01];
knockchar = '+-';

N = 5*5*(1+nchoosek(17,1)*2+nchoosek(17,2)*2^2+nchoosek(17,3)*2^3);
names = cell(N,1);
ics = zeros(m.nX,N);

index = 1;

% No knock-downs or over-expression
for iEGF = 1:length(EGFdoses)
    for iNGF = 1:length(NGFdoses)
        ic = m.ic;
        ic(1) = EGFdoses(iEGF);
        ic(2) = NGFdoses(iNGF);
        names{index} = sprintf('E=%d N=%d', EGFdoses(iEGF), NGFdoses(iNGF));
        ics(:,index) = ic;
        index = index + 1;
    end
end

% Single experiments
for iEGF = 1:length(EGFdoses)
    for iNGF = 1:length(NGFdoses)
        for i = 1:length(xchange)
            ic = m.ic;
            ic(1) = EGFdoses(iEGF);
            ic(2) = NGFdoses(iNGF);
            ic(xchange(i)) = ic(xchange(i))*knocks(1);
            names{index} = sprintf('E=%d N=%d %d+', EGFdoses(iEGF), NGFdoses(iNGF), xchange(i));
            ics(:,index) = ic;
            index = index + 1;
            
            ic = m.ic;
            ic(1) = EGFdoses(iEGF);
            ic(2) = NGFdoses(iNGF);
            ic(xchange(i)) = ic(xchange(i))*knocks(2);
            names{index} = sprintf('E=%d N=%d %d+', EGFdoses(iEGF), NGFdoses(iNGF), xchange(i));
            ics(:,index) = ic;
            index = index + 1;
        end
    end
end

% Double experiments
doubles = xchange(combinator(17, 2, 'c', 'n'));
N2 = size(doubles, 1);
for iEGF = 1:length(EGFdoses)
    for iNGF = 1:length(NGFdoses)
        for i = 1:N2
            for j = 1:2
                for k = 1:2
                    ic = m.ic;
                    ic(1) = EGFdoses(iEGF);
                    ic(2) = NGFdoses(iNGF);
                    ic(doubles(i,1)) = ic(doubles(i,1))*knocks(j);
                    ic(doubles(i,2)) = ic(doubles(i,2))*knocks(k);
                    names{index} = sprintf('E=%d N=%d %d%c %d%c', EGFdoses(iEGF), NGFdoses(iNGF), doubles(i,1), knockchar(j), doubles(i,2), knockchar(k));
                    ics(:,index) = ic;
                    index = index + 1;
                end
            end
        end
    end
end

% Triple experiments
triples = xchange(combinator(17, 3, 'c', 'n'));
N3 = size(triples, 1);
for iEGF = 1:length(EGFdoses)
    for iNGF = 1:length(NGFdoses)
        for i = 1:N3
            for j = 1:2
                for k = 1:2
                    for l = 1:2
                        ic = m.ic;
                        ic(1) = EGFdoses(iEGF);
                        ic(2) = NGFdoses(iNGF);
                        ic(triples(i,1)) = ic(triples(i,1))*knocks(j);
                        ic(triples(i,2)) = ic(triples(i,2))*knocks(k);
                        ic(triples(i,3)) = ic(triples(i,3))*knocks(l);
                        names{index} = sprintf('E=%d N=%d %d%c %d%c %d%c', EGFdoses(iEGF), NGFdoses(iNGF), triples(i,1), knockchar(j), triples(i,2), knockchar(k), triples(i,3), knockchar(l));
                        ics(:,index) = ic;
                        index = index + 1;
                    end
                end
            end
        end
    end
end

clear EGFdoses NGFdoses N2 N3 iEGF iNGF i j k l ic index xchange knocks knockchar doubles triples

posCon = struct('name', names); % Place-holder experimental conditions
posObj(N,1) = struct; % Place-holder objective functions

%% Compute Fisher information matrix for each experiment
% !!! This section is merely for reference. Do not attempt to run on a
% single machine. It would take months for a multicore desktop computer to
% run this call. In order to recreate results, you must parallelize this
% call on a cluster. !!!

Fs = cell(N,1);
for i = 1:N
    posConi = constructExperiment(m, tF, u, [], [], [], names{i}, [], [], [], ics{i});
    posObji = generateFakeData(m, posConi, outputs, lintimes, sd, true, false, seed, opts);
    Fs{i} = ObjectiveInformation(m, posConi, posObji, opts);
end

clear posConi posObji i

%% Optimal experimental design
best = BestParameterExperiment(mfit, con, obj, posCon, posObj, goal, opts, F, Fs);

%% Plot
results = [
    100, 6;
    80, 5;
    40, 6;
    20, 7;
    10, 10;
    5, 12];

plot(results(:,1), results(:,2), 'k-o', 'Linewidth', 2)
xlim([0,100])
ylim([0,max(results(:,2))])
xlabel('Number of time points')
ylabel('Number of experiments required')
