%% Tutorial 4. Fitting of NMR titration series
% by Evgenii Kovrigin, 05/17/2011
%
% In this tutorial we will explore fitting of the series of NMR line shapes
% with different models. We will synthesize data using a complex
% multi-state model and then fit data with a simpler
% one as example (we do not doubt the correct model will fit data---what
% about the wrong one?).
%%


%% 
% Clean up workspace
% (to avoid unwanted effects of old variable assignments)
close all
clear all


%% CHOOSE PARAMETERS OF THE TITRATION SERIES
% NOTE: L/R cannot be ZERO. Set some very small number instead!
Rtotal = 1e-3;
LRratio_array = [1e-6 0.2 0.4 0.6 0.8 1.0 1.2];
Y_RMSD=1e-7;  % ideal line shape
Y_offset_percent=0;  % to offset figures vertically for easier viewing


%% CREATE A SERIES OF DATASETS FOR SIMULATION OF LINE SHAPES IN TITRATION
% We will use a complex model U-R:
% Transition A: R+L<=>RL
% Transition B: R<=>R*


session_name='Tutorial_4';
session=TotalFit(session_name);
%% 
% Create a simulation series
series_name='test_series'; % arbitrary
dataset_class_name='NMRLineShapes1D'; % class name
model_name='U_R-model'; 
model_numeric_solver='analytical';
model_numeric_options='no options';
number_of_datasets=length(LRratio_array);
number_of_points_to_simulate=30;
X_min=-200;
X_max=800;
X_column=linspace(X_min, X_max, number_of_points_to_simulate)'; % use ' to transpose to a column-vector

session.create_simulation_series(series_name, dataset_class_name, model_name, ...
                 model_numeric_solver, model_numeric_options,...
                 number_of_datasets, X_column);
             
session.show_series_datasets(series_name)  % show only datasets of the series
session.initialize_relation_matrix(sprintf('%s.all_datasets.txt', session_name))
session.generate_fitting_environment(sprintf('%s.fitting_environment.txt', session_name)) 

%% CREATE PARAMETER LINKAGES
parameter_number_array=[ 3 4 5 6 7 8 9 10 11 12 13 ];
session.batch_link_parameters_in_series(series_name, parameter_number_array);

%%
% Propagate information about links across all datasets
session.generate_fitting_environment(sprintf('%s.fitting_environment.txt', session_name)) 


%% ASSIGN PARAMETERS 
% Choose some reasonable parameters (taken from Bill Hawse's spectra of MHC
% residue Q72)

setup.log10_K_A = log10(1e7) ; % Assume tight binding
setup.log10_K_B = log10(3)   ; % Assume significant population shift towards non-binding species
setup.k_2_A = 10 ; % slow off rate 
setup.k_2_B = 500 ; % fast isomerization
setup.w0_R = 0 ; % resonance frequencies in units 1/s (Hz*2pi)
setup.w0_Rstar=100; 
setup.w0_RL = 500; 
setup.FWHH_R = 150; % /s (Hz*2pi)
setup.FWHH_Rstar = 150 ; 
setup.FWHH_RL = 190 ; % complex formation increases molecular weight
setup.ScaleFactor = 1; % (no meaning in simulation context)

%%
% prepare vectors for easier handling:

FAKE_VALUE=1E-6;  % individual (unlinked) parameters: to be assigned later

ideal_parameters=[ FAKE_VALUE    FAKE_VALUE     setup.log10_K_A    setup.log10_K_B ...
    setup.k_2_A setup.k_2_B ...
    setup.w0_R  setup.w0_Rstar  setup.w0_RL...
    setup.FWHH_R setup.FWHH_Rstar setup.FWHH_RL...
    setup.ScaleFactor ] 

fake_limits=zeros(1,length(ideal_parameters));

%% 
% *Assign linked parameters*
% (I keep names of function arguments for information)
first_dataset_index=session.dataset_index(series_name,1);

Monte_Carlo_range_min=fake_limits;
Monte_Carlo_range_max=fake_limits;
parameter_limits_min=fake_limits;
parameter_limits_max=fake_limits;

