%Copyright (c) 2010, Luuk J.G.W. van Wilderen
%
%globfun_bio.m

function res = globfun_bio(start_pars_nonlin)
%This function is the engine of the fitting routine.
%The decay curves are generated here, by making use of lots of Matlab's
%built-in routines.
%The amount of parameters to be fitted (start_pars_nonlin) varies, as rate
%constants can be fixed, so they will not be called. This makes the coding
%a bit more complicated. Check globfun_bio_fix for simple version and
%obtain understanding of what is done (that function is called when all 
%parameters are fixed).
%The OD baseline is fitted with polyfit (using a straight line on the 
%residuals).
%The spectra are fitted by mldivide (or lscov), which use QR decomposition
%if the data matrix is not square. This should give a unique solution.
%The traces are fitted by feeding the changed parameters into the model,
%and letting sbiosimulate solve the differential equations that determine
%the concentration profiles.

%Defined in globfit_bio_gui
global xdata ydata y_fit res time_fix_counter
global Alg_switch start_pars start_pars_lin CSOBJ 
%Defined in fit_data
global comps obj nr_svd_comps
global Par_fix_switch Par_fix_spectra_matrix Error_switch_lin_spectra Par_fix_time_zero
global Par_fix_time Par_fix_OD_offset Par_fix_IRF IRF Alg_interpol_method Par_fix_time_mat
global weight dispersion_correction
%Defined in dispersion_IRF_calc or in load_dispersion if loaded
global delay_shift_curve
%Defined here
global decays spectra spectra_std rev_counter y_fit
global opt_pars_lin_std Par_IRF IRF_fun

%Set rate constants for each species. Importantly, it also counts the number 
%of reversible reactions, which may occur in target models with reversible 
%transitions. Without this section the used starting parameters may not 
%correspond to its identity, e.g. time zero could wrongfully be used as 
%rate constant.
rev_counter=0; %Counts amount of reversible reactions
fix_counter=0; %Counts fixed parameters as rate constants are set up
for i=1:size(obj.Reactions,1)               %Set start values to generate concentration profiles
    if obj.Reactions(i).Reversible==0
        if isequal(i+rev_counter,Par_fix_time_mat(i))==0   %Amount of rate constants to fit varies as they may be fixed
            obj.Reactions(i).KineticLaw.Parameters.value=start_pars_nonlin(i+rev_counter-fix_counter);
        elseif isequal(i+rev_counter,Par_fix_time_mat(i))==1
            obj.Reactions(i).KineticLaw.Parameters.value=start_pars(i+rev_counter);
            fix_counter=fix_counter+1;
        end
    elseif obj.Reactions(i).Reversible==1
        if isequal(i+rev_counter,Par_fix_time_mat(i))==0
            obj.Reactions(i).KineticLaw.Parameters(1).value=start_pars_nonlin(i+rev_counter-fix_counter);
        elseif isequal(i+rev_counter,Par_fix_time_mat(i))==1
            obj.Reactions(i).KineticLaw.Parameters(1).value=start_pars(i+rev_counter);
            fix_counter=fix_counter+1;
        end
        if isequal(i+1+rev_counter,Par_fix_time_mat(i))==0
            obj.Reactions(i).KineticLaw.Parameters(2).value=start_pars_nonlin(i+1+rev_counter-fix_counter);
        elseif isequal(i+1+rev_counter,Par_fix_time_mat(i))==1
            obj.Reactions(i).KineticLaw.Parameters(2).value=start_pars(i+1+rev_counter);
            fix_counter=fix_counter+1;
        end
        rev_counter=rev_counter+1;
    end
end
clear i fix_counter

%Create concentration profiles for measured timpoints and specified model
[T_fit,Conc_fit] = sbiosimulate(obj,CSOBJ);
tobj_fit = timeseries(Conc_fit,T_fit); %Create time series object only to use built-in interpolation function
tobj=tobj_fit;
if Par_fix_time_zero==1 %If time zero is fixed, do not fit it
    tobj_fit.time=tobj_fit.time+start_pars(comps+1);
elseif Par_fix_time_zero==0
        tobj_fit.time=tobj_fit.time+start_pars_nonlin(comps+1-time_fix_counter); %Added time zero
end
tobj_fit.time(end)=tobj.time(end);

