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

function fit_data
%Called by Globe_toolbox_v1_00.
%If 'FIT' button is pressed, values from GUI are passed to actual fitting
%routine (called by globfit_bio_gui). It also performs error checking
%and gives tips for proper use of the program.

global comps obj 
global filename time_choice model ydata xdata dispersfilename delay_shift_curve
%PANEL 2: Model selection
global analysis
%PANEL 3: Data selection
global wave nr_svd_comps SVDcontrib_switch calibfilename
%PANEL 4: Fitting parameters
global start_par_OD_offset start_par_time_zero Alg_choice Alg_accuracy Alg_interpol_method
global Par_fix_switch Par_fix_time Par_fix_time_zero Par_fix_OD_offset 
global IRF start_pars_nonlin Par_fix_IRF Par_fix_spectra_matrix Par_fix_disp
global start_pars start_pars_lin dispersion_correction
%PANEL 5: Data manipulation
global time_shift wave_shift time_limits wave_limits
global time_divider wave_divider sub_timepoint sub_wavelength
%PANEL 6: Plotting options
global axis_switch plot_residuals plot_fit_prog_live lin_log 
global wave_digits plot_wave_divider
%PANEL 7: Error estimation
global Error_switch_nonlin Error_switch_lin_spectra

%####################
if exist('filename')==0 || isempty(filename)==1
    errordlg('Load dataset first!')     %This only produces a box, but does not stop execution
    error('Load dataset first!')        %This does stop the execution
end

time_choice=str2num(get(findobj('Tag','time_choice_box'),'String'));   %Get rate constants from user input
if isempty(time_choice)==1
    errordlg('Define rate constants first!')
    error('Define rate constants first!')
end
comps=length(time_choice);      %Number of components (decay constants)

if get(findobj('Tag','model1'),'Value')==1
    analysis=1;
    model=cat(2,'par',num2str(comps),'.xml');
elseif get(findobj('Tag','model2'),'Value')==1
    analysis=2;
    model=cat(2,'seq',num2str(comps),'.xml');
elseif get(findobj('Tag','model3'),'Value')==1
    analysis=3;
    if get(findobj('Tag','model3_fix_toggle'),'Value')==0
        [model,model_path]=uigetfile('*.xml','Give model file');
        cd (model_path)
    end
end
set(findobj('Tag','load_model_text'),'Visible','on','String',['       Model file: ' model]);
obj = sbmlimport(model);

%If a model is used that is not supplied with the software, a special 
%consideration needs to be taken. 
%Firstly, construct your own model with each consecutive component having a 
%name that starts with a new letter of the alphabet.
%Secondly, define your rate constants as k1, k2, etc, in increasing name 
%(start from 1 to 9) and decreasing order (rate constants have units 
%1/(time units), so k1=1 1/s, k2=0.1 1/s for instance).
%Thirdly, the last defined component has zero spectral amplitude, i.e. a 
%model of two states A and B, linked by one rate constant, will result in 
%plots showing only the (spectral) amplitude of the first species A only, 
%since species B has no amplitude.
%Below the species are organised on alphabet.
species_alphabet='ABCDEFGHIJKLMNOPQRSTUVWXYZ';
species_alphabet=species_alphabet(1:length(obj.species));
for i=1:length(obj.species)
    for j=1:length(obj.species)
        test_spec=find(obj.Compartments.species(i).Name(1)==species_alphabet(j));
        if test_spec==1
            order_species(i)=j;
        end
    end