session.assign_parameter_values(first_dataset_index, ideal_parameters, ...
        Monte_Carlo_range_min, Monte_Carlo_range_max, parameter_limits_min, parameter_limits_max);

%% 
% *Assign unlinked parameters*
Rtotal_vector=ones(1,length(LRratio_array))*Rtotal; % this one is constant for all datasets
%%
% Create vectors of ranges for formal assignment (NOTE: they are different length from the above!)
unlinked_fake_ranges=zeros(1,number_of_datasets);

MonteCarlo_range_min_array=unlinked_fake_ranges;
MonteCarlo_range_max_array=unlinked_fake_ranges;
parameter_limits_min_array=unlinked_fake_ranges;
parameter_limits_max_array=unlinked_fake_ranges;


%%
% Assign Rtotal
parameter_number=1;
values_array=Rtotal_vector;

session.assign_unlinked_parameter_in_series(...
             series_name, parameter_number, values_array, ...
             MonteCarlo_range_min_array, MonteCarlo_range_max_array, ...
             parameter_limits_min_array, parameter_limits_max_array);

%%
% Assign LRratio
parameter_number=2;
values_array=LRratio_array;

session.assign_unlinked_parameter_in_series(...
             series_name, parameter_number, values_array, ...
             MonteCarlo_range_min_array, MonteCarlo_range_max_array, ...
             parameter_limits_min_array, parameter_limits_max_array);

         
         
%% CALCULATE MODEL VALUES
% Run a simulation
X_standard_dev=0; % no need to have errors of Y
session.simulate_series(series_name, X_standard_dev, Y_RMSD);


%%
% Plot the simulation result
figure_handle1=session.plot_1D_series(series_name, 'Simulation 4', Y_offset_percent);

%%
% output simulation
figure_handle=figure_handle1;
folder_name=session_name;
figure_name=series_name;
results_output.output_figure(figure_handle, folder_name, figure_name)






%% PREPARE FOR FITTING WITH SIMPLE U MODEL
% Save original simulation session
session_simulation=session.make_a_copy('Tutorial_4_simulation');
%%
% Set new session name
session_name='Tutorial_4_U_fit';
session.set_name(session_name);

%%
% Switch model
new_model_name='U-model';
model_numeric_solver='analytical';
model_numeric_options='no options';

session.change_model_in_series(series_name, new_model_name,...
                model_numeric_solver, model_numeric_options)
session.initialize_relation_matrix(sprintf('%s.all_datasets.txt', session_name))
session.generate_fitting_environment(sprintf('%s.fitting_environment.txt', session_name)) 

            
%% CREATE PARAMETER LINKAGES
parameter_number_array=[ 3 4 5 6 7 8 9 ];
session.batch_link_parameters_in_series(series_name, parameter_number_array);

%%
% Propagate information about links across all datasets
session.generate_fitting_environment(sprintf('%s.fitting_environment.txt', session_name)) 
            
            

%% ASSIGN STARTING PARAMETERS 
% 

%%
% Binding constant
Log10_K_A=log10(1e6);
Log10_K_A_AbsMin=log10(1e1);
Log10_K_A_AbsMax=log10(1e10);


%%
% Dissociation rate constant
k_2_A=100;
k_2_A_AbsMin=0.1;
k_2_A_AbsMax=3000;

%%
% Frequency of R
w0_R=50;
w0_R_AbsMin=X_min; % set these limits to the spectral window
w0_R_AbsMax=X_max;


%%
% Frequency of RL
w0_RL=500;
w0_RL_AbsMin=X_min; % set these limits to the spectral window
w0_RL_AbsMax=X_max;

%%
% Line width of R
FWHH_R=100;
FWHH_R_AbsMin=10;  
FWHH_R_AbsMax=500;



%%
% Line width of RL
FWHH_RL=100;
FWHH_RL_AbsMin=10;  
FWHH_RL_AbsMax=500;


%%
% ScaleFactor
ScaleFactor=1;
ScaleFactor_AbsMin=0.1;  
ScaleFactor_AbsMax=10;



%%
% Create a vector of parameter values
starting_parameters=[ FAKE_VALUE     FAKE_VALUE   ...
                        Log10_K_A    k_2_A  ...
                       w0_R     w0_RL      FWHH_R     FWHH_RL    ScaleFactor ];
                 
