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

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

%% Fit to nominal experiment
mfit = FitObjective(mstart, con, obj, opts);

%% Initialize data structures
% Here are some containers to hold results from each iteration
mfitall = mfit;
Fall = {ObjectiveInformation(mfit, con, obj, opts)};

%% Create a super experiment
superic = m.ic;

% Over-expression
for i = [7 11 14 16 18 24 26 28]
    superic(i) = superic(i) * 100;
end

% Under-expression
for i = [13 20 30 32]
    superic(i) = superic(i) * 0.01;
end

clear i

%% Perform the super experiment
seed = seed + 1;
con = [con; constructExperiment(m, tF, u, [], [], [], 'Super', [], [], [], superic)];
obj = [obj; generateFakeData(m, con(end), outputs, lintimes, sd, true, true, seed, opts)];

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

%% Store results from iteration
mfitall = [mfitall; mfit];
Fall = [Fall; {ObjectiveInformation(mfit, con, obj, opts)}];

%% Plot eigenspectra
subplot(3,1,3)
uncertaintySpectra(Fall, {'Nominal', 'Combined'})
hline(0.1, 'k--') % 10% error line
set(gca, 'ytick', 10.^(-3:2:8))
subfigureletter('C')

%% Plot parameter error
subplot(3,1,1)
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', 'Super'};

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
