function Simulate(setup_name, model_name)
%------------------- LineshapeKin Simulation -----------------------------
version= '4.1.5' ;
%
% 5/8/09
% Written by Evgenii Kovrigin, Medical College of Wisconsin
%
% 
%   This function runs simulation of a titration series 
%
%  NOTE:
%   for models with R2, R2L2 the population fraction is a fraction of R
%   nuclei in these species so p(R2)=2*[R2]/Rtotal
% 
% 5/10/09 - K-matrix is tested for consistency. However, manual checking
%                 must be done by the user because not all errors may be  detected
%                 automatically.
% 5/21/09 - Spectral data is saved in LineShapeKin 3.1 format for fitting
%
% 5/27/09 - determine peak linewidths and export both the figure and the data for comparison plotting
%
% 5/30/09 - a population graph bug fixed
% 6/10/09 - all models converted to a new input format
% 9/28/09 - output for BiophysicsLab added
% 10/4/10 - output of plain-text format for ITC curves added
% 5/26/2011 - calculation of K-matrix for U-R2 and U-R2L2 models
%                    corrected. The B_R2_R2L2 model was temporarily disabled (needs the same
%                    correction to be done).
% 2/25/2012 - Add checking for the spectral range to include frequencies of species
%                    


close all;



global YES NO;

YES=1;
NO=0;


%% -- Get name of the model to simulate
fprintf(1,'\n\nLineShapeKin Simulation ver. %s.\nEvgenii Kovrigin, Medical College of Wisconsin\n3/24/09\n', version);
if nargin==2
    fprintf(1,'Simulating model: %s\n using general settings from: %s\n', model_name, setup_name);
 else
    fprintf(1,'\nUsage: Simulate   general_setup_name  model_setup_name\n\n');
    return
end



%%  LOAD SETTINGS

% -- Get general settings ---
[conditions, modes, Model]=load_setup(setup_name, model_name);
Model.name=model_name;



%% Check for errors in the input data
if abs(length(conditions.Rtotal)-conditions.total_titration_points)>0
    error('Lengths of Rtotal[] and PLratio[] are different!!!');
end

if min(Model.w0)<conditions.w_min || max(Model.w0)>conditions.w_max
    error('ERROR!    Check your setup file: the spectral window [w_min, w_max] must be sufficient to include frequencies of all species!')
end




%% Make paths for results 
if ~exist(conditions.results_path,'dir')
    mkdir(conditions.results_path);
end

conditions.figures_path=fullfile(conditions.results_path,Model.name);  
if ~exist(conditions.figures_path,'dir')
     mkdir(conditions.figures_path);
end

conditions.lineshapekin31_path=fullfile(conditions.results_path,Model.name,'LineShapeKin_3');
if ~exist(conditions.lineshapekin31_path,'dir')
     mkdir(conditions.lineshapekin31_path);
end

conditions.BiophysicsLab_path=fullfile(conditions.results_path,Model.name,'BiophysicsLab');
if ~exist(conditions.BiophysicsLab_path,'dir')
     mkdir(conditions.BiophysicsLab_path);
end




%%  Begin calculatiions
% - start timer
tic;
fprintf(1,'Project:  %s\n%s \n', Model.name, Model.description);

%%  Compute equilibrium populations
% compute smooth curves
LRpoints=conditions.LRpoints;
modes.extend_factor=1.2; % how far to show beyond last LRratio point
LRrange=0: modes.extend_factor*max(conditions.LRratio)/LRpoints : modes.extend_factor*max(conditions.LRratio);
smooth_p=[];

ITC_curve=zeros(LRpoints,1); % compute relative ITC curve

% test if Rtotal was fixed or variable and give warning of inadequacy of smooth curves with variable Rtotal
for  LR_counter=1:conditions.total_titration_points
    if conditions.Rtotal(1)~=conditions.Rtotal(LR_counter)
        fprintf('\nWARNING: Rtotal is changing - the smooth curves calculated in Populations graph and ITC will be INACCURATE !!!\n');
        break
    end
end


