start_simulation.m 2.12 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;
%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]=start_simulation(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);

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

for k=1:Nr
    
  waitbar(k/Nr,h,['Realization number ',num2str(k)]);
  
  hypercube=LatinHypercube(LBpi,UBpi,Ns);
  perturbation=hypercube.generate_sample(num_parameters);
  perturbed_p=repmat(nominal_parameters,length(perturbation),1).*perturbation;
  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);  
  parfor i=1:Ns
     pi=perturbed_p(i,:);
     modelsim=@(t,x)feval(model.name,t,x,pi,model.input,model.total_proteins);
     [tsimi,psimi]=feval(model.odesolver,modelsim,model.time,model.initial_conditions);
    % for t=1:model.num_observables
     %  psim_all{i,t}=psimi(:,t);
     %end
     msgstr=lastwarn;
     if isequal(msgstr,'')        
        psim_all{i,1}=psimi;
     else 
        psim_all{i,1}=NaN(length(model.time),model.num_observables);
     end
      lastwarn('')
   end 
  
 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