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

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

%% 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
[obj data] = generateFakeData(m, con, outputs, lintimes, sd, true, true, seed, opts);
com = compareData(data.OutputList, data.TimesList); % Comparison points

clear ntimes data

%% Prepare for fitting
% Starting (that is, scrambled) parameters are the geometric mean of each
% class Bounds are 0.001X below the smallest and 1000X the largest value in
% each class
opts.LowerBound = zeros(nV,1) + min(m.p(opts.UseParams));
opts.UpperBound = zeros(nV,1) + max(m.p(opts.UseParams));

kbindindex = 1:4;
kcatindex  = 5:2:47;
kmindex    = 6:2:48;

pstart = m.p;
pstart(kbindindex) = geomean(pstart(kbindindex));
opts.LowerBound(kbindindex) = zeros(size(kbindindex)) + min(m.p(kbindindex)) / 1000;
opts.UpperBound(kbindindex) = zeros(size(kbindindex)) + max(m.p(kbindindex)) * 1000;

pstart(kcatindex) = geomean(pstart(kcatindex));
opts.LowerBound(kcatindex) = zeros(size(kcatindex)) + min(m.p(kcatindex)) / 1000;
opts.UpperBound(kcatindex) = zeros(size(kcatindex)) + max(m.p(kcatindex)) * 1000;

pstart(kmindex) = geomean(pstart(kmindex));
opts.LowerBound(kmindex) = zeros(size(kmindex)) + min(m.p(kmindex)) / 1000;
opts.UpperBound(kmindex) = zeros(size(kmindex)) + max(m.p(kmindex)) * 1000;

mstart = m.update(pstart, m.ic);

clear kbindindex kcatindex kmindex pstart

%% Compute predictions of starting model
% !!! 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 cell. In order to recreate results, you must parallelize this
% call on a cluster. !!!

predict = ComparePrediction(m, mstart, posCon, com, opts);