for LR_counter=1:LRpoints+1
    curr_LR=LRrange(LR_counter);
    
    %% Compute equilibrium concentrations
    % I assume Rtotal is not changning for smooth curve
    % calculations. Test that and give warning
    [p, K, R2_vector, v_vector, p_other, Model]=model_matrices(Model, conditions.Rtotal(1), curr_LR);
    smooth_p= [ smooth_p;  p p_other];
    
    
    %% Test consistency of the K-matrix: (1) by mass conservation law - sum of the K-columns=0 and
    %     (2) by the conditions of steady-state:  K*P=0
    unit_vector=ones(1,length(K(:,1)));
    off_equilibrium=sum(K*p');
    mass_conservation=sum(unit_vector*K);
    K_IMBALANCE_THRESHOLD=1e-6;
    if (off_equilibrium+mass_conservation)>K_IMBALANCE_THRESHOLD 
        fprintf(1,'\n\nERROR!\n K-matrix is inconsistent with equilibrium populations computation!');
        K
        p
        error('Check your K-matrix and equations for equilibrium concentrations');
    end
    
    
    %% Determine heat evolution for ITC curve
    if LR_counter>1
        dp=p-old_p;  % increment of populations of species
        dH=dp*Model.dH';  % add heat due to all increments
        if LR_counter==2
            ITC_curve(1)=dH; % assign 0-th point
        end
        ITC_curve(LR_counter)=dH;
    end
    
    % save populations for a next step
    old_p=p;
end
% show real points corresponding to line shapes
real_p=[];
for LR_counter=1:conditions.total_titration_points
    curr_LR=conditions.LRratio(LR_counter);
    curr_Rtotal= conditions.Rtotal(LR_counter);
    [p, K, R2_vector, v_vector, p_other, Model]=model_matrices(Model, curr_Rtotal, curr_LR);
    real_p= [ real_p;  p p_other];
end

% stop timer
fprintf(1,'Populations simulation done.\n');
toc;

figure_handle = PopulationsFigure(LRrange,smooth_p, real_p, Model, conditions, modes);
print(figure_handle,'-dpng','-r65',fullfile(conditions.figures_path,sprintf('%s_populations.png',Model.name)));
print(figure_handle,'-depsc2',fullfile(conditions.figures_path,sprintf('%s_populations.eps',Model.name)));
saveas(figure_handle, fullfile(conditions.figures_path,sprintf('%s_populations.fig',Model.name)));



%% Plot and save ITC curve
figure_handle = ITCFigure(LRrange, ITC_curve, Model, conditions);
print(figure_handle,'-dpng','-r65',fullfile(conditions.figures_path,sprintf('%s_ITC.png',Model.name)));
print(figure_handle,'-depsc2',fullfile(conditions.figures_path,sprintf('%s_ITC.eps',Model.name)));
saveas(figure_handle, fullfile(conditions.figures_path,sprintf('%s_ITC.fig',Model.name)));
ITCdata=[LRrange' ITC_curve];
save( fullfile(conditions.figures_path,sprintf('%s_ITC.txt',Model.name)), 'ITCdata','-ascii');





%% Compute spectra for every titration point 
tic;

w_range=conditions.w_min:  (conditions.w_max - conditions.w_min)/conditions.Spectral_Points : conditions.w_max;

spectral_series=cell(conditions.total_titration_points,1);

for point_counter=1:conditions.total_titration_points
    
    %-- Compute spectral intensity and save
    y=Spectra(Model, w_range, conditions.Rtotal(point_counter), conditions.LRratio(point_counter));
    spectral_series{point_counter}=y;

    %-- save data for import in LineShapeKin 3.1 as [ Hz, intensity]
    %   arrays
    full_spectral_slice=[(w_range/6.28)' y];
 
    
    % perform smart filtering to reduce CPU load in LineShapeKin:
    % make an array containing points in vertical intervals, dropping points if they are
    % not changing intensity - this way dropping points along the baseline
    intensity_intervals=100;
    curr_point=full_spectral_slice(1,:);
    spectral_slice=curr_point;
    interval_size=(max(y)-min(y))/intensity_intervals;
    current_level=curr_point(2);
    for raw_point_counter=2:length(y)
        curr_point=full_spectral_slice(raw_point_counter,:);
        if curr_point(2) > (current_level+interval_size) 
            % one new step up
            current_level=current_level+interval_size;
            spectral_slice=[spectral_slice; curr_point];
            continue
        end
        if curr_point(2) < (current_level-interval_size)
            % one new step down
            current_level=current_level-interval_size;
            spectral_slice=[spectral_slice; curr_point];
            continue
        end
        
    end
    
    % save for Lineshapekin 3.2
    file=fullfile(conditions.lineshapekin31_path,sprintf('%d.txt',point_counter));
    save(file,'spectral_slice','-ascii')
    % save for BiophysicsLab
    file=fullfile(conditions.BiophysicsLab_path,sprintf('%d.txt',point_counter));
    % add fake errors 0.1% of maximum intensity in the slice
    fake_errors=ones(size(spectral_slice,1),1)*max(spectral_slice(:,2))*0.001;
    spectral_slice_with_errors= [ spectral_slice fake_errors ];
    save(file,'spectral_slice_with_errors','-ascii')
     
    
end
fprintf(1,'Spectral simulation done.\n');
toc;

% Make  files needed by LineShapeKin 3
point_list=(1:1:conditions.total_titration_points)';
save(fullfile(conditions.lineshapekin31_path,'points_list.txt'), 'point_list','-ascii');
save(fullfile(conditions.BiophysicsLab_path,'points_list.txt'), 'point_list','-ascii');

a=1;
save(fullfile(conditions.lineshapekin31_path,'scale_factor.txt'), 'a','-ascii');

Rs=[];
save(fullfile(conditions.lineshapekin31_path,'manual_R2.txt'), 'Rs','-ascii');

ws=[];
save(fullfile(conditions.lineshapekin31_path,'manual_w0.txt'), 'ws','-ascii');

temp=conditions.Rtotal; % save cannot access fields in structures
save(fullfile(conditions.lineshapekin31_path,'Rtotal.txt'), 'temp','-ascii');
temp=conditions.LRratio;
save(fullfile(conditions.lineshapekin31_path,'LRratio.txt'), 'temp','-ascii');



%% Plot spectra and a chemical shift titration curve



[ spectral_figure, titration_curve_figure, titration_curve, line_width_curve_figure, line_width_curve]  = SpectralFigure(w_range, spectral_series, Model, conditions, modes);
% save spectra
print(spectral_figure,'-dpng','-r65',fullfile(conditions.figures_path,sprintf('%s_spectra.png',Model.name)));
print(spectral_figure,'-depsc2',fullfile(conditions.figures_path,sprintf('%s_spectra.eps',Model.name)));
saveas(spectral_figure, fullfile(conditions.figures_path,sprintf('%s_spectra.fig',Model.name)));
% save titration curve
print(titration_curve_figure,'-dpng','-r65',fullfile(conditions.figures_path,sprintf('%s_titration_curve.png',Model.name)));
print(titration_curve_figure,'-depsc2',fullfile(conditions.figures_path,sprintf('%s_titration_curve.eps',Model.name)));
saveas(titration_curve_figure, fullfile(conditions.figures_path,sprintf('%s_titration_curve.fig',Model.name)));
save( fullfile(conditions.figures_path,sprintf('%s_titration_curve.txt',Model.name)), 'titration_curve','-ascii');

% save line_width_curve_figure
print(line_width_curve_figure,'-dpng','-r65',fullfile(conditions.figures_path,sprintf('%s_line_width_curve.png',Model.name)));
print(line_width_curve_figure,'-depsc2',fullfile(conditions.figures_path,sprintf('%s_line_width_curve.eps',Model.name)));
saveas(line_width_curve_figure, fullfile(conditions.figures_path,sprintf('%s_line_width_curve.fig',Model.name)));
save( fullfile(conditions.figures_path,sprintf('%s_line_width_curve.txt',Model.name)), 'line_width_curve','-ascii');

 

%% Save model parameters

outstring=sprintf('%s\n', Model.name);
outstring=sprintf('%s%s\n', outstring,Model.description);
% output independent equilibrium constants
outstring=sprintf('%s\nEquilibrium constants:\n',outstring);
for counter=1:length(Model.Ka)
    outstring=sprintf('%sKa(%s)=%.2e\n', outstring,Model.Ka_names{counter},Model.Ka(counter));
end
% output dependent equilibrium constants
for counter=1:length(Model.Ka_dependent)
    outstring=sprintf('%sKa(%s)=%.2e (dependent)\n', outstring,Model.Ka_dependent_names{counter},Model.Ka_dependent(counter));
end

% output kinetic constants
outstring=sprintf('%s\nRate constants (1-forward, 2-reverse):\n',outstring);
for counter=1:length(Model.k2)
    outstring=sprintf('%sk1(%s)=%.2e \n', outstring,Model.k_names{counter},Model.k1(counter));
    outstring=sprintf('%sk2(%s)=%.2e \n', outstring,Model.k_names{counter},Model.k2(counter));
end

% output species chemical shifts
outstring=sprintf('%s\nChemical shifts:\n',outstring);
for counter=1:length(Model.w0)
    species=Model.species_names{counter};
    outstring=sprintf('%sw0(%s)=%.1f /s (%.1f Hz)\n', outstring,species,Model.w0(counter),Model.w0(counter)/6.28);
end

% output R2s
outstring=sprintf('%s\nBase relaxation rates:\n',outstring);
for counter=1:length(Model.R2)
    species=Model.species_names{counter};
    outstring=sprintf('%sR2(%s)=%.1f /s\n', outstring,species,Model.R2(counter));
end

% output dH
outstring=sprintf('%s\nEnthalpy difference from the base state:\n',outstring);
for counter=1:length(Model.dH)
    species=Model.species_names{counter};
    outstring=sprintf('%sdH(%s)=%.1f \n', outstring,species,Model.dH(counter));
end
% output Rtotal and LRratio
outstring=sprintf('%s\nTotal R concentration (*1000):\n', outstring);
for counter=1:length(conditions.Rtotal)
    outstring=sprintf('%s%.2f  ', outstring,1000*conditions.Rtotal(counter));
end
outstring=sprintf('%s\n\nRatio of total L to total R:\n', outstring);
for counter=1:length(conditions.LRratio)
    outstring=sprintf('%s%.2f  ', outstring,conditions.LRratio(counter));
end

%-- save in a text file
fid=fopen(fullfile(conditions.figures_path,sprintf('%s.txt',Model.name)),'w');
fprintf(fid,'%s',outstring);
fclose(fid);

%-- save as image for HTML display
%figure_handle= TextFigure(outstring);
%print(figure_handle,'-dpng','-r75',fullfile(conditions.figures_path,sprintf('%s_parameters.png',Model.name)));

%-- generate HTML report
system(sprintf('lineshapekin_simulation_report.py %s', conditions.results_path));
close all;