%Time points generated by sbiosimulate do not match xdata => interpolate
%timesdecays_start. Data now contains resampled data == measured data
%Although the linear option is the preferred choice (fast and accurate), it
%cannot extrapolate. This may give problems when time zero is fitted (or
%IRF). Therefore, manually switch to any of the other options, and save the
%optimised values for a new 'linear' run.
if Alg_interpol_method==1 %linear==default
    decays_fit=resample(tobj_fit,xdata,'linear');
    decays=decays_fit.Data;
elseif Alg_interpol_method==2
    decays_fit=resample(tobj_fit,xdata,'zoh'); %Time points generated by sbiosimulate do not match xdata => interpolate timesdecays_start. Data now contains resampled data == measured data
    decays=decays_fit.Data;
elseif Alg_interpol_method==3
    decays = interp1(tobj_fit.time,tobj_fit.data,xdata,'cubic'); 
elseif Alg_interpol_method==4
    decays = interp1(tobj_fit.time,tobj_fit.data,xdata,'spline');
end
clear CSOBJ decays_fit tobj tobj_fit T_fit Conc_fit tobj_fit;

%Instrument response function convolution with data (according to BBA van Stokkum 2004
%1657 82-104).
%Dispersion function could be included, as in van Global and target analysis of time-
%resolved spectra Lecture notes Troisième Cycle de la Physique en Suisse Romande, March
%14-24, 2005 Ivo H.M. van Stokkum
%IRF corrected concentration profiles (decays) and IRF profile only are
%generated.
%If IRF function is generated, calculate also its spectrum. This parameter
%appears a lot, as its concentration profile and spectrum are added as last 
%column (but only used when IRF is needed).
Par_IRF=0; 
if length(IRF)==2
    [decays,IRF_fun]=glob_fun_IRF(start_pars,start_pars_nonlin,Par_fix_time_mat,decays,time_fix_counter);
    Par_IRF=1;
    decays=cat(2,decays,IRF_fun);
end

%Fitting of linear parameters by Matlab built-in function of mlr (matrix
%left division, mldivide), which is very fast.
%If the errors need to be calculated, lscov is used. Both mldivide and lscov
%use QR decomposition, and give identical results. However 'lscov' is
%slightly slower, so therefore this is only used when the errors are
%requested. It is also possible that mldivide convergences to a slightly
%less accurate solution, but these difference are minimal (of the order of 1e-4).
%NB. Last decay curve not fitted, as its spectrum is zero! This is because
%the before-last species needs to decay into another species (without
%amplitude).
[spectra]=globfun_spectra(decays,start_pars_lin,rev_counter);

%Spectral constraints; Equalise pairs of spectra.
%It cycles over multiple rows for more constraints.
if size(Par_fix_spectra_matrix,2)==2 && Par_fix_switch==1              
    for s=1:size(Par_fix_spectra_matrix,1)                             
        spectra(Par_fix_spectra_matrix(s,2),:)=spectra(Par_fix_spectra_matrix(s,1),:);
    end
    clear s
end

%Data reconstruction from fitted parameters
y_fit=decays(:,1:comps+Par_IRF-rev_counter)*spectra;  

%Dispersion correction applied to reconstructed data (y_fit)
%For every channel, each time point is shifted with a specific delay time
%as determined by the dispersion curve (can be loaded or generated).
if dispersion_correction==1
    for i=1:size(ydata,2)
        %Fast interpolator that does no input checking is used because time 
        %points are most likely irregularly spaced.
        yi = interp1q(xdata,y_fit(:,i),xdata-delay_shift_curve(i));
        y_fit(:,i)=yi;
    end
    y_fit(isnan(y_fit)) = 0; %Remove any NaN's ('interp1q' above cannot extrapolate)
end

%Generate the residual value (res) to be minimised by the function.
if Alg_switch==1
    if nr_svd_comps==0
        res=(y_fit-ydata).*weight;
    elseif nr_svd_comps>0
        res=(y_fit-ydata); %Weight already applied to ydata
    end
elseif Alg_switch==2
    if nr_svd_comps==0
        res=sum(sum(((y_fit-ydata).*weight).^2));
    elseif nr_svd_comps>0
        res=sum(sum((y_fit-ydata).^2));  %Weight already applied to ydata
    end
end