function ex12
close all

%% read data
data.t = xlsread('Data_fit.xlsx', 'A2:A19');
data.m = xlsread('Data_fit.xlsx', 'B2:B19');
data.e = xlsread('Data_fit.xlsx', 'C2:C19');

data.t_complete = xlsread('Data_complete.xlsx', 'A2:A23');
data.m_complete = xlsread('Data_complete.xlsx', 'B2:B23');

%% declare variables
SumResid = nan(6,1); % vector to be filled below with RSS values
popt = cell(6,1); % cell with the optimized parameter values to be filled below
n = size(data.t,1); % number of samples
k = [2,3,2,3,4,3]'; % number of free parameters in each model, order of models same as in exercise sheet
     
SumResid_complete = SumResid;% vector to be filled below with chi^2 values
popt_complete = popt;


%% fit growth laws with a conti. declining growth rate and chi^2

%%% choose in which way you want to fit the model 
%%% RSS: deviation = 1;
%%% chi^2: deviation = 2;
    deviation = 2;

for growth_law = 1:6
    if growth_law == 1 %%% CST
        %%% k0, T0
            parguess = [0.05, 1E3 ];
    elseif growth_law == 2 %%% EXP
        %%% k0, delta, T0
            parguess = [0.07, 0.01, 1E3]; 
    elseif growth_law == 3 %%% PA
        %%% k0, T0
            parguess = [1, 1E3]; 
    elseif growth_law == 4 %%% POW
        %%% k0, delta, T0
            parguess = [0.07, 0.001, 1E2];
    elseif growth_law == 5 %%% EC
        %%% k0, delta, kbasic, T0
            parguess = [0.07, 0.01, 0.01, 1E3];
    elseif growth_law == 6 %%% LOG
        %%% k0, Tset, T0
            parguess = [4E-7, 1E5, 200];
    end
    

%%% lsqnonlin - comment out if you don't want to use this one
    optimizer = 1;
    optimfunc = @(pars) optimfunc_lsqnonlin(pars,data.t, data.m, data.e, growth_law, optimizer, deviation);
    [popt{growth_law},resnorm,~,~,~,~,~] = ...
            lsqnonlin(optimfunc,parguess, zeros(size(parguess)));

%%% particleswarm - comment out if you don't want to use this one
%     optimizer = 2; 
%     optimfunc = @(pars) optimfunc_lsqnonlin(pars,data.t, data.m, data.e, growth_law, optimizer, deviation);
%     options = optimoptions('particleswarm', 'InitialSwarmMatrix', parguess);
%         ub = ones(size(parguess));
%         ub(end) = 1E4;
%     if growth_law == 6; ub(end-1) = 1E6; end
%     [popt{growth_law},resnorm,~] = particleswarm(optimfunc,length(parguess),zeros(size(ub)), ub , options);
        
%%% simulate results and plot
    [tsim,ysim] = simulate(popt{growth_law}, data.t, popt{growth_law}(end), growth_law);
    
    figure(1)
    subplot(3,3,growth_law)
    hold on
    errorbar(data.t,data.m,data.e, '*k')
    plot(tsim,ysim, '-r')
    set(gca,'Yscale', 'log')
    xlabel('time [h]')
    ylabel('area [um^2]')
    
    SumResid(growth_law) = resnorm;
end

    
%% model comparison
subplot(3,3,7);
plot(1:6, SumResid,'*')
ylabel('RSS')
 
%%% AIC
if deviation == 1
    AIC = 2.*k + n.*log(SumResid);
elseif deviation == 2
    AIC = 2.*k + SumResid;
end
subplot(3,3,8); 
plot(1:6, AIC,'*')
ylabel('AIC')
xlabel('models')

%%% BIC
if deviation == 1
	BIC = n.*log(SumResid./n) + k.*log(n);
elseif deviation == 2
    BIC = SumResid + k.*log(n);
end
subplot(3,3,9);
plot(1:6, BIC,'*')
ylabel('BIC')
xlabel('models')

%%% print out
disp('[RSS, AIC, BIC]')
disp([SumResid, AIC, BIC])

[~,I_RSS] = sort(SumResid);
[~,I_AIC] = sort(AIC);
[~,I_BIC] = sort(BIC);
disp('Model ranking [RSS, AIC, BIC]')
disp([I_RSS, I_AIC, I_BIC])


%% comapre to earlier/later time points
for growth_law = 1:6
    [tsim,ysim] = simulate(popt{growth_law}, [data.t;170], popt{growth_law}(end), growth_law);
    [tsim2,ysim2] = simulate(popt{growth_law}, data.t(1):-1:1, popt{growth_law}(end), growth_law);

    tsim = [flip(tsim2);tsim];
    ysim = [flip(ysim2);ysim];

    figure(2)
    subplot(2,3,growth_law)
    hold on
    plot(data.t_complete,data.m_complete,'*k')
    plot(tsim,ysim, '-r')
    set(gca,'Yscale', 'log')
    xlabel('time [h]')
    ylabel('area [um^2]')
    
    ysim_init(growth_law) = ysim(1);