end
[~,order_species_ind]=sort(order_species);
reorder(obj.Compartments,obj.Compartments.species(order_species_ind'));
%Because the reorder command shows inconsistent behavior, the rate
%constants need to be arranged separately according to (increasing) number 
%k#, with # up to 99:
for k=1:length(obj.reactions)
    for m=1:1:9
        test_rates=strcmp(obj.reactions(k).reactionRate(2:3),[num2str(m) '*']);
        if test_rates==1
            order_reaction_rates(k)=m;
        end
    end
    for m=10:1:99
        test_rates=strcmp(obj.reactions(k).reactionRate(2:3),num2str(m));
        if test_rates==1
            order_reaction_rates(k)=m;
        end
    end
end
reorder(obj,obj.reactions(order_reaction_rates'));
clear order_species_ind test_spec species_alphabet order_reaction_rates test_rates

if analysis==3
    rev_counter=0;
    for i=1:size(obj.Reactions,1)
        if obj.Reactions(i).Reversible==1
            rev_counter=rev_counter+1;
        end
    end
    if length(time_choice)~=size(obj.Reactions,1)+rev_counter
        errordlg('Number of rate constants not compatible with model file!')
        error('Number of rate constants not compatible with model file!')
    end
end

data=load (filename);
ydata=data(:,2:end);
xdata=data(:,1);

%Only load dispersion curve file if name exists and is not empty
%Don't do anything if it is not required. The 'Apply' dispersion toggle
%button decides if dispersion file is applied to data or not.
if exist('dispersfilename')==1 && isempty(dispersfilename)==0
    delay_shift_curve=load (dispersfilename);
    if length(delay_shift_curve)~=size(ydata,2)
        errordlg('Size dispersion curve data file does not match size of datafile!');
        error('Size dispersion curve data file does not match size of datafile!');
    end
    
end

%PANEL 2: Model selection
%NB. 'analysis' already defined in load_model

%PANEL 3: Data selection
if get(findobj('String','SVD'),'Value')==1
    nr_svd_comps=str2double(get(findobj('Tag','nr_SVD_comps_box'),'string'));
    SVDcontrib_switch= get(findobj('Tag','SVD_contrib_toggle'),'Value');
elseif get(findobj('String','All data'),'Value')==1
    nr_svd_comps=0;
    SVDcontrib_switch=0;
end
if get(findobj('String','SVD'),'Value')==1 && str2double(get(findobj('Tag','nr_SVD_comps_box'),'string'))==0
    errordlg('Set number of SVD components > 0!')
    error('Set number of SVD components > 0!');
end
if exist('calibfilename')==1 && isempty(calibfilename)==0
    wave=load (calibfilename);
    if length(wave)~=size(data,2)-1
        errordlg('Calibration file does not match datafile!');
        error('Calibration file does not match datafile!');
    end
    if size(wave,2)~=length(wave)
        wave=wave';
    end
elseif exist('calibfilename')==0 || isempty(calibfilename)==1
    wave=1:size(ydata,2);
end

%PANEL 4: Fitting parameters
start_par_OD_offset=str2double(get(findobj('Tag','OD_offset_box'),'String'));
start_par_time_zero=str2double(get(findobj('Tag','time_zero_box'),'String'));
Alg_choice=get(findobj('Tag','Alg_choice_toggle'),'Value');
Alg_accuracy=get(findobj('Tag','Alg_accuracy_toggle'),'Value');
Alg_interpol_method=str2num(get(findobj('Tag','Alg_interpol_box'),'String'));
dispersion_correction=get(findobj('Tag','dispersion_correction_toggle'),'Value');

if Alg_interpol_method ~= 1 && Alg_interpol_method ~= 2 && Alg_interpol_method ~= 3 && Alg_interpol_method ~= 4
    errordlg('Interpolation method can only be integer from 1 to 4.')
    error('Interpolation method can only be integer from 1 to 4.');
end

if str2num(get(findobj('Tag','IRF_box'),'String'))==0
    IRF=0;
elseif length(str2num(get(findobj('Tag','IRF_box'),'String')))==2
    IRF=str2num(get(findobj('Tag','IRF_box'),'String'));
end
if exist('delay_shift_curve','var')==0 || isempty(delay_shift_curve)==1 && dispersion_correction==1
    errordlg({'Dispersion file does not exist!';'Uncheck ''Apply'' button in dispersion parameters panel.'})
    error('Dispersion file does not exist! Uncheck ''Apply'' button in dispersion parameters panel.')
end
if exist('delay_shift_curve','var')==1 && length(delay_shift_curve)~=length(wave) && dispersion_correction==1
    errordlg('Size of dispersion correction does not match size of data file!')
    error('Size of dispersion correction does not match size of data file!')
end

%Start parameters for non-linear fit.
%NB Spectra and OD offset are actually fitted linearly (which is much faster).
%Dispersion correction is fitted separately
if IRF(1)==0
    start_pars_nonlin=cat(2,time_choice,start_par_time_zero); %decay times, tzero
elseif length(IRF)==2
    start_pars_nonlin=cat(2,time_choice,start_par_time_zero,IRF); %decay times, tzero, FWHM and mu0 of IRF
end
start_pars_lin=start_par_OD_offset;

%Constraints
Par_fix_switch=get(findobj('Tag','par_fix_switch_toggle'),'Value');
if Par_fix_switch==0
    Par_fix_time=0;
    Par_fix_time_zero=0;
    Par_fix_OD_offset=0;
    Par_fix_IRF=0;
    Par_fix_disp=0;
    Par_fix_spectra_matrix=0;
    start_pars=start_pars_nonlin;
elseif Par_fix_switch==1
    Par_fix_time=str2num(get(findobj('Tag','Par_fix_time_box'),'String'));
    Par_fix_time_zero=get(findobj('Tag','Par_fix_time_zero_toggle'),'Value');
    Par_fix_OD_offset=get(findobj('Tag','Par_fix_OD_offset_toggle'),'Value');
    Par_fix_IRF=get(findobj('Tag','Par_fix_IRF_toggle'),'Value');
    start_pars=start_pars_nonlin;
    %Strip fixed parameters from end if fixed
    if Par_fix_IRF==1 && length(IRF)==2 %Do not fit IRF parameter if they are fixed
        start_pars_nonlin(comps+2:comps+3)=[];
    end
    if Par_fix_time_zero==1 %Do not fit time-zero parameter if it is fixed
        start_pars_nonlin(comps+1)=[];
    end
    if Par_fix_time~=0
        for i=length(Par_fix_time):-1:1
            start_pars_nonlin(i)=[];
        end
    end
     
    start_pars_lin=start_par_OD_offset; %Linear parameter
    Par_fix_spectra_matrix=str2num(get(findobj('Tag','Par_fix_spectra_matrix_box'),'String'));
    if Par_fix_spectra_matrix~=0
        if size(Par_fix_spectra_matrix,2)~=2
            errordlg('To fix two spectra, give them in pairs. If multiple pairs need to be fixed, seperate the pairs with a semicolon.')
            error('To fix two spectra, give them in pairs. If multiple pairs need to be fixed, seperate the pairs with a semicolon.');
        end
        if isempty(find(Par_fix_spectra_matrix==0))==0
            errordlg('Cannot find spectrum zero. Check ''Fix spectral pairs'' box.')
            error('Cannot find spectrum zero. Check ''Fix spectral pairs'' box.')
        end
        if numel(find(abs(mod( Par_fix_spectra_matrix,round(Par_fix_spectra_matrix)))==0)) ~=2*size(Par_fix_spectra_matrix,1)
            errordlg('Spectral component not found, should be index.  Check ''Fix spectral pairs'' box.')
            error('Spectral component not found, should be index.  Check ''Fix spectral pairs'' box.')
        end
    end
end

if IRF(1)==0 && Par_fix_switch==1 && Par_fix_IRF ==1
    errordlg('FWHM of IRF is set to 0 with the fixed IRF constraints activated. Set FWHM of IRF to value > 0, or release IRF values in constraints.')
    error('FWHM of IRF is set to 0 with the fixed IRF constraints activated. Set FWHM of IRF to value > 0, or release IRF values in constraints.');
end


%*****PANEL 5: Data manipulation*****
time_shift=str2double(get(findobj('Tag','time_shift_box'),'String'));
wave_shift=str2double(get(findobj('Tag','wave_shift_box'),'String'));
if str2double(get(findobj('Tag','time_lim_box1'),'String'))==0 && str2double(get(findobj('Tag','time_lim_box2'),'String'))==0
    time_limits=0;
elseif str2double(get(findobj('Tag','time_lim_box1'),'String'))~=0 && str2double(get(findobj('Tag','time_lim_box2'),'String'))~=0
    time_limits(1)=str2double(get(findobj('Tag','time_lim_box1'),'String'));
    time_limits(2)=str2double(get(findobj('Tag','time_lim_box2'),'String'));
end
if str2double(get(findobj('Tag','wave_lim_box1'),'String'))==0 && str2double(get(findobj('Tag','wave_lim_box2'),'String'))==0
    wave_limits=0;
elseif str2double(get(findobj('Tag','wave_lim_box1'),'String'))~=0 && str2double(get(findobj('Tag','wave_lim_box2'),'String'))~=0
    wave_limits(1)=str2double(get(findobj('Tag','wave_lim_box1'),'String'));
    wave_limits(2)=str2double(get(findobj('Tag','wave_lim_box2'),'String'));
end
time_divider=str2double(get(findobj('Tag','time_divider_box'),'String'));
wave_divider=str2double(get(findobj('Tag','wave_divider_box'),'String'));
sub_timepoint=str2num(get(findobj('Tag','time_subtract_box'),'String'));
sub_wavelength=str2num(get(findobj('Tag','wave_subtract_box'),'String'));

if nr_svd_comps>0 && (time_divider>1 || wave_divider>1)
    SVD_error=errordlg('Are you sure you want to calculate SVD for less than the whole dataset? Otherwise, set Time/Wavelength divider to 0.');
    uiwait(SVD_error);
    uiresume;
end

%*****PANEL 6: Plotting options*****
axis_switch=get(findobj('Tag','axis_switch_toggle'),'Value');
plot_residuals=get(findobj('Tag','plot_residuals_toggle'),'Value');
plot_fit_prog_live=get(findobj('Tag','plot_fit_prog_live_toggle'),'Value');
lin_log=get(findobj('Tag','lin_log_toggle'),'Value');
wave_digits=str2num(get(findobj('Tag','wave_digits_box'),'String'));
plot_wave_divider=str2num(get(findobj('Tag','plot_wave_divider_box'),'String'));

%*****PANEL 7: Error estimation*****
Error_switch_nonlin=get(findobj('Tag','Error_switch_nonlin_toggle'),'Value');
Error_switch_lin_spectra=get(findobj('Tag','Error_switch_lin_spectra_toggle'),'Value');

%Checks
%Because the 'resample' function (only used for 'timeseries objects') is
%used, the time zero offset parameter can only be negative or zero. The
%reason is that interpolation is not possible for negative time points.
% if start_par_time_zero>0
%     errordlg('Time offset parameter must be =<0. Instead, use ''Time shift'' to compensate for positive time zero.');
%     error('Time offset parameter must be =<0. Instead, use ''Time shift'' to compensate for positive time zero.');
% end

if IRF(1)<0
    errordlg('IRF has zero or negative width. Set FWHM>0.');
    error('IRF has zero or negative width. Set FWHM>0.');
elseif IRF(1)/2 >= 1/max(start_pars(1:comps)) && Par_fix_IRF==0
    errordlg('To be able to fit IRF, FWHM/2 must be shorter than fastest rate constant.');
    error('To be able to fit IRF, FWHM/2 must be shorter than fastest rate constant.');
end

%Error estimation not possible for 'patternsearch'
if Alg_choice==1 && Error_switch_nonlin==1
    errordlg('Error estimation for decay rates not possible when selecting Pattern Search!')
    error('Error estimation for decay rates not possible when selecting Pattern Search!')
end

%Fitting routine called
globfit_bio_gui;