%%
% Create vectors of absolute limits. Parameters, which will not be
% optimized do not need meaningful limits (use 0).
min_limits=[ 0 0 Log10_K_A_AbsMin k_2_A_AbsMin w0_R_AbsMin w0_RL_AbsMin  FWHH_R_AbsMin  FWHH_RL_AbsMin ScaleFactor_AbsMin ];
max_limits=[ 0 0 Log10_K_A_AbsMax k_2_A_AbsMax w0_R_AbsMax w0_RL_AbsMax  FWHH_R_AbsMax  FWHH_RL_AbsMax ScaleFactor_AbsMax ];


%% 
% *Assign linked parameters*
% (I keep names of function arguments for information)
% I am formally setting Monte-Carlo range to the same arrays as absolute
% limits. In reality, they will not be used because we are using the default mode
% vary_starting_parameters='no' meaning that starting parameters are not
% varied for new Monte-Carlo runs.

first_dataset_index=session.dataset_index(series_name,1);

Monte_Carlo_range_min=min_limits;
Monte_Carlo_range_max=max_limits;
parameter_limits_min=min_limits;
parameter_limits_max=max_limits;

session.assign_parameter_values(first_dataset_index, starting_parameters, ...
        Monte_Carlo_range_min, Monte_Carlo_range_max, parameter_limits_min, parameter_limits_max);



%% 
% *Assign unlinked parameters*
% Create vectors of ranges for formal assignment (different length from the above!)
unlinked_fake_ranges=zeros(1,number_of_datasets);
%%
% Assign Rtotal --- this will not be fit so no need for ranges
parameter_number=1;
values_array=Rtotal_vector;  % this remains the same
MonteCarlo_range_min_array=unlinked_fake_ranges;
MonteCarlo_range_max_array=unlinked_fake_ranges;
parameter_limits_min_array=unlinked_fake_ranges;
parameter_limits_max_array=unlinked_fake_ranges;

session.assign_unlinked_parameter_in_series(...
             series_name, parameter_number, values_array, ...
             MonteCarlo_range_min_array, MonteCarlo_range_max_array, ...
             parameter_limits_min_array, parameter_limits_max_array);
    
    
%%
% Assign LRratio
parameter_number=2;
values_array=LRratio_array;

session.assign_unlinked_parameter_in_series(...
             series_name, parameter_number, values_array, ...
             MonteCarlo_range_min_array, MonteCarlo_range_max_array, ...
             parameter_limits_min_array, parameter_limits_max_array);

    
%%                          
% Inspect initial guess
figure_handle1=session.plot_1D_series(series_name, 'U_fitting', Y_offset_percent);
%%
% Good terrible guess.

%%
% Save initial state of fitting session on a disk for future inspection and reuse
session.export_results(sprintf('%s.initial_guess', session_name), {figure_handle1})

         
         
         
         
%% FIT

%%
% Set which parameters of the model will be fixed and variable
parameter_modes=[ 0 0 1 1 1 1 1 1 0];
session.set_parameters_modes(first_dataset_index, parameter_modes);


%%
% Select fitting algorithm. Choices: 'simplex', 'MultiStart'
% 'GlobalSearch', 'DirectSearch'. Do not forget to start your parallel
% processors:
% *|matlabpool open|*
fit_only_algorithm='GlobalSearch';


%%
% Report on start-up and show variable parameter names
fprintf(1,'\n\n==============\nFitting Algorithm: %s\nVariable parameters:', fit_only_algorithm);
session.fit_parameter_names{1:end}

%%
% Fit without error determination
tic
session.fit_only(fit_only_algorithm);
toc

%% 
% Report on finishing
session.show_fit_statistics()
session.transfer_all_parameters_to_individual_datasets(); 

%%
% Inspect final result
figure_handle2=session.plot_1D_series(series_name, 'U_final', Y_offset_percent);

%%
% Export final result
session.export_results(sprintf('%s.Fit_Only.%s', session_name, fit_only_algorithm), {figure_handle2})

         
%% CONCLUSION
% In slow exchange, two-state model easily can be fit to the three-state
% data: we are not sensitive to the model structure enough.
         
         
         
         
 