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

function globfit_bio_gui
%This script is called by fit_data. It uses variables from that script,
%and performs data manipulation, singular value decomposition if required,
%calls the fitting function, and sets its (algorithm options).
%The fittin results are displayed in the Command window, and plotting
%functions are called.

%Defined in fit_data
global comps obj
global analysis
global nr_svd_comps SVDcontrib_switch start_par_time_zero
global Alg_choice Alg_accuracy IRF start_pars_nonlin start_pars_lin start_pars start_pars_dispers
global Par_fix_switch Par_fix_time Par_fix_time_zero Par_fix_OD_offset Par_fix_IRF Par_fix_spectra_matrix
global time_shift wave_shift time_limits wave_limits
global time_divider wave_divider sub_timepoint sub_wavelength
global plot_fit_prog_live
global Error_switch_nonlin Error_switch_lin_spectra

%Defined here 
%Needed for glob_fun_bio/data_props)
global Alg_switch weight CSOBJ Par_fix_time time_fix_counter Par_fix_time_mat
%Needed for plot_figs (y_fit also needed for fit_explorer_plot)
global nr_waves y_fit res RESNORM RESIDUAL S V opt_pars opt_specs opt_times sing_val
%Needed for save_settings
global ci_other ci_times_errorbars 
%Needed for save_figs
global OUTPUT

%Defined in fit_data
global xdata ydata wave plot_residuals
%Defined in glob_fun_bio
global spectra decays rev_counter
global opt_pars_lin opt_pars_lin_std Par_IRF
%Defined in plot_figs
global svd_comp_perc
%Defined in dispersion_IRF_calc
global delay_shift_curve
%Defined in disp_cor_data
global ydata_cor opt_disp_pars 
%%%%%%%%%%%%%%%%

%Data manipulation
xdata=xdata+time_shift; %Shift time axis with specified amount
wave=wave+wave_shift;       %Shift wavelength axis with specified amount

%Decrease the amount of data to be analysed, by user-set limits or by increasing stepsize
if length(wave_limits)==2 && length(time_limits)==1
    ydata=ydata(:,wave_limits(1):wave_limits(2));
    wave=wave(wave_limits(1):wave_limits(2));
    weight=zeros(size(ydata));
elseif length(wave_limits)==1 && length(time_limits)==2
    xdata=xdata(time_limits(1):time_limits(2));
    ydata=ydata(time_limits(1):time_limits(2),:);
    weight=zeros(size(ydata));
elseif length(wave_limits)==2 && length(time_limits)==2
    xdata=xdata(time_limits(1):time_limits(2));
    ydata=ydata(time_limits(1):time_limits(2),wave_limits(1):wave_limits(2));
    weight=zeros(size(ydata));
    wave=wave(wave_limits(1):wave_limits(2));
end
clear wave_shift time_shift;

if wave_divider>1;
    ydata=ydata(:,1:wave_divider:end);
    wave=wave(1:wave_divider:end);
    weight=zeros(size(ydata));
end
if time_divider>1;
    xdata=xdata(1:time_divider:end,1);
    ydata=ydata(1:time_divider:end,:);
    weight=zeros(size(ydata));
end

%Data weighting
%Works by multiplying the selected range of the residual with the desired
%weight. It cycles over multiple rows for more weights.
time_weight=str2num(get(findobj('Tag','time_weight_box'),'String'));
if time_weight~=0
    if size(time_weight,2)~=3
        errordlg('To give a time range a weight, give weight, then range (from, to), like ''0.1 1 20''). If multiple time weights are needed, seperate them with a semicolon.')
        error('To give a time range a weight, give weight, then range (from, to), like ''0.1 1 20''). If multiple time weights are needed, seperate them with a semicolon.')
    end
    if isempty(find(time_weight(:,2:end)==0))==0
            errordlg('Cannot find time index zero. Check ''Time weight''.')
            error('Cannot find time index zero. Check ''Time weight''.')
    end
    if numel(find(abs(mod( time_weight(:,2:end),round(time_weight(:,2:end))))==0)) ~= 2*size(time_weight,1)
        errordlg('Time weight range should be in indices.')
        error('Time weight range should be in indices.')
    end 