end



%% repeat fits with RSS and complete data set
    deviation = 1;
    
for growth_law = 1:6
    if growth_law == 1 %%% CST
        %%% k0, T0
%             parguess = [0.05, 50 ];  
            parguess = popt{growth_law};
            parguess(end) = ysim_init(growth_law);
    elseif growth_law == 2 %%% EXP
        %%% k0, delta, T0
%             parguess = [0.07, 0.01, 50]; 
            parguess = popt{growth_law};
            parguess(end) = ysim_init(growth_law);
    elseif growth_law == 3 %%% PA
        %%% k0, T0
%             parguess = [1, 50]; 
            parguess = popt{growth_law};
            parguess(end) = ysim_init(growth_law);
    elseif growth_law == 4 %%% POW
        %%% k0, delta, T0
%             parguess = [0.07, 0.001, 50];
            parguess = popt{growth_law};
            parguess(end) = ysim_init(growth_law);
    elseif growth_law == 5 %%% EC
        %%% k0, delta, kbasic, T0
%             parguess = [0.07, 0.01, 0.01, 50];
            parguess = popt{growth_law};
            parguess(end) = ysim_init(growth_law);
    elseif growth_law == 6 %%% LOG
        %%% k0, Tset, T0
%             parguess = [4E-7, 1E5, 50];
%             parguess = [4.2269e-07 1.1778E5    81.651];
            parguess = popt{growth_law};
            parguess(end) = ysim_init(growth_law);
    end
    
%%% lsqnonlin - comment out if you don't want to use this one
%     optimizer = 1;
%     optimfunc = @(pars) optimfunc_lsqnonlin(pars,data.t_complete, data.m_complete, [], growth_law, optimizer, deviation);
%     [popt_complete{growth_law},resnorm,~,~,~,~,~] = ...
%             lsqnonlin(optimfunc,parguess, zeros(size(parguess)), [], []);
        
%%% particleswarm - comment out if you don't want to use this one
    optimizer = 2; 
    optimfunc = @(pars) optimfunc_lsqnonlin(pars,data.t_complete, data.m_complete, [], growth_law, optimizer, deviation);
    options = optimoptions('particleswarm', 'InitialSwarmMatrix', parguess);
        ub = ones(size(parguess));
        ub(end) = 1E4;
    if growth_law == 6; ub(end-1) = 1E6; end
    [popt_complete{growth_law},resnorm,~] = particleswarm(optimfunc,length(parguess),zeros(size(ub)), ub , options);
        
%%% simulate results and plot
    [tsim,ysim] = simulate(popt_complete{growth_law}, data.t_complete, popt_complete{growth_law}(end), growth_law);
    
    figure(3)
    subplot(3,3,growth_law)
    hold on
    plot(data.t_complete,data.m_complete,'*k')
    plot(tsim,ysim, '-r')
    set(gca,'Yscale', 'log')
    xlabel('time [h]')
    ylabel('area [um^2]')
    
    SumResid_complete(growth_law) = resnorm;
end

end

function deviation = optimfunc_lsqnonlin(pars,tspan, Data, Std, growth_law, optimizer, deviation)  
T0 = pars(end);
    [t,y] = simulate(pars, tspan, T0, growth_law);
    
    if deviation == 1
        deviation = (y - Data);
    elseif deviation == 2
        deviation = (y - Data) ./ Std;
    end
    
    if optimizer == 2
        deviation = sum(deviation.^2);
    end
end 

function [t,y] = simulate(pars, tspan, T0,growth_law)
    [t,y] = ode15s(@ode_WingLimb,tspan,T0, [], pars, growth_law);
end

function dy = ode_WingLimb(t,y, pars, growth_law)
    kobs = get_k(t,y, pars, growth_law);
    dy = kobs*y;
end

function kt = get_k(t,y, pars, growth_law)
    if growth_law == 1; % const growth rate
        k0 = pars(1);
        kt = k0;
    elseif growth_law == 2; % exp declining growth rate
        k0 = pars(1);
        delta = pars(2);
        kt = k0 .* exp(-delta .* t);
    elseif growth_law == 3; % 1/T
        k0 = pars(1);
        y0 = pars(2);
        kt = k0 .* y0./(y);
    elseif growth_law == 4; % Power law
        k0 = pars(1);
        delta = pars(2);
        kt = k0 .* 1./t.^delta;
    elseif growth_law == 5; % exp + const
        k0 = pars(1);
        delta = pars(2);
        k1 = pars(3);
        kt = k0 .* exp(-delta .* t) + k1;
    elseif growth_law == 6; % logistic
        k0 = pars(1);
        yset = pars(2);
        kt = k0.*(yset - y);
    end
end


