% Copyright (c) 2015, Domenico L. Gatti
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without 
% modification, are permitted provided that the following conditions are 
% met:
% 
%     * Redistributions of source code must retain the above copyright 
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright 
%       notice, this list of conditions and the following disclaimer in 
%       the documentation and/or other materials provided with the 
%       distribution
%       
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
% 

%% Two Substrates Bi-Uni Michaelis Menten Kinetics: Simulation and Analysis.

% Add to the path the enhanced scripts for drawing lines at specified
% coordinates.
addpath(genpath('../GENERAL_SCRIPTS_FUNCTIONS'));
addpath(genpath('../TOOLBOXES/ENZYME_KINETICS'));

global m1
global Time_points ntimes
global s1_cons s1_conc_vec s2_cons s2_conc_vec
global C C_init E E_init S1 S2 ES1 ES2 P kon1 koff1 kcat kon2 koff2
global Product_vector Time_vector
global Product_mat Time_mat

% Load the project. Here we can use:
% 1. The correct model:
sbioloadproject('../TOOLBOXES/ENZYME_KINETICS/Bi_Uni_Michaelis_Menten');
% 2. The incorrect model:
% Load the alternative project with switched order of substrate binding
% sbioloadproject('../TOOLBOXES/ENZYME_KINETICS/Bi_Uni_Michaelis_Menten_Switch_Substrates');

%%
% Get information about the model.
sbioselect(m1,'Type','compartment')
sbioselect(m1,'Type','species')
sbioselect(m1,'Type','parameter')
sbioselect(m1,'Type','reaction')

% We also want to know the differential equations used in the model. This
% information may be very usefult for special applications (e.g., global
% fitting of multiple curves)
equations = getequations(m1)

% First, let's get the configuration of the simulation.
configset_m1 = getconfigset(m1);
get(configset_m1)

%%

Time_mat = dlmread('../TOOLBOXES/ENZYME_KINETICS/Bi_Uni_Michaelis_Menten_Time_mat.txt');
Substrate_mat = dlmread('../TOOLBOXES/ENZYME_KINETICS/Bi_Uni_Michaelis_Menten_Substrate_mat.txt');
Product_mat = dlmread('../TOOLBOXES/ENZYME_KINETICS/Bi_Uni_Michaelis_Menten_Product_mat.txt');

%% 

% set(configset_m1, 'StopTime', stop_time);
set(configset_m1,'SolverType','ode15s');
set(configset_m1.SolverOptions, 'AbsoluteTolerance', 1.0e-9);

% And we configure the solver so that states are logged only at
% specific times.

Time_points = Time_mat(:,1);
ntimes = size(Time_points,1);

set(configset_m1.SolverOptions, 'OutputTimes',Time_points);
get(configset_m1)

% It looks OK. Now we want to see what happens when we change some of the
% parameters. First, we will extract the parameters from the model and save
% them as variables in the Workspace.
C = sbioselect(m1,'Name','Cell');
S1 = sbioselect(m1,'Name','S1');
E = sbioselect(m1,'Name','E');
S2 = sbioselect(m1,'Name','S2');
ES1 = sbioselect(m1,'Name','ES1');
ES2 = sbioselect(m1,'Name','ES2');
P = sbioselect(m1,'Name','P');
kon1 = sbioselect(m1,'Name','kon1');
koff1 = sbioselect(m1,'Name','koff1');
kon2 = sbioselect(m1,'Name','kon2');
koff2 = sbioselect(m1,'Name','koff2');
kcat = sbioselect(m1,'Name','kcat');

% We store the values in the model, which were used to generate the
% simulated experimental data with noise.
koff1_orig = koff1;
koff2_orig = koff2;
kcat_orig = kcat;

% We store the initial value for the enzyme and the cell volume.
C_init = C.capacity;
E_init = E.InitialAmount;

%% Global fit
close all

s1_conc_vec = Substrate_mat(:,1);
s2_conc_vec = Substrate_mat(1,2:2:20);
s1_cons = length(s1_conc_vec);
s2_cons = length(s2_conc_vec);

% Alternatively we could fit globally all the progress curves with a single
% set of parameters. For this purpose we will first convert all the
% progress curves into a single consecutive progress curve. 

Product_vector = Product_mat(:);
Time_vector = repmat(Time_mat(:),s2_cons,1);
plot(Time_vector,Product_vector,'+')

% Check a few curves: i.e. 5th S2 concentration
figure;plot(Time_mat(:,1),Product_mat(:,41:50),'--+');
xlabel('Time (s)  ')
ylabel('[Product]  ')

