clear();
funcprot(0);

exp_data_file = './1782_data.txt' // Experimental data 
disp("Data: " + exp_data_file)

[exp_data, text] = fscanfMat(exp_data_file);
t_exp = exp_data(:,1);
y_exp = exp_data(:,2:6);
t0 = t_exp(1);
t_end = t_exp($);
C0 = [1 0 0 0 0]; // initial concentrations of the reagents
reagent_name = ['I','II','III','IV','tBuBr']

exp_err = 0.01*ones(y_exp);

exp_data_HBr = exp_data;
t_exp_HBr = exp_data_HBr(:,1);
HBr_exp = exp_data_HBr(:,7);
t0_HBr = t_exp_HBr(1);
t_end_HBr = t_exp_HBr($);
HBr_name = 'HBr'

rate_start = [2.7e-4  8.1e-2  4.8e-4  4.1e-5  1.5e-4  3.8e-3]; //1782
//rate_start = [5e-5 6e-2 4e-5 3e-6 2e-5 9e-4] //1777
//rate_start = [5.4e-05 1.5e-01 4.2e-05 4.2e-06 8.2e-05 1.3e-03] //1780
//rate_start = [1.7e-05 2.7e-04 3.5e-05 2.7e-06 2.3e-05 1.3e-03] //1781
//rate_start = [1e-3 0.13 4e-3 4e-4 7e-4 1e-2] //1784
//rate_start = [0.011  0.43  0.050  0.0056  0.0048  0.046] //1790
rate_name = ['k13','K','k34','k12','k24','kBu']
//rate_name = ['k13','k31','k34','k12','k24','kBu']
rate_fix = [0 0 0 0 0 0]


// Errors handling:
// W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling. 
// Numerical Recipes in C. Cambridge University Press, Cambridge: 1988

// Constant chi-square (parameters variation): p.551
// Monte Carlo (bootstrap): p. 548

nsigma = 1; // 68.3%
//nsigma = 2; // 95.4%
nsigma = 3; // 99.7%
//hessian = 1;
variation = 1;
montecarlo = 200;

///////////////////////////////////////////////////////////////////////////////
m = length(y_exp)
N = length(t_exp) // Number of measurements
N = length(y_exp) // Number of measurements
Y_EXP = y_exp;
rate_n = length(rate_start)
rate_n_opt = 0
rate_opt_i = []
for i=1:rate_n
    if rate_fix(i)==0 then
        rate_n_opt = rate_n_opt + 1
        rate_opt_i = [rate_opt_i, i]
    end
end
M = rate_n_opt

global count;
count = 0;

if isdef("HBr_exp") then
    // HBr ////////////////////////////////////////////////////////////////////////
    HBr_n = 6; // number of breakpoints
    HBr_x = linspace(t0_HBr,t_end_HBr,HBr_n)';
    [HBr_y, HBr_d] = lsq_splin(t_exp_HBr, HBr_exp, HBr_x);  // compute the spline use equal weights
    //HBr_ys = interp(t, HBr_x, HBr_y, HBr_d);
    //HBr_SSD = sum((HBr_exp-HBr_ys).^2);
    ///////////////////////////////////////////////////////////////////////////////
    
    function HBr=C_HBr(t)
        //HBr = -2.770733E-16*t.^4 + 6.554356E-12*t.^3 - 5.872748E-08*t.^2 + 3.172389E-04*t - 3.979956E-03;
        //HBr = interpln([t_exp';HBr_exp'],t);
        //HBr = interp1(t_exp',HBr_exp',t,'spline');
        HBr = interp(t, HBr_x, HBr_y, HBr_d);
    endfunction
    HBr_SSD = sum((HBr_exp-C_HBr(t_exp_HBr)).^2)
    mprintf("HBr SSD: %e\n\n", HBr_SSD);
end

function dy=myModel(t,y,rate)
    HBr = C_HBr(t);
    
    kk = rate_start
    for i=1:length(rate_opt_i)
        kk(rate_opt_i(i)) = rate(i)
    end

    k13 = kk(1) 
    //b31 = kk(2)*HBr  // kk(2)=k31
    b31 = k13/kk(2)*HBr  // kk(2)=k13/k31 equilibrium constsnt as optimized parameter
    k34 = kk(3) 
    k12 = kk(4) 
    k24 = kk(5) 
