start_simulation_v2_apr2020.m 2.61 KB
%Function for simulating the ODE model using the perturbed parameter
%vector. It takes in input:
%1. model which is the struct containing all the characheristics of the
%model
%2. Nr is the number of independent realizations to perform;
%3. LBpi and UBpi are the lower and upper boundaries of the Latin
%Hypercube. They are arrays since each parameter may have its own range of
%variation;
%4. Ns is the number of samples to generate for creating the Latin
%Hypercube;
%It returns in output the array AllResults containing all the output
%variables simulated using the parameter vectors in AllPerturbations

function [AllResults, AllPerturbations,AllDerivedParam]=start_simulation_v2_apr2020(model,Nr,LBpi,UBpi,Ns)

nominal_parameters=model.nominal_parameters;
num_parameters=length(nominal_parameters);

AllResults=cell(Nr,model.num_observables);

AllPerturbations=cell(Nr,1);

AllDerivedParam=cell(Nr,1);

h = waitbar(0,'Please wait...');

for k=1:Nr
    
  waitbar(k/Nr,h,['Realization number ',num2str(k)]);
  
  perturbation=zeros(Ns,num_parameters);
  for i=1:length(LBpi)
     hypercube=LatinHypercube(LBpi(i),UBpi(i),Ns);
     perturbation_i=hypercube.generate_sample(1);
     perturbation(:,i)=perturbation_i;
  end  
  perturbed_p=repmat(nominal_parameters,length(perturbation),1).*perturbation;
  size(perturbed_p)
  AllPerturbations{k,1}=perturbation;
  disp(strcat('Generated a Latin Hypercube of size',{' '},sprintf('%.0f',size(perturbation,1)),'x',sprintf('%.0f',size(perturbation,2))));
    %psim_all=cell(Ns,model.num_observables);
  psim_all=cell(Ns,1);  
  R0_all=zeros(Ns,1);
  parfor i=1:Ns
     pi=perturbed_p(i,:);
     %modelsim=@(t,x)feval(model.name,t,x,pi,model.input,model.total_proteins);
     modelsim=@(t,x)feval(model.name,t,x,pi,model.Tlock1,model.s01,model.s11);
     [tsimi,psimi]=feval(model.odesolver,modelsim,model.time,model.initial_conditions);
     R0i=compute_R0(pi,model.Tlock1,model.s01,model.s11);
    % for t=1:model.num_observables
     %  psim_all{i,t}=psimi(:,t);
     %end
     msgstr=lastwarn;
     if isequal(msgstr,'')        
        psim_all{i,1}=psimi;
        R0_all(i,1)=R0i;
     else 
        psim_all{i,1}=NaN(length(model.time),model.num_observables);
        R0_all{i,1}=NaN;
     end
      lastwarn('')
  end 
   
 AllDerivedParam{k,1}=R0_all;
  
 for j=1:model.num_observables
    for w=1:Ns
     AllResults{k,j}{w,1}=psim_all{w,1}(:,j);
    end
  end
  %AllResults{k,1:end}=psim_all;
  disp(strcat('Completed one realization. The array of results has size',{' '},sprintf('%.0f',size(psim_all,1)),'x',sprintf('%.0f',size(psim_all,2))));
 end

delete(gcp) %close parallel pool
close (h)  %close the waitbar
end