%% Choice of parameters to refine.
% We need to set some initial value for the refinement of parameters: we
% will refine koff1, kcat, and koff2. 
p1 = 50 % koff1.Value;
p2 = 0.3 % kcat.Value;
p3 = 150 % koff2.Value;

pin=[p1;p2;p3];

% We also need to reinitialize the initial value for the enzyme
% concentration and the cell volume.
E.InitialAmount = E_init;
C.capacity = C_init;

%% Non-linear global fit of all the progress curves.

% lsqnonlin (Optimization Toolbox)

% Here we show how to set the plotting of the uptated fit for one set of
% reactions. Global values are passed to the plotting function
% 'optimplot_tracefit.m'.
global trace lasttrace
trace = 4; lasttrace = 10;
% Here we set the options
options = ...
    optimoptions('lsqnonlin','Display','iter','FinDiffType','central',...
    'TolFun',1e-6,'TolX',1e-6,'PlotFcns',@optimplot_tracefit);
% Here we run the global fit
fun = @(pin) (bi_uni_fit(pin,Time_vector) - Product_vector);
[FP1,sos1,R1,exitflag,output,lambda,J1] = lsqnonlin(fun,pin,[],[],options); 

%% Alternative Matlab solutions
% lsqcurvefit (Optimization Toolbox)
options = ...
    optimoptions('lsqcurvefit','Display','iter','FinDiffType','central',...
    'TolFun',1e-6,'TolX',1e-6,'PlotFcns',@optimplot_tracefit);
[FP1,sos1,R1,exitflag,output,lambda,J1] = ...
    lsqcurvefit(@bi_uni_fit,pin,Time_vector,Product_vector,[],[],options);

% nlinfit (Statistics Toolbox) [very fast, but smal radius of convergence 
% and prone to fail if the guess is very far from real values]
options=statset('TolX',1e-6,'TolFun',1e-6,'Display','iter',...
    'RobustWgtFun','','UseParallel',true);
[FP1,R1,J1,Cov1,MSE1] = ...
    nlinfit(Time_vector,Product_vector,'bi_uni_fit',pin,options);
[Corr1,sigma1] = corrcov(Cov1);
sos1 = R1'*R1;
[fcurve1,delta1] = nlpredci('bi_uni_fit',Time_vector,FP1,R1,'Covar',Cov1);
FP1ci = nlparci(FP1,R1,'covar',Cov1)

% The following alternative syntax gives the same result as the previous 
% two lines.
% [fcurve1,delta1] = nlpredci('fit5fitM',lt,FP1,R1,'jacobian',J1);
% FPci = nlparci(FP1,R1,'jacobian',J1);

% fminsearch (Optimization Toolbox) [widest possible radius of convergence,
% but potentially very slow to converge, as it does not use derivatives. 
% Consequently there is no Jacobian. Use as first step if the other methods
% fail.]
modelfun = @(pin) sum((bi_uni_fit(pin,Time_vector) - Product_vector).^2);
options = optimset('Display','iter', 'TolFun',1e-6, 'TolX',1e-6,...
    'PlotFcns',@optimplot_tracefit);
[FP1,fval,exitflag,output] = fminsearch(modelfun,pin,options);

FP1

%% Regeneration of the calculated individual progress curves.

calc_product_mat = reshape(fcurve1,ntimes,s1_cons,s2_cons);

% Check a few curves.
close all
figure;plot(Time_mat(:,1),Product_mat(:,41:50),'+');
colors = get(gca,'ColorOrder');
hold on
for i = 1:10
plot(Time_mat(:,1),calc_product_mat(:,40+i),'-','Color', colors(i,:));
end
xlabel('Time (s)  ')
ylabel('[Product]  ')

%% Step-by-step non-linear least-squares

p1 = 50 % koff1.Value;
p2 = 0.3 % kcat.Value;
p3 = 150 % koff2.Value;
pin=[p1;p2;p3];

nobs = length(Product_vector);
 
u = pin;
nvars = size(u,1);
u1 = u(1);
u2 = u(2);
u3 = u(3);
 
tolerance = 1E-3;
delta_sos = tolerance*10;
diff_g_sos = tolerance*10;

Product_sim_vector = bi_uni_fit(pin,Time_vector);

g = Product_sim_vector - Product_vector;

trust = 4E-1;
niter = 0;

% for housekeeping table
grad = zeros(nobs,nvars,1000);
grad_sos = zeros(1000,nvars);
res_sos = zeros(1000,1);
delta_u_sos = zeros(1000,1);
delta_g_sos = zeros(1000,1);

du = eps^(1/3);