end
wave_weight=str2num(get(findobj('Tag','wave_weight_box'),'String'));
if wave_weight~=0
    if size(wave_weight,2)~=3
        errordlg('To give a time range a weight, give weight, then range (from, to), like ''0.1 1 20''). If multiple time weights are needed, seperate them with a semicolon.')
        error('To give a time range a weight, give weight, then range (from, to), like ''0.1 1 20''). If multiple time weights are needed, seperate them with a semicolon.')
    end
    if isempty(find(wave_weight(:,2:end)==0))==0
            errordlg('Cannot find wavelength index zero. Check ''Time weight''.')
            error('Cannot find wavelength index zero. Check ''Time weight''.')
    end
    if numel(find(abs(mod( wave_weight(:,2:end),round(wave_weight(:,2:end))))==0)) ~= 2*size(wave_weight,1)
        errordlg('Wavelength weight range should be in indices.')
        error('Wavelength weight range should be in indices.')
    end 
end

if size(time_weight,2)==3  || size(wave_weight,2)==3
    weight=ones(size(ydata));
    if max(time_weight)>size(weight,1)
        errordlg('Check that time range selected for weighting is smaller than the number of time points to be analysed.')
        error('Check that time range selected for weighting is smaller than the number of time points to be analysed.')
    end 
    if max(wave_weight)>size(weight,2)
        errordlg('Check that wavelength range selected for weighting is smaller than the number of wavelengths to be analysed.')
        error('Check that wavelength range selected for weighting is smaller than the number of wavelengths to be analysed.')
    end 
    if size(time_weight,2)==3
        for w_t=1:size(time_weight,1)
            weight(time_weight(w_t,2):time_weight(w_t,3),:)=time_weight(w_t,1);
        end
        clear w_t
    end
    
    if size(wave_weight,2)==3
        for w_w=1:size(wave_weight,1)
            weight(:,wave_weight(w_w,2):wave_weight(w_w,3))=wave_weight(w_w,1);
        end
        clear w_w
    end
elseif time_weight==0 && wave_weight==0
    weight=1;
end

%Parameter initialisation
nr_waves=length(wave);          %Number of wavelengths
nr_times=size(xdata,1);         %Number of time points
y_fit=zeros(size(ydata));
res=zeros(size(ydata));

%Count the number of fixed rate constants
Par_fix_time_mat=zeros(comps+1,1); %Create matrix containing only zeros if all rate constant are free, but with number of rate constant in the same index if it's fixed
if isequal(Par_fix_time,0)==1
    time_fix_counter=0;
elseif isequal(Par_fix_time,0)==0
    time_fix_counter=length(Par_fix_time);
    for i=1:comps+1
        for j=1:time_fix_counter
            if isequal(i,Par_fix_time(j))==1
                Par_fix_time_mat(i,1)=Par_fix_time(j);
            end
        end
    end
    clear i
end

%Set model parameters
%Define start values concentration profile
CSOBJ = getconfigset(obj, 'active');
set(CSOBJ, 'StopTime', max(xdata));
set(CSOBJ.CompileOptions, 'DimensionalAnalysis', false);