%% Prediction error
diff = zeros(N,3);
ind = 20; % ERK index
for i = 1:N
    truemat = max(reshape(predict{i}, 100,32).', 1); % Floor species concentrations at 1 molecule
    truemat = truemat([1:12,14:29],:); % Remove unchangable species
    fitmat = max(reshape(predict{i,2}, 100,32).',1); % Floor species concentrations at 1 molecule
    fitmat = fitmat([1:12,14:29],:); % Remove unchangable species
    
    nonblank = (truemat > 1) & (fitmat > 1); % Those predictions greater than 1 molecule
    
    % Max relative error in species
    diff(i,1) = max(abs(log(vec(truemat(nonblank) ./ fitmat(nonblank)))));
    
    % ERK
    temp = abs(log(vec(truemat(ind, :) ./ fitmat(ind, :)))); % Extract ERK row
    
    % Max relative error in ERK
    diff(i,2) = max(temp);
    
    % Median relative error in ERK
    diff(i,3) = median(temp(nonblank(ind,:)));
end

% Remove experiments that are all zeros
diff = diff(~isnan(diff(:,3)),:);

clear truemat fitmat i ind nonblank temp

%% Summarize predictions of starting model
fract10all = sum(diff < 0.1, 1) ./ sum(diff > 0, 1);
medianall = median(diff, 1);

%% Fit to nominal experiment
mfit = FitObjective(mstart, 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 predictions of model fit to nominal
% !!! 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 cell. In order to recreate results, you must parallelize this
% call on a cluster. !!!

predict = ComparePrediction(m, mfit, posCon, com, opts);

%% Prediction error
diff = zeros(N,3);
ind = 20; % ERK index
for i = 1:N
    truemat = max(reshape(predict{i}, 100,32).', 1); % Floor species concentrations at 1 molecule
    truemat = truemat([1:12,14:29],:); % Remove unchangable species
    fitmat = max(reshape(predict{i,2}, 100,32).',1); % Floor species concentrations at 1 molecule
    fitmat = fitmat([1:12,14:29],:); % Remove unchangable species
    
    nonblank = (truemat > 1) & (fitmat > 1); % Those predictions greater than 1 molecule
    
    % Max relative error in species
    diff(i,1) = max(abs(log(vec(truemat(nonblank) ./ fitmat(nonblank)))));
    
    % ERK
    temp = abs(log(vec(truemat(ind, :) ./ fitmat(ind, :)))); % Extract ERK row
    
    % Max relative error in ERK
    diff(i,2) = max(temp);
    
    % Median relative error in ERK
    diff(i,3) = median(temp(nonblank(ind,:)));
end

% Remove experiments that are all zeros
diff = diff(~isnan(diff(:,3)),:);

clear truemat fitmat i ind nonblank temp

%% Initialize data structures
% Here are some containers to hold results from each iteration
mfitall = mfit;
Fall = {ObjectiveInformation(mfit, con, obj, opts)};
bestall = zeros(0,1);
fract10all = [fract10all; sum(diff < 0.1, 1) ./ sum(diff > 0, 1)];
medianall = [medianall; median(diff, 1)];

% while(we.FeelLikeIt)

%% Select the next experiment randomly
RandStream.setDefaultStream(RandStream('mt19937ar', 'seed', seed));
best = randi(N);

%% Perform the random experiment
seed = seed + 1;
con = [con; constructExperiment(mfit, tF, u, [], [], [], names{best}, [], [], [], ics(:,best))];
obj = [obj; generateFakeData(m, con(end), outputs, lintimes, sd, true, true, seed, opts)];

%% Fit to all experiments
mfit = FitObjective(mfit, con, obj, opts);

%% Compute predictions
% !!! This is another section only for reference. Simulating 150000
% experiments on a single computer would take a very long time. !!!

predict = ComparePrediction(m, mfit, posCon, com, opts);

%% Prediction error
diff = zeros(N,3);
ind = 20; % ERK index
for i = 1:N
    truemat = max(reshape(predict{i}, 100,32).', 1); % Floor species concentrations at 1 molecule
    truemat = truemat([1:12,14:29],:); % Remove unchangable species
    fitmat = max(reshape(predict{i,2}, 100,32).',1); % Floor species concentrations at 1 molecule
    fitmat = fitmat([1:12,14:29],:); % Remove unchangable species
    
    nonblank = (truemat > 1) & (fitmat > 1); % Those predictions greater than 1 molecule
    
    % Max relative error in species
    diff(i,1) = max(abs(log(vec(truemat(nonblank) ./ fitmat(nonblank)))));
    
    % ERK
    temp = abs(log(vec(truemat(ind, :) ./ fitmat(ind, :)))); % Extract ERK row
    
    % Max relative error in ERK
    diff(i,2) = max(temp);
    
    % Median relative error in ERK
    diff(i,3) = median(temp(nonblank(ind,:)));
end

% Remove experiments that are all zeros
diff = diff(~isnan(diff(:,3)),:);

clear truemat fitmat i ind nonblank temp

%% Store results from iteration
mfitall = [mfitall; mfit];
Fall = [Fall; {ObjectiveInformation(mfit, con, obj, opts)}];
bestall = [bestall; best];
fract10all = [fract10all; sum(diff < 0.1, 1) ./ sum(diff > 0, 1)];
medianall = [medianall; median(diff, 1)];

% endwhile

%%%%% Repeat until we feel like stopping %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%% Plot parameter error
figure % Figure S6
n = length(Fall);
pall = [mstart.p, cat(2, mfitall.p)];
pall = pall(1:48,:);
pall = pall ./ repmat(m.p(1:48),1,n+1); % Compute ratio
labels = [{'Start', 'Nominal'}, cellfun(@num2str, num2cell(1:n-1), 'UniformOutput', false)];

semilogy(1:n+1, pall)
xlabel('Experiment')
set(gca, 'xtick', 1:n+1)
set(gca, 'xticklabel', labels)
ylabel('Error ratio')
yMin = 10.^(floor(log10(min(min(pall(:)), 1/max(pall(:)))) - 0.09));
yMax = 1/yMin;
ylim([yMin yMax])
set(gca, 'ytick', unique([10.^(-(0:2:log10(yMax))), 10.^(0:2:log10(yMax))]))

hline([1.1 0.9], 'k--') % 10% error line

%insetaxes([0.65 0.65 0.33 0.33])
%semilogy(0:1, pall(:,end-1:end))

%hline([1.1 0.9], 'k--')
%hline([1.2 0.8], 'g--')
%hline([1.3 0.7], 'r--')

%xlim([0.7 1.03])
%ylim([0.68 1.33])
%xtick([])
%ytick([])

clear n pall labels yMin yMax

%% Plot prediction scatter plots
subplot(1,3,1)
semilogy(1:size(diff,1), diff(:,1), 'k.', 'MarkerSize', 0.5, 'Color', 'r')
xlim([1, size(diff,1)])
xlabel('Experiments')
set(gca, 'xtick', [])
ylabel('Max relative error')
axis([0 N 0 2])
set(gca, 'ytick', unique([10.^(-(0:7)), 10.^(0:7)]))

hline(0.1, 'k--') % 10% error line

subplot(1,3,2)
semilogy(1:size(diff,1), diff(:,2), 'k.', 'MarkerSize', 0.5, 'Color', 'r')
xlim([1, size(diff,1)])
xlabel('Experiments')
set(gca, 'xtick', [])
ylabel('Max error in ERK')
axis([0 N 0 2])
set(gca, 'ytick', unique([10.^(-(0:7)), 10.^(0:7)]))

hline(0.1, 'k--') % 10% error line

subplot(1,3,3)
semilogy(1:size(diff,1), diff(:,3), 'k.', 'MarkerSize', 0.5, 'Color', 'r')
xlim([1, size(diff,1)])
xlabel('Experiments')
set(gca, 'xtick', [])
ylabel('Median error in ERK')
axis([0 N 0 2])
set(gca, 'ytick', unique([10.^(-(0:7)), 10.^(0:7)]))

hline(0.1, 'k--') % 10% error line

%% Plot prediction histograms
bins = logspace(-8, 1, 50);
count = histc(diff, bins, 1);

subplot(1,3,1) % Figure S7A
h = bar(bins, count(:,1), 'histc');
set(h,'FaceColor', 'r', 'LineStyle', 'none')
set(gca, 'xscale', 'log')
delete(findobj('Marker', '*'))
xlabel('Max relative error')
ylabel('Number of experiments')
xlim([bins(1) bins(end)])
set(gca, 'xtick', unique([10.^(-(0:7)), 10.^(0:7)]))
set(gca, 'ytick', 0:10000:100000)
subfigureletter('A')

vline(0.1, 'k--') % 10% error line

subplot(1,3,2) % Figure S7B
h = bar(bins, count(:,2), 'histc');
set(h,'FaceColor', 'r', 'LineStyle', 'none')
set(gca, 'xscale', 'log')
delete(findobj('Marker', '*'))
xlabel('Max error in ERK')
ylabel('Number of experiments')
xlim([bins(1) bins(end)])
set(gca, 'xtick', unique([10.^(-(0:7)), 10.^(0:7)]))
set(gca, 'ytick', 0:10000:100000)
subfigureletter('B')

vline(0.1, 'k--') % 10% error line

subplot(3,3,9) % Figure S7C
h = bar(bins, count(:,3), 'histc');
set(h,'FaceColor', 'r', 'LineStyle', 'none')
set(gca, 'xscale', 'log')
delete(findobj('Marker', '*'))
xlabel('Median error in ERK')
ylabel('Number of experiments')
xlim([bins(1) bins(end)])
set(gca, 'xtick', unique([10.^(-(0:7)), 10.^(0:7)]))
set(gca, 'ytick', 0:10000:100000)
subfigureletter('C')

vline(0.1, 'k--') % 10% error line

clear bins count h

%% Plot prediction improvements
n = length(Fall);
labels = [{'Start', 'Nominal'}, cellfun(@num2str, num2cell(1:n-1), 'UniformOutput', false)];

subplot(1,2,1) % Figure S8A
plot(0:n, fract10all, 'Linewidth', 2)
axis([0, n, 0, 1])
ylabel('Fraction')
xlabel('Experiments completed')
set(gca, 'xticklabel', labels)
ylim([0 1])
subfigureletter('A')

subplot(1,2,2) % Figure S8B
semilogy(0:n, medianall, 'Linewidth', 2)
ylabel('Median')
xlabel('Experiments completed')
set(gca, 'xticklabel', labels)
ylim([1e-4 1e2])
subfigureletter('B')

clear labels n
