% Copyright (c) 2017, 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.
%
%% General dependencies
% We always start from the CODE directory and we add to the path
% subdirectories containing various tools described in the book chapters.
addpath(genpath('../GENERAL_SCRIPTS_FUNCTIONS'));
% addpath(genpath('../DATABASE'));
% addpath(genpath('../TOOLBOXES'));

%% CHAPTER 8: Non-linear least squares
 
t = [0 5 10 15 20 25 30 35 40 45 50 55 60]'
nt = length(t)
 
yvec = [0.9669 0.5627 0.4608 0.2979 0.3493 0.4414 0.2387 0.2586 ...
           0.0988 0.0896 0.1247 0.0378 0.03031]';
       
OLS_plot = figure;       
set(gcf,'Unit','Normalized','Position',[0 0.2 0.5 0.8])
plot(t,yvec,'sb','MarkerSize',30,'LineWidth',2)
xlim([-1 1.02*max(t)])
xlabel('Time (minutes)')
ylabel('[compound X] (microM)')
hold on

u1 = yvec(1);
u2 = log(1/2)/10;

yvec1 = u1*exp(u2*t);
resid_vec1 = yvec1 - yvec;
sse1 = resid_vec1'*resid_vec1
plot(t,yvec1,'-r','LineWidth',3)

u1 = .95*u1;
u2 = .95*u2;
yvec2 = u1*exp(u2*t);
resid_vec2 = yvec2 - yvec;
sse2 = resid_vec2'*resid_vec2
plot(t,yvec2,'-g','LineWidth',3)

u2 = .90*u2;
yvec3 = u1*exp(u2*t);
resid_vec3 = yvec3 - yvec;
sse3 = resid_vec3'*resid_vec3
plot(t,yvec3,'-c','LineWidth',3)

u2 = .90*u2;
yvec4 = u1*exp(u2*t);
resid_vec4 = yvec4 - yvec;
sse4 = resid_vec4'*resid_vec4
plot(t,yvec4,'-m','LineWidth',3)

legend({'yvec','yvec_1','yvec_2','yvec_3','yvec_4'},'FontSize',30)

% Grid search
u1 = [0.7:0.02:1.0]; 
u2 = [-.07:0.001:-.03];
gsq = zeros(length(u1),length(u2));
for i = 1:length(u1)
    for j = 1:length(u2)
        g = yvec - u1(i)*exp(u2(j)*t);
        gsq(i,j) = g'*g;
    end