//    b42 = kk(6)*HBr  // kk(6)=k42
//    b42 = k24/kk(6)*HBr  // kk(6)=k24/k42
    kBu = kk(6)

    k = [-k13-k12,      0,       +b31,  0,    0
             +k12,   -k24,          0,  0,    0
             +k13,      0,   -b31-k34,  0,    0
                0,   +k24,       +k34,  0,    0
             +k12,      0,       +k34,  0, -kBu]
    dy = k*y
endfunction

function f = myDifferences ( rate, m ) 
    global count;
    count = count + 1;
    // Returns the difference between the simulated differential 
    // equation and the experimental data.
    y_calc=ode(C0',t0,t_exp,list(myModel,rate)) 
    diffmat = y_calc' - y_exp
    // Make a column vector
    f = diffmat(:)
endfunction 

function [rate,SSD,diffopt]=mySolve(rate_start,rate_fix)
    rate0 = []
    rate_opt_i = []
    for i=1:rate_n
        if rate_fix(i)==0 then
           rate_opt_i = [rate_opt_i,i] 
           rate0 = [rate0,rate_start(i)]
//           rate_opt_i($+1) = i 
//           rate0($+1) = rate_start(i)
        end
    end
    //disp(rate_start)
    //disp(rate_opt_i)
    //disp(rate0)
    if length(rate_opt_i)==0 then
        rate1_nonfixed = []
        diffopt = myDifferences(rate0,m)
    else
        [rate1_nonfixed,diffopt]=lsqrsolve(rate0,myDifferences,m)
    end
    //disp(rate1_nonfixed)
    SSD = sum(diffopt.^2)
    rate = rate_start
    for i=1:length(rate_opt_i)
//        if rate1_nonfixed(i)<0 then
//            rate1_nonfixed(i) = 0
//        end
        rate(rate_opt_i(i)) = rate1_nonfixed(i)
    end
    //disp(rate)
    
endfunction

[rate1,SSD,diffopt] = mySolve(rate_start,rate_fix)

[my,ny] = size(y_exp)
if isdef("exp_err") then
//if isdef("exp_err") & exp_err ~= zeros(y_exp) then
    chi_square = sum(matrix(diffopt, my, ny).^2 ./ exp_err.^2)
else
   chi_square = N-M
end
//disp(chi_square)

STD = sqrt(SSD/(N-M))

for i=1:length(rate1)
    if rate_fix(i) == 0 then
        fff='  *' 
    else
        fff=''
    end
    mprintf("%s\t%e%s\n",rate_name(i),rate1(i),fff)
end
mprintf("Irerations: %d\nSSD (sum of squares): %e\nStandard deviation: %e\n", count, SSD, STD);

t_plot = t_exp(1):(t_exp($)-t_exp(1))/256:t_exp($);
rate_opt_i = 1:rate_n
y_plot = ode(C0',t0,t_plot,list(myModel,rate1));
    


f0=scf(0);
clf(f0)
plot(t_plot',y_plot')
if isdef("HBr_exp") then
    plot(t_plot',C_HBr(t_plot)','m')
end
plot(t_exp,y_exp,'o')
if isdef("HBr_exp") then
    plot(t_exp_HBr,HBr_exp,'m*')
end
xtitle('','Time, s', 'Mole fraction')
legend(reagent_name, HBr_name, -1)
//legend('I','II','III','IV','HBr',-1)

///////////////////////////////////////////////////////////////////////////////
if isdef("hessian") & rate_n_opt>=1 & hessian>0 then
    mprintf("\n%s\n", 'Use of Hessian matrix to obtain the standard deviations of the fitting parameters')
    mprintf("nsigma=%d\n\n",nsigma)
    
    if rate_n_opt<rate_n then
        mprintf("Warning! There is fixed perameters. Bad results are expected.\n\n")
    end
    
    function ssd=return_ssd(rate)
        y_calc=ode(C0',t0,t_exp,list(myModel,rate)) 
        diffmat = (y_calc' - y_exp)/STD
        ssd = sum(diffmat.^2)
    endfunction
    //ssd_ttt = return_ssd(rate1)
    
    [J, H] = numderivative(return_ssd, rate1, [], [], "blockmat")
    //disp(H)
    H_inv = inv(H)
    rate_std = sqrt(2*diag(H_inv)) * nsigma
    for i=1:rate_n
        if rate_fix(i) ~= 0 then
            continue
        end
        mprintf("%s\t%e +/- %.1e (%.1f%%)\n",...
          rate_name(i),rate1(i),rate_std(i),rate_std(i)/rate1(i)*100)
    end

end
///////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////
if isdef("variation") & rate_n_opt>=1 & variation>0 then
    mprintf("\n%s\n", 'Use of constant chi-square boundaries as confidence limits')
    mprintf("nsigma=%d\n\n",nsigma)
    //target_SSD = (nsigma^2 + 1) * SSD  // doubling of SSD
    // delta_chi_square = 1 where chi_square = SSD/(SSD/(N-M))
    target_SSD = SSD + SSD/(N-M)*nsigma^2
    var_step = 0.001

    f1=scf(1);
    clf(f1);
    
    pos = 1
    better_solve = []
    prev_ssd = SSD
    for ind=1:rate_n
        if rate_fix(ind)~=0 then
            continue
        end
        var_mat = [SSD,rate1]
        rate_fix_var = rate_fix
        rate_fix_var(ind) = 1
        ssd = 0
        count_var = 1
        rate_var = rate1
        too_big_variation = %f
        
        while ssd < target_SSD
            if rate_var(ind) > rate1(ind)*2 then
                too_big_variation = %t
                break
            end
            rate_var(ind) = rate_var(ind)*(1+var_step*count_var)
            //disp(1+step*count_var)
            [rate_var,ssd]=mySolve(rate_var,rate_fix_var)
            var_mat = cat(1, var_mat,[ssd,rate_var])
            mprintf("%8.3f%8.3f\r", rate_var(ind)/rate1(ind), (ssd/SSD-1)*(N-M))
            if ssd<prev_ssd then
                better_solve = rate_var
                prev_ssd = ssd
            end
            //mprintf("+%.3f\r",ssd/SSD)
            //mprintf("%6d\r",count_var)
            count_var = count_var + 1
        end
        ssd = 0
        count_var = 1
        rate_var = rate1
        while ssd < target_SSD
            if too_big_variation then
                break
            end
            if rate_var(ind) < rate1(ind)/2 then
                break
            end
            rate_var(ind) = rate_var(ind)*(1-var_step*count_var)
            if rate_var(ind) < rate1(ind)*var_step then
                too_big_variation = %t
                break
            end
            //disp(1+step*count_var)
            [rate_var,ssd]=mySolve(rate_var,rate_fix_var)
            var_mat = cat(1, [ssd,rate_var], var_mat)
            mprintf("%8.3f%8.3f\r", rate_var(ind)/rate1(ind), (ssd/SSD-1)*(N-M))
            if ssd<prev_ssd then
                better_solve = rate_var
                prev_ssd = ssd
            end
            //mprintf("-%.3f\r",ssd/SSD)
            count_var = count_var + 1
        end
        //disp(var_mat)
        
        if too_big_variation then
            mprintf("%s\t%e +/- >100%%\n",rate_name(ind),rate1(ind))
        else
            rate_minus = var_mat(1,2:rate_n+1) - ...
              (var_mat(1,2:rate_n+1)-var_mat(2,2:rate_n+1)) * ...
              (var_mat(1,1)-target_SSD)/(var_mat(1,1)-var_mat(2,1))
            //disp(rate_minus)
            rate_plus = var_mat($,2:rate_n+1) - ...
              (var_mat($,2:rate_n+1)-var_mat($-1,2:rate_n+1)) * ...
              (var_mat($,1)-target_SSD)/(var_mat($,1)-var_mat($-1,1))
            //disp(rate_plus)
            rate_err = ((rate_plus-rate1)+(rate1-rate_minus))/2
            mprintf("%s\t%e +/- %.1e (%.1f%%)\n",...
              rate_name(ind),rate1(ind),rate_err(ind),rate_err(ind)/rate1(ind)*100)
            //disp(var_mat)
        end
        
        rrr_var_mat = var_mat(:,2:rate_n+1)
        for i=1:rate_n
            if rate_fix(i)==1 then
               rrr_var_mat(:,i) = ones(rrr_var_mat(:,i))
            elseif rate1(i)==0 then
               continue
            else
               rrr_var_mat(:,i) = rrr_var_mat(:,i)/rate1(i)
            end
        end
        //disp(rrr_var_mat)
        rrr_min = min(rrr_var_mat)
        rrr_max = max(rrr_var_mat)
        
        subplot(1,rate_n_opt,pos)
        plot(rrr_var_mat(:,ind),rrr_var_mat)
        dc=gca();
        dc.axes_visible =["off", "off"]
        dc.data_bounds=[rrr_min rrr_max rrr_min rrr_max]
        xtitle(rate_name(ind))
        pos = pos + 1
    end
//    subplot(1,rate_n_opt+1,pos)
    legend(rate_name)
    if ~isempty(better_solve) then
        mprintf("%s\n", "Warning! There is a better solution:")
        for i=1:length(better_solve)
            mprintf("%e ", better_solve(i))
        end
        mprintf("%s\n\n","Please, try it as a starting approximation...")
    end
    
end
/////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////
if isdef("montecarlo") & rate_n_opt>=1 & montecarlo>0 then
    //mark_size=int(1+6/log(montecarlo))
    if montecarlo>300 then
        mark_size=2
    elseif montecarlo>100 then
        mark_size=3
    elseif montecarlo>30 then
        mark_size=4
    else
        mark_size=6
    end
    mprintf("\n%s\n", 'Confidence limits by Monte Carlo simulation (montecarlo method)');
    mprintf("%d synthetic data sets, nsigma = %d\n\n", montecarlo, nsigma);
    [my,ny] = size(y_exp)
    resid = matrix(diffopt, my, ny)
    
    mean_res = zeros(1,ny)
    //mean_res = mean(resid,'r')
    
    std_res = stdev(resid) //Use total chi-square for all experimental curves
    //std_res = stdev(resid,'r') //Use own chi-square for each experimental curve
    
    if length(std_res)==1 then  // if total chi-square
       mean_res = zeros(1,ny)
       std_res = ones(1,ny) * std_res
    end
    
    std_res = std_res * nsigma
    rate_bs = rate1
    
//    disp(rate_name)
//    disp(rate1')
    for i=1:montecarlo
        resid1 = []
        for j=1:ny
          resid1 = cat(2,resid1,grand(my, 1, "nor", mean_res(j), std_res(j)))
        end
        y_exp = Y_EXP + resid1
        [rrr,ssd,diffopt]=mySolve(rate1,rate_fix)
        //[rrr,diffopt]=lsqrsolve(rate1,myDifferences,m)
        //rate_bs($+1,:) = rrr'
        rate_bs = cat(1,rate_bs,rrr)
//        disp(rate_new')
        mprintf("\r%4d %8.3f", i, sqrt(ssd/SSD-1))
    end
    mprintf("\r%s",'')
    rate_std = stdev(rate_bs, 'r');
    //rate_av = rate1  //Use rates optimized from experimental data
    rate_av = mean(rate_bs, 'r')  //Use averaged rates from synthetic data sets
    for i=1:rate_n
        if rate_fix(i) ~= 0 then
            continue
        end
        mprintf("%s\t%e +/- %.1e (%.1f%%)\n",...
          rate_name(i),rate_av(i),rate_std(i),rate_std(i)/rate_av(i)*100)
    end
    
    // rrr_ - data of optimized parameters only (for plot)
    rrr_name = []
    rrr_bs = []
    rrr_n = rate_n_opt
    for i=1:rate_n
        if rate_fix(i)==0 then
            rrr_name = [rrr_name,rate_name(i)]
            rrr_bs = [rrr_bs,rate_bs(:,i)]
        end
    end
    
    if rrr_n > 1 then
        f2=scf(2);
        clf(f2);
        for i=2:rrr_n
            for j=1:i-1
                pos =((i-2)*(rrr_n-1)+j)
                //disp(pos)
                subplot(rrr_n-1,rrr_n-1,pos)
                R = correl(rrr_bs(:,j),rrr_bs(:,i))
                plot(rrr_bs(:,j),rrr_bs(:,i),'.')
                dc=gca();
                dc.axes_visible =["off", "off"]
                a=get("current_axes");
                p1=a.children.children(1);
                //set(p1,'mark_mode',"on");
                set(p1,'mark_size',mark_size);
                xtitle(rrr_name(j)+'-'+rrr_name(i)+'  correl '+msprintf("%.3f",R))
            end
        end
    end
end
///////////////////////////////////////////////////////////////////////////////