%Baseline correction(s): Or whole spectrum at selected time point (or average
%of range of time points), or time trace at selected wavelength (or average
%of range of wavelengths
if length(sub_timepoint)==1 && sub_timepoint>0
    ydata_sub_timepoint=repmat(0,nr_times,nr_waves);    %Parameter initilisation
    for i=1:size(ydata,1)
        ydata_sub_timepoint(i,:)=ydata(i,:)-ydata(sub_timepoint,:);
    end
    clear i
    ydata=ydata_sub_timepoint;
elseif length(sub_timepoint)==2 && sub_timepoint(1)>0 && sub_timepoint(2)>0
    ydata_sub_timepoint=repmat(0,nr_times,nr_waves);    %Parameter initilisation
    ydata_sub_timepoint_ave=sum(ydata(sub_timepoint(1):sub_timepoint(1),:),1)/(abs(sub_timepoint(2)-sub_timepoint(1))+1);
    for i=1:size(ydata,1)
        ydata_sub_timepoint(i,:)=ydata(i,:)-ydata_sub_timepoint_ave;
    end
    clear i
    ydata=ydata_sub_timepoint;
end
clear ydata_sub_timepoint ydata_sub_timepoint_ave
if length(sub_wavelength)==1 && sub_wavelength>0
    ydata_sub_wavelength=repmat(0,nr_times,nr_waves);    %Parameter initilisation
    for j=1:size(ydata,2)
        ydata_sub_wavelength(:,j)=ydata(:,j)-ydata(:,sub_wavelength);
    end
    clear j
    ydata=ydata_sub_wavelength;
elseif length(sub_wavelength)==2 && sub_wavelength(1)>0 && sub_wavelength(2)>0
    ydata_sub_wavelength=repmat(0,nr_times,nr_waves);    %Parameter initilisation
    ydata_sub_wavelength_ave=sum(ydata(:,sub_wavelength(1):sub_wavelength(2)),2)/(abs(sub_wavelength(2)-sub_wavelength(1))+1);
    for j=1:size(ydata,2)
        ydata_sub_wavelength(:,j)=ydata(:,j)-ydata_sub_wavelength_ave;
    end
    clear j
    ydata=ydata_sub_wavelength;
end
clear ydata_sub_wavelength ydata_sub_wavelength_ave

%%
%Using only positive times
%ADD CONSTANT TO 'XDATA' IF PART OF NEGATIVE TIMES NEED TO BE FITTED AS WELL
pos_test=find(xdata<=0);
if isempty(pos_test)==0
    xdata=xdata(max(pos_test)+1:end);
    ydata=ydata(max(pos_test)+1:end,:);
    res=zeros(size(ydata));
    if isequal(weight,1)==0 %If no weight is applied, then next line makes no sense
        weight=weight(max(pos_test)+1:end,:);
    end
end
clear pos_test;

%%%%%%%%%%%%%%%
%For SVD analysis only requested number of SVD components is fitted
%Preparing data for svd analyis.
%If desired (determined by 'SVDcontrib_switch'), an additional calculation
%can be done, meaning that are calculatd:
%left singular values (LSV)*singular values (S)
%This can be used as measure of the relative time component contribution to
%whole dataset
%If required, weights can be applied as well to full data set
if nr_svd_comps>0
    [U,S,V]=svd(ydata.*weight);
    if SVDcontrib_switch==0
        ydata=U(:,1:nr_svd_comps);
    elseif SVDcontrib_switch==1
        clear US;
        for i=1:nr_svd_comps
            US(:,i)=U(:,i)*S(i,i);
        end
        clear i
        ydata=US;
    end
    sing_val=diag(S);
    nr_waves=size(ydata,2);
end
clear wave_divider time_divider;

%#################################################
%%
%Fitting routine
%#################
disp('Starting optimisation, please wait...');
tic;    %Start recording elapsed time of fit

if isempty(start_pars_nonlin)==0 %At least one parameter must be passed to fitting routine. If all parameters are fixed, use different routine.
      
%     if isequal(Par_fix_time,0)==1   %Determine how many rate constants need to be fitted
%         time_fix_counter=0;         %Number that represents the amount of fixed rate constants
%     elseif isequal(Par_fix_time,0)==0
%         time_fix_counter=length(Par_fix_time);
%     end
    LB(1:comps-time_fix_counter)=0; %Set lower bound for rate constants zero, so that they can only be positive
    UB=[];
    if length(IRF)==2 && Par_fix_IRF==0
        LB(comps+2-time_fix_counter-Par_fix_time_zero)=0; %FWHM of IRF must be greater than zero
    end

    if Alg_choice==0
        Alg_switch=1; %Algorithm used: lsqnonlin
        options=optimset('Display','iter'); %Setting the algorithm options to display in Command Window
        if Alg_accuracy==0
            options.TolFun=1e-6;
        elseif Alg_accuracy==1
            options.TolFun=1e-12;
        end
        if plot_fit_prog_live==0
            options.PlotFcns=[];
        elseif plot_fit_prog_live==1
            options.PlotFcns={@plot_fit_prog_live_fun_lsqnonlin,@optimplotresnorm};
        end
        
        if Error_switch_nonlin==0 %Different output options if errors in non-linear parameters need to be calculated or not
            [opt_pars,RESNORM,RESIDUAL,~,OUTPUT]= lsqnonlin(@globfun_bio, start_pars_nonlin, LB,UB,options);
            %Jacobian not needed, so not generated
        elseif Error_switch_nonlin==1
            [opt_pars,RESNORM,RESIDUAL,~,OUTPUT,~,jacobian_fit]= lsqnonlin(@globfun_bio, start_pars_nonlin, LB,UB,options);
            jacobian_fit=full(jacobian_fit); %Makes error estimation MUCH faster (due to smaller matrix size and method) with identical results (f.i. 0.1 vs 210 sec.!)
        end
    elseif Alg_choice==1
        Alg_switch=2; %Algorithm used: patternsearch
        options = psoptimset('Display','iter','CompletePoll','on');
        
        if Alg_accuracy==0
            options.TolFun=1e-6;
        elseif Alg_accuracy==1
            options.TolFun=1e-12;
        end
        if plot_fit_prog_live==0
            options.PlotFcns=[];
        elseif plot_fit_prog_live==1
            options.PlotFcns={@plot_fit_prog_live_fun_pattern,@psplotbestf};
        end
        
        [opt_pars,RESNORM,~,OUTPUT] = patternsearch(@globfun_bio,start_pars_nonlin,[],[],[],[],LB,UB,options);
        %Note that different 'res' must be defined in object function: res=sum(sum((y_fit-ydata).^2));;
        %Note also that the 'Vectorized' option is apparently much faster, but
        %this is not implemented
        
        %Generating RESIDUAL as it is not generated automatically as for
        %algorithm 1 (lsqnonlin). y_fit is generated by globfun_bio
        RESIDUAL=y_fit-ydata;
    end
elseif isempty(start_pars_nonlin)==1
    if Error_switch_nonlin==1 || Error_switch_lin_spectra==1
        errordlg('All parameters are fixed, so errors cannot be calculated. Uncheck boxes in Error estimation panel (both rate constants/time zero and spectra)!')
        error('All parameters are fixed, so errors cannot be calculated. Uncheck boxes in Error estimation panel (both rate constants/time zero and spectra)!')
    end
    %Generate 'fit' with all parameters fixed
    RESIDUAL=globfun_bio_fix(start_pars);
    opt_pars=start_pars;
    OUTPUT.iterations=0;
    time_fix_counter=comps;
end

disp('Optimisation completed');
toc;

%#################
%%
%Preparing output data for visualization
%#################################################
%Generating time constants (opt_times) from fitted decay parameters, adding 
%the fixed ones
fix_counter=0; %Counts fixed parameters as rate constants are set up
if isequal(Par_fix_time,0)==1
    opt_times=1./opt_pars(1:comps);
elseif isequal(Par_fix_time,0)==0
   if isequal(time_fix_counter,comps)==1 %If all rate constants are fixed
        opt_times=1./start_pars(1:comps);
    elseif isequal(time_fix_counter,comps)==0
        opt_times=zeros(1,comps);           %If one or more rate constants are fixed
        for i=1:comps
            for j=1:time_fix_counter
                if isequal(i,Par_fix_time(j))==1
                    opt_times(1,i)=1/start_pars(i); %First fill in the fixed ones
                    fix_counter=fix_counter+1;
                end
            end
        end
        free_pars_find=find(opt_times==0);
        if isempty(free_pars_find)==0
            for k=1:length(free_pars_find)
                opt_times(free_pars_find(k))=1/opt_pars(k); %Then fill in the fitted ones
            end
        end
    end
    clear i j k
end
clear fix_counter

if nr_svd_comps==0
    disp('Global fit of all data');
elseif nr_svd_comps~=0
    disp(['Global fit of ',num2str(nr_svd_comps),' SVD components'])
end
if size(Par_fix_spectra_matrix,2)==2 && Par_fix_switch==1
    for s=1:size(Par_fix_spectra_matrix,1)
        disp(['Spectrum of state ',num2str(Par_fix_spectra_matrix(s,1)),' is equal to ',num2str(Par_fix_spectra_matrix(s,2))])
    end
    clear s
end

%Type of reaction rate law used is displayed. 'k*species' corresponds to 
%mass-action kinetics
clear rate_law_display
for i=1:size(obj.reactions,1)
    rate_law_display(i,:)=obj.reactions(i).ReactionRate;
end

if analysis==1
    disp(' Reaction rates used:')
    disp(rate_law_display);
    disp('Optimised parameters PARALLEL model')
elseif analysis==2
    disp(' Reaction rates used:')
    disp(rate_law_display);
    disp('Optimised parameters SEQUENTIAL model')
elseif analysis==3
    obj.reactions   %Used target model
    disp(' Reaction rates used:')
    disp(rate_law_display);
    obj.species     %Input matrix
    disp('Optimised parameters TARGET model')
end

if Par_fix_OD_offset==0
    disp([        '      OD offset:   ',num2str(opt_pars_lin)])
elseif Par_fix_OD_offset==1
    disp([        '      OD offset:   ',num2str(start_pars_lin),' (fixed)'])
end
if Par_fix_time_zero==0
    disp([        '      time zero:   ',num2str(opt_pars(comps+1-time_fix_counter))])
elseif Par_fix_time_zero==1
    disp([        '      time zero:   ',num2str(start_pars(comps+1)),' (fixed)'])
end

if sum(Par_fix_time)==0
    disp([        ' Time constants:   ',num2str(opt_times)])
elseif  sum(Par_fix_time)>0 && length(Par_fix_time)<comps
    disp([        ' Time constants:   ',num2str(opt_times),' (Index fixed component(s): ',num2str(Par_fix_time),')'])
elseif length(Par_fix_time)==comps && isequal(Par_fix_time,0)==0
    disp([        ' Time constants:   ',num2str(opt_times),' (All fixed)'])
end
if length(IRF)==2
    if Par_fix_IRF==0
        disp([    'IRF (FWHM, t0):   ',num2str(opt_pars(comps+2-time_fix_counter-Par_fix_time_zero:comps+3-time_fix_counter-Par_fix_time_zero))]);
    elseif Par_fix_IRF==1
        disp([    'IRF (FWHM, t0):   ',num2str(start_pars(comps+2:comps+3)),' (fixed)']);
    end
end

%Generating confidence intervals of coefficients
if Error_switch_nonlin == 1
    if Par_fix_switch==1
        disp('Errors are generated, but there are fixed parameters')
    end
    if isempty(jacobian_fit)==1 %jacobian_fit is not generated if user interrupts fit, therefore ci is initiated so that program finishes normally
        ci=(opt_pars'*ones(1,2));
        disp('Errors not generated due to user-aborted fit.');
    elseif isempty(jacobian_fit)==0
        ci = nlparci(opt_pars,RESIDUAL,'jacobian',jacobian_fit);
        disp('Errors calculated');
    end
elseif Error_switch_nonlin == 0
    if Par_fix_IRF==1 || IRF(1)==0
        ci=(opt_pars(1:comps+1-Par_fix_time_zero-time_fix_counter)'*ones(1,2));
    elseif Par_fix_IRF==0 && length(IRF)==2
        if comps>time_fix_counter
            ci=(opt_pars(1:comps+2-Par_fix_time_zero-time_fix_counter)'*ones(1,2));
        elseif comps==time_fix_counter
            if Par_IRF==0
                ci=zeros(size(opt_pars));
            elseif Par_IRF==1
                ci=zeros(length(start_pars),2);
            end
        end
    end
end
ci_times=ci(1:comps-time_fix_counter,:);
if Par_fix_time_zero==0
    ci_other=ci(size(ci_times,1)+1,:); %Confidence intervals of time zero
elseif Par_fix_time_zero==1
    ci_other=[0 0];
end
if Par_fix_IRF==0 && length(IRF)==2
    ci_other=cat(1,ci_other,ci(end-1:end,:));
elseif Par_fix_IRF==1
    ci_other=cat(1,ci_other,zeros(2,2));
end

%Generating confidence interval for optimised times; Errors propagate
%because decay constants are inversed to obtain decay times
%Only for fitted rate constants: d(tau) = (d(k)/k) * tau, with k the rate
%constant, which is the inverse of the time constant tau. 
if comps>time_fix_counter
    for c=1:comps-time_fix_counter
        ci_times_inv(c,1)=ci_times(c,2)-opt_pars(c);
        ci_times_inv(c,2)=opt_pars(c)-ci_times(c,1);
    end
    ci_times_errorbars=(ci_times_inv(:,1)./opt_pars(1:comps-time_fix_counter)').*opt_times(1:comps-time_fix_counter)';
elseif comps-time_fix_counter==0 %If all rate constants are fixed
    ci_times_inv=zeros(comps,2);
    ci_times_errorbars=zeros(comps,1);
end

if eq(time_fix_counter,comps)==0 && isequal(Par_fix_time,0)==0
    ci_times_errorbars_temp=zeros(comps+Par_fix_time,1); %Create temporary error matrix for all time constants (all zeros)
    for k=1:length(free_pars_find)
            ci_times_errorbars_temp(free_pars_find(k))=ci_times_errorbars(k); %Then fill in the fitted ones
    end
    ci_times_errorbars=ci_times_errorbars_temp;
end
clear ci_times ci_times_inv c free_pars_find ci_times_errorbars_temp

%Generating spectra from the optimised parameter matrix opt_pars
disp('Preparing plots');

opt_specs=spectra';         %Fitted spectra
%%
%On-screen display of fit results
%#################
% disp(OUTPUT);

if Error_switch_nonlin==1
    if Par_fix_switch==1
        disp(     'Note that there are fixed parameters (i.e. can be other(s) than the rate constants)')
        disp(    ['   Error (+/-STD): ',num2str(ci_times_errorbars(1:comps)')])
        if Par_fix_time_zero==1 && Par_fix_OD_offset==0
            disp(['    Error (OD,to): ',num2str(opt_pars_lin_std),' 0'])
        elseif Par_fix_time_zero==0 && Par_fix_OD_offset==0
            disp(['    Error (OD,t0): ',num2str(opt_pars_lin_std),' ',num2str(ci_other(1,2)'-opt_pars(comps+1-time_fix_counter))])
        elseif Par_fix_time_zero==0 && Par_fix_OD_offset==1
            disp(['    Error (OD,t0): 0 ',num2str(ci_other(1,2)'-opt_pars(comps+1-time_fix_counter))])
        elseif Par_fix_time_zero==1 && Par_fix_OD_offset==1
            disp( '    Error (OD,t0): 0 0')
        end
        if length(IRF)==2 && Par_fix_IRF==0
            disp(['Err.IRF(FWHM, t0):  ',num2str(ci_other(2:3,2)'-opt_pars(comps+2-time_fix_counter-Par_fix_time_zero:comps+3-time_fix_counter-Par_fix_time_zero))])
        elseif length(IRF)==2 && Par_fix_IRF==1
            disp(['Er.IRF(FWHM, t0):  0 0'])
        end
    elseif Par_fix_switch==0
        disp([    '   Error (+/-STD): ',num2str(ci_times_errorbars')])
        disp([    '    Error (OD,t0): ',num2str(opt_pars_lin_std),' ',num2str(ci_other(1,2)'-opt_pars(comps+1-time_fix_counter))])
        if length(IRF)==2 && Par_fix_IRF==0
            disp(['Er.IRF(FWHM, t0):  ',num2str(ci_other(2:3,2)'-opt_pars(comps+2-time_fix_counter-Par_fix_time_zero:comps+3-time_fix_counter-Par_fix_time_zero))])
        end
    end
elseif Error_switch_nonlin==0
    disp(         '   Error (+/-STD): Not calculated')
end
disp(' ')
disp([            '         Residual: ',num2str(RESNORM)])
disp([            '       Iterations: ',num2str(OUTPUT.iterations)])
%%
%#################
%Plotting
plot_figs;
%%
%#################
if nr_svd_comps>0 %Needs 'svd_comp_perc' calculated by plot_figs
    disp(['First ',num2str(nr_svd_comps),' singular values : ',num2str(sing_val(1:nr_svd_comps)')]);
    if SVDcontrib_switch==1
        disp(['Contribution per species to data (%): ',num2str(svd_comp_perc)]);
    end
end
%#################
if OUTPUT.iterations==0
    msgbox('Iterations = 0!','Warning')     %This only produces a box, but does not stop execution
    disp('Iterations = 0!')
end