end
residfun = figure;
[X,Y] = meshgrid(u1,u2);
surf(X,Y,gsq')
ylabel('u2')
xlabel('u1')
zlabel('g^2')
box on
grid on
% colormap default
% colormap jet
% colorbar

% Or alternatively
heatmap_resid_funct = figure'
imagesc(gsq) 
set(gca,'YDir','normal');

%%
clear, clc, close all

t = [0 5 10 15 20 25 30 35 40 45 50 55 60]';
nt = length(t);

xvec = t; % time points
yvec = [0.9669 0.5627 0.4608 0.2979 0.3493 0.4414 0.2387 0.2586 ...
           0.0988 0.0896 0.1247 0.0378 0.03031]'; % data points
nobs = length(yvec);

lin_nonlin_lsq = figure;
set(gcf,'Unit','Normalized','Position',[0 0.2 0.5 0.8])
plot(t,yvec,'sb','MarkerSize',30,'LineWidth',2)
xlim([-1 1.02*max(t)])
xlabel('Time (minutes)')
ylabel('[compound X] (microM)')
hold on


% Linear least squares
A = [ones(nt,1) xvec];
u_l = A\log(yvec);
u_l(1) = exp(u_l(1));
yvec_l = u_l(1)*exp(u_l(2)*t);
plot(t,yvec_l,'--r','LineWidth',3);
sse_l = (yvec_l-yvec)'*(yvec_l-yvec);
mse_l = sse_l/(nobs-length(u_l));
sigma_l = sqrt(mse_l);
c_yvec = yvec - mean(yvec);
sst = c_yvec'*c_yvec;
rsquare_l = 1-sse_l/sst;
u_cov_l = mse_l*inv(A'*A);
[Corr_l,sigma_u_l] = corrcov(u_cov_l);   
Conf_95_l = [u_l-1.96*sigma_u_l u_l+1.96*sigma_u_l];
legend({'Exper. data','linear fit'},'FontSize',30)

%%
u = [yvec(1);log(1/2)/xvec(3)];
nvars = size(u,1);
u1 = u(1);
u2 = u(2);
 
tolerance = 1E-8;
delta_sos = tolerance*10;

g = u1*exp(u2*xvec) - yvec;
niter = 0;

% for housekeeping
grad = zeros(nobs,nvars,1000);
grad_norm = zeros(1000,1);
res_sos = zeros(1000,1);
res_norm = zeros(1000,1);
delta_u_norm = zeros(1000,1);
iter_error = zeros(1000,1);

while delta_sos > tolerance 
 
niter = niter + 1;
  
g_old = g;
  
du1 = u1*eps^(1/3);
du2 = u2*eps^(1/3);

g1 = @(u1) u1*exp(u2*xvec) - yvec;
g2 = @(u2) u1*exp(u2*xvec) - yvec;
J1 = (g1(u1+du1)-g1(u1-du1))/(2*du1);
J2 = (g2(u2+du2)-g2(u2-du2))/(2*du2); 
J = [J1 J2]; 

% Alternative Jacobian using complex step differentiation: in this case
% the stepsize can be machine tolerance.
%
% du = eps;
% g1 = @(u1) u1*exp(u2*xvec) - yvec;
% g2 = @(u2) u1*exp(u2*xvec) - yvec;
% J1 = imag(g1(u1+du*1i)/du);
% J2 = imag(g2(u2+du*1i)/du);
% J = [J1 J2]; 

h = -g;
 
[Q,R] = qr(J,0);
delta_u = R\Q'*h;
 
u = u + delta_u;
u1 = u(1);
u2 = u(2);
 
g = u1*exp(u2*xvec) - yvec;
  
sos_old = g_old'*g_old;
sos = g'*g;
delta_sos = abs(sos-sos_old);

% Table
grad(:,:,niter) = J;
grad_norm(niter) = norm(J'*g);
res_sos(niter) = sos;
res_norm(niter) = sqrt(sos);
delta_u_norm(niter) = norm(delta_u);
iter_error(niter) = (J*delta_u + g_old)'*(J*delta_u + g_old);
     
end

yvec_nl = u1*exp(u2*xvec);
plot(xvec,yvec_nl,'-b','LineWidth',3)
legend({'Exper. data','linear fit','non-linear fit'},'FontSize',30)
 
mse = sos/(nobs-nvars);
sigma = sqrt(mse);
u_cov = mse*inv(J'*J);
[Corr,sigma_u] = corrcov(u_cov);
Conf_95 = [u-1.96*sigma_u u+1.96*sigma_u];
rsquare = 1-sos/sst;

% Plot confidence interval
u_minus = u-1.96*sigma_u;
u_plus = u+1.96*sigma_u;
yvec_nl_lci = u_minus(1)*exp(u_minus(2)*xvec);
yvec_nl_uci = u_plus(1)*exp(u_plus(2)*xvec);

% yvec_nl_luci = u_minus(1)*exp(u_plus(2)*xvec);
% yvec_nl_ulci = u_plus(1)*exp(u_minus(2)*xvec);
% plot(xvec,yvec_nl_luci,'-c','LineWidth',2)
% plot(xvec,yvec_nl_ulci,'--c','LineWidth',2)

% Plot strict interpretation of confidence intervals
yvec_nl_lci1m = u(1)*exp(u_minus(2)*xvec);
yvec_nl_lci2m = u(1)*exp(u_plus(2)*xvec);
% plot(xvec,yvec_nl_lci1m,'-g','LineWidth',2)
% plot(xvec,yvec_nl_lci2m,'--g','LineWidth',2)
yvec_nl_uci1m = u_minus(1)*exp(u(2)*xvec);
yvec_nl_uci2m = u_plus(1)*exp(u(2)*xvec);
% plot(xvec,yvec_nl_uci1m,'-r','LineWidth',2)
% plot(xvec,yvec_nl_uci2m,'--r','LineWidth',2)

yvec_high = max([yvec_nl_lci2m';yvec_nl_uci2m']);
yvec_low = min([yvec_nl_lci1m';yvec_nl_uci1m']);
plot(xvec,yvec_nl_lci,':g','LineWidth',2)
plot(xvec,yvec_low,':k','LineWidth',2)
plot(xvec,yvec_nl_uci,':g','LineWidth',2)
plot(xvec,yvec_high,':k','LineWidth',2)
legend({'Exper. data','linear fit','non-linear fit',...
    '95% conf. int','strict 95% c.i.'},'FontSize',30)

%% Table 1
cycles = [1 : niter]
cnames = {'Iterations' 'Residual_sos' 'Step_norm' 'Gradient_norm'};
NLLS_T1 = table((cycles)',res_sos(cycles),delta_u_norm(cycles),grad_norm(cycles),...
    'VariableNames',cnames)

%% Table 2
cnames = {'Parameter' 'LLS' 'NLLS' };
NLLS_T2 = table(categorical({'U1';'U2';'SSE';'MSE';'SIGMA';'RSQUARE'}),...
    [u_l(1);u_l(2);sse_l;mse_l;sigma_l;rsquare_l],...
    [u(1);u(2);sos;mse;sigma;rsquare],'VariableNames',cnames)

%% Built-in MATLAB functions for non-linear least-squares
u_zero = [yvec(1);log(1/2)/xvec(3)];

options = ...
    optimoptions('lsqnonlin','Display','iter','FinDiffType','central',...
    'TolFun',1e-8,'TolX',1e-8);
[u,sos,res,flag,output,lambda,J] = ...
    lsqnonlin(@(u)u(1)*exp(u(2)*t) - yvec,u_zero,[],[],options); 
[u,sos,res,~,~,~,J] = ...
    lsqnonlin(@(u)u(1)*exp(u(2)*t) - yvec,u_zero,[],[],options); 

options = ...
    optimoptions('lsqcurvefit','Display','iter','FinDiffType','central',...
    'TolFun',1e-8,'TolX',1e-8);
[u,sos,res,~,~,~,J] = ...
    lsqcurvefit(@(u,t)u(1)*exp(u(2)*t),u_zero,t,yvec,[],[],options);

modelfun = @(u,t)u(1)*exp(u(2)*t);
options = statset('Display','iter','TolTypeFun','abs',...
    'TolFun',1e-8,'TolTypeX','abs','TolX',1e-8);
 [u,R,J,CovB,MSE,ErrorModelInfo] = ...
     nlinfit(t,yvec,modelfun,u_zero,options);    
ci = nlparci(u,R,'Jacobian',J,'alpha',0.05);
[Ypred,delta] = nlpredci(modelfun,t,u,R,'Jacobian',J, 'alpha',0.05)
% plot(xvec,yvec_nl+delta,':y','LineWidth',2)
% plot(xvec,yvec_nl-delta,':y','LineWidth',2)
    nlintool(t,yvec,modelfun,u_zero,0.05,'time','[Y]');    
    
fo = fitoptions('Method','NonlinearLeastSquares','StartPoint',u_zero, ...
    'Algorithm','Trust-Region','Display','iter','TolFun',1e-8,'TolX',1e-8);
f = fittype('a*exp(b*x)','options',fo);
[decay_model,GOF,output] = fit(t,yvec,f);

modelfun = @(u) sum((u(1)*exp(u(2)*xvec)-yvec).^2);
options = optimset('Display','iter', 'TolFun',1e-8, 'TolX',1e-8);
[u,fval,exitflag,output] = fminsearch(modelfun,u_zero,options);

