%% 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;

seed = 1337;

%% 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 = [];
con1 = constructExperiment(m, tF, u, [], [], [], 'nominal');
con2 = con1;
con3 = con1;

%% Generate data
outputs = (1:m.nY).';
ntimes = 100;
lintimes = linspace(tF/ntimes, tF, ntimes).'; % 100 evenly spaced time points
sd = @(t,yInd,yVal)(max(0.1*yVal,1)); % Standard deviation of measurement error is 10% floored at 1
obj1 = generateFakeData(m, con1, outputs, lintimes, sd, true, true, seed, opts);
obj2 = obj1;
obj3 = obj1;

clear ntimes data

%% 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

%% Load best experiments from three goal functions
bestall = [
      138082,  25666, 127008;
      101133, 105465,  62575;
      139093,  60244, 104686;
       43764, 138877,  16709;
      127964,  28510, 111644;
        1004,  74045,  10560];
    
Fall1 = {ObjectiveInformation(m, con1, obj1, opts)};
Fall2 = Fall1;
Fall3 = Fall1;

clear Fall

%% Mix and match the experiments
for i = 1:6
    seed = seed + 1;
    con1 = [con1; constructExperiment(m, tF, u, [], [], [], names{bestall(i,mod(i-1,3)+1)}, [], [], [], ics(:,bestall(i,mod(i-1,3)+1)))];
    obj1 = [obj1; generateFakeData(m, con1(end), outputs, lintimes, sd, true, true, seed, opts)];
    Fall1{i+1,1} = ObjectiveInformation(m, con1, obj1, opts);

    con2 = [con2; constructExperiment(m, tF, u, [], [], [], names{bestall(i,mod(i,3)+1)}, [], [], [], ics(:,bestall(i,mod(i,3)+1)))];
    obj2 = [obj2; generateFakeData(m, con2(end), outputs, lintimes, sd, true, true, seed, opts)];
    Fall2{i+1,1} = ObjectiveInformation(m, con2, obj2, opts);

    con3 = [con3; constructExperiment(m, tF, u, [], [], [], names{bestall(i,mod(i+1,3)+1)}, [], [], [], ics(:,bestall(i,mod(i+1,3)+1)))];
    obj3 = [obj3; generateFakeData(m, con3(end), outputs, lintimes, sd, true, true, seed, opts)];
    Fall3{i+1,1} = ObjectiveInformation(m, con3, obj3, opts);
end

clear i

%% Plot eigenspectra
subplot(3,1,1)
uncertaintySpectra(Fall1, [{'Nominal'}, cellfun(@num2str, num2cell(1:length(Fall1)-1), 'UniformOutput', false)])
hline(0.1, 'k--') % 10% error line
set(gca, 'ytick', 10.^(-3:2:8))
subfigureletter('A')

subplot(3,1,2)
uncertaintySpectra(Fall2, [{'Nominal'}, cellfun(@num2str, num2cell(1:length(Fall2)-1), 'UniformOutput', false)])
hline(0.1, 'k--') % 10% error line
set(gca, 'ytick', 10.^(-3:2:8))
subfigureletter('B')

subplot(3,1,3)
uncertaintySpectra(Fall3, [{'Nominal'}, cellfun(@num2str, num2cell(1:length(Fall3)-1), 'UniformOutput', false)])
hline(0.1, 'k--') % 10% error line
set(gca, 'ytick', 10.^(-3:2:8))
subfigureletter('C')