figure;plot(Time_mat(:,1),Product_mat(:,41:50),'+');
colors = get(gca,'ColorOrder');
colors = [colors ; colors(1:3,:)];
xlabel('Time (s)  ')
ylabel('[Product]  ')
hold on

while diff_g_sos > tolerance 
 
niter = niter + 1;
  
g_minus = g;
  
 
u1_plus = (1+du)*u1;
u2_plus = (1+du)*u2;
u3_plus = (1+du)*u3;

pin_init = pin;
pin(1) = u1_plus;
Product_sim_vector = bi_uni_fit(pin,Time_vector);
g1_plus = Product_sim_vector - Product_vector;
pin = pin_init;
pin(2) = u2_plus;
Product_sim_vector = bi_uni_fit(pin,Time_vector);
g2_plus = Product_sim_vector - Product_vector;
pin = pin_init;
pin(3) = u3_plus;
Product_sim_vector = bi_uni_fit(pin,Time_vector);
g3_plus = Product_sim_vector - Product_vector;


u1_minus = (1-du)*u1;
u2_minus = (1-du)*u2;
u3_minus = (1-du)*u3;

pin_init = pin;
pin(1) = u1_minus;
Product_sim_vector = bi_uni_fit(pin,Time_vector);
g1_minus = Product_sim_vector - Product_vector;
pin = pin_init;
pin(2) = u2_minus;
Product_sim_vector = bi_uni_fit(pin,Time_vector);
g2_minus = Product_sim_vector - Product_vector;
pin = pin_init;
pin(3) = u3_minus;
Product_sim_vector = bi_uni_fit(pin,Time_vector);
g3_minus = Product_sim_vector - Product_vector;

 
J = [(g1_plus-g1_minus)/(u1_plus-u1_minus) ...
     (g2_plus-g2_minus)/(u2_plus-u2_minus) ...
     (g3_plus-g3_minus)/(u3_plus-u3_minus)];

h = -g;

[Q,R] = qr(J,0);
delta_u = R\Q'*h;
 
delta_u = delta_u*trust;
u = u + delta_u;
u1 = u(1);
u2 = u(2);
u3 = u(3);

pin = [u1;u2;u3];
Product_sim_vector = bi_uni_fit(pin,Time_vector);
g = Product_sim_vector - Product_vector;
  
sos_minus = g_minus'*g_minus;
sos = g'*g;
delta_sos = abs(sos-sos_minus);
diff_g_sos = (g-g_minus)'*(g-g_minus);

% Housekeeping table
grad(:,:,niter) = J;
grad_sos(niter,:) = diag(J'*J);
res_sos(niter) = sos;
delta_u_sos(niter) = delta_u'*delta_u;
delta_g_sos(niter) = diff_g_sos;

% Follow one set of curves
calc_product_mat = reshape(Product_sim_vector,ntimes,s1_cons,s2_cons);

if niter > 1
for ntrace = 1:10
    fieldname = ['f' num2str(ntrace)];
    delete(eval(['p.' fieldname])); 
end
end

for ntrace = 1:10
    fieldname = ['f' num2str(ntrace)];
    p.(fieldname) = plot(Time_mat(:,1),calc_product_mat(:,ntrace,5),'-','Color', colors(ntrace,:));
    drawnow expose; 
end

% Display table
disp('   iter   res_sos      delta_u_sos    delta_sos      delta_g_sos')
fprintf('%6i %10.4f %14.4f %14f %14.4f\n',...
    [niter res_sos(niter) delta_u_sos(niter) delta_sos delta_g_sos(niter)])

end

grad_table = grad(:,:,1:niter);
sos_table = [(1:niter)' grad_sos(1:niter,:) res_sos(1:niter) ...
    delta_u_sos(1:niter) delta_g_sos(1:niter)./delta_u_sos(1:niter)];

figure; plot(Product_vector,'+b')
hold on
plot(Product_sim_vector,'or')

%%
Params = pin

sse = sos
MSE = sos/(nobs-nvars)
RMSE = sqrt(MSE)
 
u_cov = inv(J'*J).*MSE;

[Corr,sigma] = corrcov(u_cov);
  
cov_inv = (J'*J)./MSE;
 
Rho = zeros(nvars,nvars);
 
    for i = 1:nvars
        for j = i:nvars
            Rho(i,j) = cov_inv(i,j)/sqrt(cov_inv(i,i)*cov_inv(j,j));
            Rho(j,i) = Rho(i,j);
        end
    end
     
Corr
Rho
 
Conf_95 = [u-1.96*sigma u+1.96*sigma]
 
c_Product_vector = Product_vector - mean(Product_vector);
sst = c_Product_vector'*c_Product_vector;
rsquare = 1-sse/sst

niter
 
