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

%This function is the engine of the fitting routine, based on globfun_bio.
%However, this script does very little, because it is only called when all 
%parameters are fixed. A quick glance of the 'fit' to your data can be
%obtained.

function res = globfun_bio_fix(start_pars)

%Defined in globfit_bio_gui
global xdata ydata Alg_switch
global start_pars_lin
%Defined in fit_data
global comps obj
global Par_fix_switch Par_fix_spectra_matrix
global Par_fix_IRF IRF
%Defined here
global decays spectra rev_counter
global Par_IRF IRF_fun y_fit

%All models; Setup start values concentration profile
CSOBJ = getconfigset(obj, 'active');
set(CSOBJ, 'StopTime', max(xdata));

%Count number of reversible reactions
%Important for 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
for i=1:size(obj.Reactions,1)               %Set start values to generate concentration profiles
    if obj.Reactions(i).Reversible==0
        obj.Reactions(i).KineticLaw.Parameters.value=start_pars(i+rev_counter);
    elseif obj.Reactions(i).Reversible==1
        obj.Reactions(i).KineticLaw.Parameters(1).value=start_pars(i+rev_counter);
        obj.Reactions(i).KineticLaw.Parameters(2).value=start_pars(i+1+rev_counter);
        rev_counter=rev_counter+1;
    end
end
clear i

%Create concentration profiles for user-determined timpoints (i.e. measured) 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;
tobj_fit.time=tobj_fit.time+start_pars(comps+1);
tobj_fit.time(end)=tobj.time(end);
decays_fit=resample(tobj_fit,xdata,'zoh'); %Interpolate timesdecays_start.Data now contains resampled data
decays=decays_fit.Data;
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 are generated
Par_IRF=0; %If IRF function is generated, calculate also its spectrum
if (Par_fix_IRF==1 && length(IRF)==2) %If IRF parameters are present
    FWHM_gauss_tilde=start_pars(comps+2)/(2*sqrt(2*log(2)));
    for i=1:comps
        IRF_comps(:,i)=0.5*exp(start_pars(i)*(start_pars(comps+3)+0.5*start_pars(i)*FWHM_gauss_tilde^2))*...
            (1+erf((xdata-(start_pars(comps+3)+start_pars(i)*FWHM_gauss_tilde^2))/sqrt(2)/FWHM_gauss_tilde));
        decays_IRF(:,i)=decays(:,i).*IRF_comps(:,i);
    end
    clear i
    IRF_fun=(1/(FWHM_gauss_tilde*sqrt(2*pi)))*exp(-log(2)*(2*(xdata-start_pars(comps+3))/start_pars(comps+2)).^2);
    Par_IRF=1;
    decays=cat(2,decays_IRF,IRF_fun);
    clear decays_IRF
end


%Fitting of linear parameters by Matlab built-in function of mlr (matrix
%left division, mldivide), which is very fast.
y_fit_temp=ydata+start_pars_lin;
spectra=decays(:,1:comps+Par_IRF-rev_counter)\y_fit_temp;   %Last decay curve not fitted, as its spectrum is zero!
clear y_fit_temp

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

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

if Alg_switch==1
    res=y_fit-ydata;
elseif Alg_switch==2
    res=sum(sum((y_fit-ydata).^2));
end

clear y_fit