% data=dsSimulate(model,'option',value,...)
Purpose: manage simulation of a DynaSim model. This high-level function
offers many options to control the number of simulations, how the model is
optionally varied across simulations, and how the numerical integration
is performed. It can optionally save the simulated data and create/submit
simulation jobs to a compute cluster.
Inputs:
model: DynaSim model structure or equations (see dsGenerateModel and
dsCheckModel for more details)
solver options (provided as key/value pairs: 'option1',value1,'option2',value2,...):
'solver' : solver for numerical integration (see dsGetSolveFile)
{'euler','rk2','rk4', or any built-in matlab solver} (default: 'rk4')
'tspan' : time limits of simulation [begin,end] (default: [0 100]) [ms]
note: units must be consistent with dt and model equations
'dt' : time step used for DynaSim solvers (default: .01) [ms]
'downsample_factor': downsampling applied during simulation (default: 1, no downsampling)
(only every downsample_factor-time point is stored in memory and/or written to disk)
'ic' : numeric array of initial conditions, one value per state
variable (default: all zeros). overrides definition in model structure
'random_seed' : seed for random number generator (default: 'shuffle', set randomly) (usage: rng(options.random_seed))
'compile_flag': whether to compile simulation using coder instead of
interpreting Matlab {0 or 1} (default: 0)
'sparse_flag' : whether to convert numeric fixed variables to sparse matrices {0 or 1} (default: 0)
options for running sets of simulations:
'vary' : (default: [], vary nothing): cell matrix specifying model
components to vary across simulations (see NOTE 1 and dsVary2Modifications)
options to control saved data:
'matCompatibility_flag': whether to save mat files in compatible mode, vs to prioritize > 2GB VARs {0 or 1} (default: 1)
'save_results_flag': whether to save results of analysis and plotting
'save_data_flag': whether to save simulated data to disk after completion {0 or 1} (default: 0)
'overwrite_flag': whether to overwrite existing data files {0 or 1} (default: 0)
'study_dir' : relative or absolute path to output directory (default: current directory)
'prefix' : string to prepend to all output file names (default: 'study')
'disk_flag' : whether to write to disk during simulation instead of storing in memory {0 or 1} (default: 0)
'precision' : {'single','double'} precision of simulated data saved to disk (default: 'single')
options for cluster computing:
'cluster_flag' : whether to run simulations on a cluster submitted
using qsub (see dsCreateBatch) {0 or 1} (default: 0)
'sims_per_job' : number of simulations to run per batch job (default: 1)
'memory_limit' : memory to allocate per batch job (default: '8G')
'qsub_mode' : whether to use SGE -t array for 1 qsub, mode: 'array'; or
qsub in csh for loop, mode: 'loop'. (default: 'loop').
'one_solve_file_flag': only use 1 file of each time when solving (default: 0)
'optimize_big_vary': Select best options for doing many sims {0 or 1} (default: 0)
options for parallel computing: (requires Parallel Computing Toolbox)
'parallel_flag' : whether to use parfor to run simulations {0 or 1} (default: 0)
'num_cores' : number of cores to specify in the parallel pool
*note: parallel computing has been disabled for debugging...
options for post-processing:
'analysis_functions': cell array of analysis function handles
'analysis_options' : cell array of option cell arrays {'option1',value1,...}
'plot_functions' : cell array of plot function handles
'plot_options' : cell array of option cell arrays {'option1',value1,...}
other options:
'verbose_flag' : whether to display informative messages/logs (default: 0)
'modifications' : how to modify DynaSim specification structure component before simulation (see dsApplyModifications)
'experiment' : function handle of experiment function (see NOTE 2)
'experiment_options' : single cell array of key/value options for experiment function
'optimization' : function handle of optimization function (see NOTE 2)
'debug_flag' : set to debug mode
'benchmark_flag': set to benchmark mode. will add tic/toc to sims.
Outputs:
DynaSim data structure:
data.labels : list of state variables and monitors recorded
data.(state_variables): state variable data matrix [time x cells]
data.(monitors) : monitor data matrix [time x cells]
data.time : time vector [time x 1]
data.simulator_options: simulator options used to generate simulated data
data.model : model used to generate simulated data
[data.varied] : list of varied model components (present only if anything was varied)
DynaSim studyinfo structure (only showing select fields, see dsCheckStudyinfo for more details)
studyinfo.study_dir
studyinfo.base_model (=[]): original model from which a set of simulations was derived
studyinfo.base_simulator_options (=[])
studyinfo.base_solve_file (='')
studyinfo.simulations(k): metadata for each simulation in a set of simulations
.sim_id : unique identifier in study
.modifications : modifications made to the base model during this simulation
.data_file : full filename of eventual output file
.batch_dir (=[]): directory where batch jobs were saved (if cluster_flag=1)
.job_file (=[]) : m-file batch job that runs this simulation (if cluster_flag=1)
.simulator_options: simulator options for this simulation
.solve_file : full filename of m- or mex-file that numerically integrated the model
NOTE 1: 'vary' indicates the variable to vary, the values
it should take, and the object whose variable should be varied.
Syntax: vary={object, variable, values; ...}. For instance, to vary
parameter 'gNa', taking on values 100 and 120, in population 'E', set
vary={'E','gNa',[100 120]}. To additionally vary 'gSYN' in the connection
mechanism from 'E' to 'I', set vary={'E','gNa',[100 120];'E->I','gSYN',[0 1]}.
Mechanism lists and equations can also be varied. (see dsVary2Modifications
for more details and examples).
EXAMPLES:
Example 1: Lorenz equations with phase plot
eqns={
's=10; r=27; b=2.666';
'dx/dt=s*(y-x)';
'dy/dt=r*x-y-x*z';
'dz/dt=-b*z+x*y';
};
data=dsSimulate(eqns,'tspan',[0 100],'ic',[1 2 .5]);
plot(data.pop1_x,data.pop1_z); title('Lorenz equations'); xlabel('x'); ylabel('z')
Example 2: Leaky integrate-and-fire with spike monitor
eqns={
'tau=10; R=10; E=-70; I=1.55; thresh=-55; reset=-75';
'dV/dt=(E-V+R*I)/tau; if(V>thresh)(V=reset)';
'monitor V.spikes(thresh)';
};
data=dsSimulate(eqns,'tspan',[0 200],'ic',-75);
data.pop1_V(data.pop1_V_spikes==1)=20; % insert spike
plot(data.time,data.pop1_V); xlabel('time (ms)'); ylabel('V'); title('LIF with spikes')
Example 3: Hodgkin-Huxley-type Intrinsically Bursting neuron
eqns='dv/dt=5+@current; {iNaF,iKDR,iM}; gNaF=100; gKDR=5; gM=1.5; v(0)=-70';
data=dsSimulate(eqns,'tspan',[0 200]);
figure; plot(data.time,data.(data.labels{1}))
xlabel('time (ms)'); ylabel('membrane potential (mV)'); title('Intrinsically Bursting neuron')
Example 4: varying max Na+ conductance in Hodgkin-Huxley neuron
eqns='dv/dt=@current+10; {iNa,iK}; v(0)=-60';
data=dsSimulate(eqns,'vary',{'','gNa',[50 100 200]});
% plot how mean firing rate varies with parameter
dsPlotFR(data,'bin_size',30,'bin_shift',10); % bin_size and bin_shift in [ms]
Example 5: simulate predefined Hodgkin-Huxley neuron with Iapp=10
data = dsSimulate('HH','vary',{'HH','Iapp',10});
figure; plot(data.time,data.HH_V);
Example 6: Sparse Pyramidal-Interneuron-Network-Gamma rhythm with rastergram
% define equations of cell model (same for E and I populations)
eqns={
'dv/dt=Iapp+@current/Cm+noise*randn(1,N_pop)*sqrt(dt)/dt';
'monitor v.spikes, iGABAa.functions, iAMPA.functions'
};
% define specification for two-population network model
s=[];
s.populations(1).name='E';
s.populations(1).size=80;
s.populations(1).equations=eqns;
s.populations(1).mechanism_list={'iNa','iK'};
s.populations(1).parameters={'Iapp',5,'gNa',120,'gK',36,'Cm',1,'noise',4};
s.populations(2).name='I';
s.populations(2).size=20;
s.populations(2).equations=eqns;
s.populations(2).mechanism_list={'iNa','iK'};
s.populations(2).parameters={'Iapp',0,'gNa',120,'gK',36,'Cm',1,'noise',4};
s.connections(1).source='I';
s.connections(1).target='E';
s.connections(1).mechanism_list={'iGABAa'};
s.connections(1).parameters={'tauD',10,'gSYN',.1,'netcon','ones(N_pre,N_post)'};
s.connections(2).source='E';
s.connections(2).target='I';
s.connections(2).mechanism_list={'iAMPA'};
s.connections(2).parameters={'tauD',2,'gSYN',.1,'netcon',ones(80,20)};
% simulate model
data=dsSimulate(s);
% plot voltages and rastergram
figure;
subplot(2,1,1); % voltage traces
plot(data.time,data.E_v,'b-',data.time,data.I_v,'r-')
title('Sparse Pyramidal-Interneuron-Network-Gamma (sPING)'); ylabel('membrane potential (mV)');
subplot(2,1,2); % rastergram
E_spikes=nan(size(data.E_v_spikes)); E_spikes(data.E_v_spikes==1)=1;
I_spikes=nan(size(data.I_v_spikes)); I_spikes(data.I_v_spikes==1)=1;
plot(data.time,E_spikes+repmat(1:80,[length(data.time) 1]),'bo'); hold on
plot(data.time,I_spikes+repmat(80+(1:20),[length(data.time) 1]),'ro'); axis([0 100 0 100]);
title('rastergram'); xlabel('time (ms)'); ylabel('cell index');
% simulate model varying two parameters (Iapp and tauD in sPING)
% warning: this may take up to a minute to complete:
vary={
'E' ,'Iapp',[0 10 20]; % amplitude of tonic input to E-cells
'I->E','tauD',[5 10 15] % inhibition decay time constant from I to E
};
data=dsSimulate(s,'vary',vary);
% plot firing rates calculated from spike monitor in both populations
dsPlotFR(data,'variable','*_spikes','bin_size',30,'bin_shift',10);
See also: dsGenerateModel, dsCheckModel, dsGetSolveFile, dsCheckData,
dsVary2Modifications, dsCheckStudyinfo, dsCreateBatch
Author: Jason Sherfey, PhD <jssherfey@gmail.com>
Copyright (C) 2016 Jason Sherfey, Boston University, USA0001 function [data,studyinfo,result] = dsSimulate(model,varargin) 0002 %% data=dsSimulate(model,'option',value,...) 0003 % Purpose: manage simulation of a DynaSim model. This high-level function 0004 % offers many options to control the number of simulations, how the model is 0005 % optionally varied across simulations, and how the numerical integration 0006 % is performed. It can optionally save the simulated data and create/submit 0007 % simulation jobs to a compute cluster. 0008 % Inputs: 0009 % model: DynaSim model structure or equations (see dsGenerateModel and 0010 % dsCheckModel for more details) 0011 % 0012 % solver options (provided as key/value pairs: 'option1',value1,'option2',value2,...): 0013 % 'solver' : solver for numerical integration (see dsGetSolveFile) 0014 % {'euler','rk2','rk4', or any built-in matlab solver} (default: 'rk4') 0015 % 'tspan' : time limits of simulation [begin,end] (default: [0 100]) [ms] 0016 % note: units must be consistent with dt and model equations 0017 % 'dt' : time step used for DynaSim solvers (default: .01) [ms] 0018 % 'downsample_factor': downsampling applied during simulation (default: 1, no downsampling) 0019 % (only every downsample_factor-time point is stored in memory and/or written to disk) 0020 % 'ic' : numeric array of initial conditions, one value per state 0021 % variable (default: all zeros). overrides definition in model structure 0022 % 'random_seed' : seed for random number generator (default: 'shuffle', set randomly) (usage: rng(options.random_seed)) 0023 % 'compile_flag': whether to compile simulation using coder instead of 0024 % interpreting Matlab {0 or 1} (default: 0) 0025 % 'sparse_flag' : whether to convert numeric fixed variables to sparse matrices {0 or 1} (default: 0) 0026 % 0027 % options for running sets of simulations: 0028 % 'vary' : (default: [], vary nothing): cell matrix specifying model 0029 % components to vary across simulations (see NOTE 1 and dsVary2Modifications) 0030 % 0031 % options to control saved data: 0032 % 'matCompatibility_flag': whether to save mat files in compatible mode, vs to prioritize > 2GB VARs {0 or 1} (default: 1) 0033 % 'save_results_flag': whether to save results of analysis and plotting 0034 % 'save_data_flag': whether to save simulated data to disk after completion {0 or 1} (default: 0) 0035 % 'overwrite_flag': whether to overwrite existing data files {0 or 1} (default: 0) 0036 % 'study_dir' : relative or absolute path to output directory (default: current directory) 0037 % 'prefix' : string to prepend to all output file names (default: 'study') 0038 % 'disk_flag' : whether to write to disk during simulation instead of storing in memory {0 or 1} (default: 0) 0039 % 'precision' : {'single','double'} precision of simulated data saved to disk (default: 'single') 0040 % 0041 % options for cluster computing: 0042 % 'cluster_flag' : whether to run simulations on a cluster submitted 0043 % using qsub (see dsCreateBatch) {0 or 1} (default: 0) 0044 % 'sims_per_job' : number of simulations to run per batch job (default: 1) 0045 % 'memory_limit' : memory to allocate per batch job (default: '8G') 0046 % 'qsub_mode' : whether to use SGE -t array for 1 qsub, mode: 'array'; or 0047 % qsub in csh for loop, mode: 'loop'. (default: 'loop'). 0048 % 'one_solve_file_flag': only use 1 file of each time when solving (default: 0) 0049 % 'optimize_big_vary': Select best options for doing many sims {0 or 1} (default: 0) 0050 % 0051 % options for parallel computing: (requires Parallel Computing Toolbox) 0052 % 'parallel_flag' : whether to use parfor to run simulations {0 or 1} (default: 0) 0053 % 'num_cores' : number of cores to specify in the parallel pool 0054 % *note: parallel computing has been disabled for debugging... 0055 % 0056 % options for post-processing: 0057 % 'analysis_functions': cell array of analysis function handles 0058 % 'analysis_options' : cell array of option cell arrays {'option1',value1,...} 0059 % 'plot_functions' : cell array of plot function handles 0060 % 'plot_options' : cell array of option cell arrays {'option1',value1,...} 0061 % 0062 % other options: 0063 % 'verbose_flag' : whether to display informative messages/logs (default: 0) 0064 % 'modifications' : how to modify DynaSim specification structure component before simulation (see dsApplyModifications) 0065 % 'experiment' : function handle of experiment function (see NOTE 2) 0066 % 'experiment_options' : single cell array of key/value options for experiment function 0067 % 'optimization' : function handle of optimization function (see NOTE 2) 0068 % 'debug_flag' : set to debug mode 0069 % 'benchmark_flag': set to benchmark mode. will add tic/toc to sims. 0070 % 0071 % Outputs: 0072 % DynaSim data structure: 0073 % data.labels : list of state variables and monitors recorded 0074 % data.(state_variables): state variable data matrix [time x cells] 0075 % data.(monitors) : monitor data matrix [time x cells] 0076 % data.time : time vector [time x 1] 0077 % data.simulator_options: simulator options used to generate simulated data 0078 % data.model : model used to generate simulated data 0079 % [data.varied] : list of varied model components (present only if anything was varied) 0080 % 0081 % DynaSim studyinfo structure (only showing select fields, see dsCheckStudyinfo for more details) 0082 % studyinfo.study_dir 0083 % studyinfo.base_model (=[]): original model from which a set of simulations was derived 0084 % studyinfo.base_simulator_options (=[]) 0085 % studyinfo.base_solve_file (='') 0086 % studyinfo.simulations(k): metadata for each simulation in a set of simulations 0087 % .sim_id : unique identifier in study 0088 % .modifications : modifications made to the base model during this simulation 0089 % .data_file : full filename of eventual output file 0090 % .batch_dir (=[]): directory where batch jobs were saved (if cluster_flag=1) 0091 % .job_file (=[]) : m-file batch job that runs this simulation (if cluster_flag=1) 0092 % .simulator_options: simulator options for this simulation 0093 % .solve_file : full filename of m- or mex-file that numerically integrated the model 0094 % 0095 % NOTE 1: 'vary' indicates the variable to vary, the values 0096 % it should take, and the object whose variable should be varied. 0097 % Syntax: vary={object, variable, values; ...}. For instance, to vary 0098 % parameter 'gNa', taking on values 100 and 120, in population 'E', set 0099 % vary={'E','gNa',[100 120]}. To additionally vary 'gSYN' in the connection 0100 % mechanism from 'E' to 'I', set vary={'E','gNa',[100 120];'E->I','gSYN',[0 1]}. 0101 % Mechanism lists and equations can also be varied. (see dsVary2Modifications 0102 % for more details and examples). 0103 % 0104 % EXAMPLES: 0105 % Example 1: Lorenz equations with phase plot 0106 % eqns={ 0107 % 's=10; r=27; b=2.666'; 0108 % 'dx/dt=s*(y-x)'; 0109 % 'dy/dt=r*x-y-x*z'; 0110 % 'dz/dt=-b*z+x*y'; 0111 % }; 0112 % data=dsSimulate(eqns,'tspan',[0 100],'ic',[1 2 .5]); 0113 % plot(data.pop1_x,data.pop1_z); title('Lorenz equations'); xlabel('x'); ylabel('z') 0114 % 0115 % Example 2: Leaky integrate-and-fire with spike monitor 0116 % eqns={ 0117 % 'tau=10; R=10; E=-70; I=1.55; thresh=-55; reset=-75'; 0118 % 'dV/dt=(E-V+R*I)/tau; if(V>thresh)(V=reset)'; 0119 % 'monitor V.spikes(thresh)'; 0120 % }; 0121 % data=dsSimulate(eqns,'tspan',[0 200],'ic',-75); 0122 % data.pop1_V(data.pop1_V_spikes==1)=20; % insert spike 0123 % plot(data.time,data.pop1_V); xlabel('time (ms)'); ylabel('V'); title('LIF with spikes') 0124 % 0125 % Example 3: Hodgkin-Huxley-type Intrinsically Bursting neuron 0126 % eqns='dv/dt=5+@current; {iNaF,iKDR,iM}; gNaF=100; gKDR=5; gM=1.5; v(0)=-70'; 0127 % data=dsSimulate(eqns,'tspan',[0 200]); 0128 % figure; plot(data.time,data.(data.labels{1})) 0129 % xlabel('time (ms)'); ylabel('membrane potential (mV)'); title('Intrinsically Bursting neuron') 0130 % 0131 % Example 4: varying max Na+ conductance in Hodgkin-Huxley neuron 0132 % eqns='dv/dt=@current+10; {iNa,iK}; v(0)=-60'; 0133 % data=dsSimulate(eqns,'vary',{'','gNa',[50 100 200]}); 0134 % % plot how mean firing rate varies with parameter 0135 % dsPlotFR(data,'bin_size',30,'bin_shift',10); % bin_size and bin_shift in [ms] 0136 % 0137 % Example 5: simulate predefined Hodgkin-Huxley neuron with Iapp=10 0138 % data = dsSimulate('HH','vary',{'HH','Iapp',10}); 0139 % figure; plot(data.time,data.HH_V); 0140 % 0141 % Example 6: Sparse Pyramidal-Interneuron-Network-Gamma rhythm with rastergram 0142 % % define equations of cell model (same for E and I populations) 0143 % eqns={ 0144 % 'dv/dt=Iapp+@current/Cm+noise*randn(1,N_pop)*sqrt(dt)/dt'; 0145 % 'monitor v.spikes, iGABAa.functions, iAMPA.functions' 0146 % }; 0147 % % define specification for two-population network model 0148 % s=[]; 0149 % s.populations(1).name='E'; 0150 % s.populations(1).size=80; 0151 % s.populations(1).equations=eqns; 0152 % s.populations(1).mechanism_list={'iNa','iK'}; 0153 % s.populations(1).parameters={'Iapp',5,'gNa',120,'gK',36,'Cm',1,'noise',4}; 0154 % s.populations(2).name='I'; 0155 % s.populations(2).size=20; 0156 % s.populations(2).equations=eqns; 0157 % s.populations(2).mechanism_list={'iNa','iK'}; 0158 % s.populations(2).parameters={'Iapp',0,'gNa',120,'gK',36,'Cm',1,'noise',4}; 0159 % s.connections(1).source='I'; 0160 % s.connections(1).target='E'; 0161 % s.connections(1).mechanism_list={'iGABAa'}; 0162 % s.connections(1).parameters={'tauD',10,'gSYN',.1,'netcon','ones(N_pre,N_post)'}; 0163 % s.connections(2).source='E'; 0164 % s.connections(2).target='I'; 0165 % s.connections(2).mechanism_list={'iAMPA'}; 0166 % s.connections(2).parameters={'tauD',2,'gSYN',.1,'netcon',ones(80,20)}; 0167 % % simulate model 0168 % data=dsSimulate(s); 0169 % % plot voltages and rastergram 0170 % figure; 0171 % subplot(2,1,1); % voltage traces 0172 % plot(data.time,data.E_v,'b-',data.time,data.I_v,'r-') 0173 % title('Sparse Pyramidal-Interneuron-Network-Gamma (sPING)'); ylabel('membrane potential (mV)'); 0174 % subplot(2,1,2); % rastergram 0175 % E_spikes=nan(size(data.E_v_spikes)); E_spikes(data.E_v_spikes==1)=1; 0176 % I_spikes=nan(size(data.I_v_spikes)); I_spikes(data.I_v_spikes==1)=1; 0177 % plot(data.time,E_spikes+repmat(1:80,[length(data.time) 1]),'bo'); hold on 0178 % plot(data.time,I_spikes+repmat(80+(1:20),[length(data.time) 1]),'ro'); axis([0 100 0 100]); 0179 % title('rastergram'); xlabel('time (ms)'); ylabel('cell index'); 0180 % % simulate model varying two parameters (Iapp and tauD in sPING) 0181 % % warning: this may take up to a minute to complete: 0182 % vary={ 0183 % 'E' ,'Iapp',[0 10 20]; % amplitude of tonic input to E-cells 0184 % 'I->E','tauD',[5 10 15] % inhibition decay time constant from I to E 0185 % }; 0186 % data=dsSimulate(s,'vary',vary); 0187 % % plot firing rates calculated from spike monitor in both populations 0188 % dsPlotFR(data,'variable','*_spikes','bin_size',30,'bin_shift',10); 0189 % 0190 % 0191 % See also: dsGenerateModel, dsCheckModel, dsGetSolveFile, dsCheckData, 0192 % dsVary2Modifications, dsCheckStudyinfo, dsCreateBatch 0193 % 0194 % Author: Jason Sherfey, PhD <jssherfey@gmail.com> 0195 % Copyright (C) 2016 Jason Sherfey, Boston University, USA 0196 0197 % TODO: rename 'disk_flag' to something more descriptive 0198 0199 % dependencies: dsWriteDynaSimSolver, dsWriteMatlabSolver, dsPropagateFunctions, dsCheckModel, 0200 % dsCheckOptions, dsOptions2Keyval, displayError, DynaSim2Odefun 0201 0202 % <-- temporarily removed from help section --> 0203 % NOTE 2: special functions that recursively call dsSimulate: 0204 % - "Experiments" are ways of hacking the ODE system to incorporate additional 0205 % models (e.g., controlled inputs) and use them to simulate experimental 0206 % protocols by systematically varying Model Components across simulations in 0207 % prescribed ways. Technically, an Experiment could be any function that takes 0208 % a DynaSim model structure as its first input followed by key/value options. 0209 % Ideally, Experiments represent standardized procedural methods (experimental 0210 % protocols) for studying the modeled system. Experiment functions typically 0211 % involve applying a set of Modifications to a Base Model and varying the 0212 % modified model in prescribed ways. 0213 % - during "Optimization", each iteration involves a single Study producing 0214 % modified models, their simulated data sets and analysis results (e.g., 0215 % cost functions) that shape the Base Model for a subsequent iteration and 0216 % its Study. Hence, a closed-loop optimization protocol produces a set of 0217 % evolving Studies. Technically, an Optimization could be any function 0218 % that takes a DynaSim model structure as its first input followed by 0219 % key/value options. Optimization functions will typically involve 0220 % while-looping through multiple Studies and analyzing data sets to update 0221 % Model Components on each iteration until some stop condition is reached. 0222 0223 %% 0.0 Outputs & inputs. 0224 % Initialize outputs 0225 data=[]; 0226 studyinfo=[]; 0227 0228 % check path 0229 dynasim_path=fileparts(which(mfilename)); 0230 onPath=~isempty(strfind(path,[dynasim_path, pathsep])); 0231 if ~onPath 0232 if 1 0233 fprintf('adding dynasim and sub-directory to Matlab path: %s\n',dynasim_path); 0234 end 0235 addpath(genpath(dynasim_path)); % necessary b/c of changing directory for simulation 0236 end 0237 0238 % Check inputs 0239 varargin = backward_compatibility(varargin); 0240 options=dsCheckOptions(varargin,{... 0241 'tspan',[0 100],[],... % [beg,end] (units must be consistent with dt and equations) 0242 'ic',[],[],... % initial conditions (overrides definition in model structure; can input as IC structure or numeric array) 0243 'solver','rk4',{'euler','rk1','rk2','rk4','modified_euler','rungekutta','rk',... 0244 'ode23','ode45','ode113','ode15s','ode23s','ode23t','ode23tb'},... % DynaSim and built-in Matlab solvers 0245 'matlab_solver_options',[],[],... % options from odeset for use with built-in Matlab solvers 0246 'dt',.01,[],... % time step used for fixed step DynaSim solvers 0247 'downsample_factor',1,[],... % downsampling applied during simulation (only every downsample_factor-time point is stored in memory or written to disk) 0248 'reduce_function_calls_flag',1,{0,1},... % whether to eliminate internal (anonymous) function calls 0249 'save_parameters_flag',1,{0,1},... 0250 'random_seed','shuffle',[],... % seed for random number generator (usage: rng(random_seed)) 0251 'data_file','data.csv',[],... % name of data file if disk_flag=1 0252 'precision','single',{'single','double'},... 0253 'logfid',1,[],... 0254 'store_model_flag',1,{0,1},... % whether to store model structure with data 0255 'verbose_flag',0,{0,1},... 0256 'modifications',[],[],... % *DynaSim modifications structure 0257 'vary',[],[],... % specification of things to vary or custom modifications_set 0258 'experiment',[],[],... % experiment function. func(model,args) 0259 'experiment_options',[],[],... 0260 'optimization',[],[],... 0261 'cluster_flag',0,{0,1},... % whether to run simulations on a cluster 0262 'sims_per_job',1,[],... % how many sims to run per batch job 0263 'memory_limit','8G',[],... % how much memory to allocate per batch job 0264 'qsub_mode','loop',{'loop','array'},... % whether to submit jobs as an array using qsub -t or in a for loop 0265 'one_solve_file_flag',0,{0,1},... % use only 1 solve file of each type, but can't vary mechs yet 0266 'parallel_flag',0,{0,1},... % whether to run simulations in parallel (using parfor) 0267 'num_cores',4,[],... % # cores for parallel processing (SCC supports 1-12) 0268 'compile_flag',0,{0,1},... % exist('codegen')==6, whether to compile using coder instead of interpreting Matlab 0269 'sparse_flag',0,{0,1},... % whether to sparsify fixed variables before simulation 0270 'disk_flag',0,{0,1},... % whether to write to disk during simulation instead of storing in memory 0271 'matCompatibility_flag',1,{0,1},... % whether to save mat files in compatible mode, vs to prioritize > 2GB VARs 0272 'save_data_flag',0,{0,1},... % whether to save simulated data 0273 'save_results_flag',0,{0,1},... % whether to save results from simulated data 0274 'project_dir',pwd,[],... 0275 'study_dir',[],[],... % study directory 0276 'prefix','study',[],... % prefix prepended to all output files 0277 'overwrite_flag',0,{0,1},... % whether to overwrite existing data 0278 'solve_file',[],[],... % m- or mex-file solving the system 0279 'sim_id',[],[],... % sim id in an existing study 0280 'studyinfo',[],[],... 0281 'email',[],[],... % email to send notification upon study completion 0282 'analysis_functions',[],[],... 0283 'analysis_options',[],[],... 0284 'plot_functions',[],[],... 0285 'plot_options',[],[],... 0286 'optimize_big_vary',0,{0,1},... 0287 'mex_dir_flag',1,{0,1},... % Flag to tell whether or not to search in mex_dir for pre-compiled solve files (solve*_mex*). 0288 'mex_dir',[],[],... % Directory to search for pre-compiled mex files. Can be relative to 'study_dir' or absolute path. 0289 'in_parfor_loop_flag',0,{0,1},... % if inside parfor loop 0290 'debug_flag',0,{0,1},... 0291 'auto_gen_test_data_flag',0,{0,1},... 0292 'unit_test_flag',0,{0,1},... 0293 'benchmark_flag',0,{0,1},... 0294 },false); 0295 % more options: remove_solve_dir, remove_batch_dir, post_downsample_factor 0296 0297 %% 0.1 Auto_gen_test_data_flag save argin (arguments passed to dsSimulate). 0298 if options.auto_gen_test_data_flag 0299 varargs = varargin; 0300 varargs{find(strcmp(varargs, 'auto_gen_test_data_flag'))+1} = 0; 0301 varargs(end+1:end+2) = {'unit_test_flag',1}; 0302 argin = [{model}, varargs]; % specific to this function 0303 end 0304 0305 %% 0.2 Unit testing. 0306 if options.unit_test_flag 0307 options.study_dir = getAbsolutePath(options.study_dir); 0308 end 0309 0310 %% 0.3 Prepare solve options. 0311 0312 if options.compile_flag && ~strcmp(reportUI,'matlab') 0313 fprintf('Setting ''compile_flag'' to 0 in Octave.\n') 0314 options.compile_flag = 0; 0315 end 0316 0317 if isempty(options.mex_dir) 0318 options.mex_dir = dsGetConfig('mex_path'); 0319 end 0320 0321 if options.compile_flag && options.sparse_flag 0322 error('The Matlab Coder toolbox does not support sparse matrices. Choose either ''compile_flag'' or ''sparse_flag''.'); 0323 end 0324 0325 if options.parallel_flag && (strcmp(reportUI,'matlab') && feature('numCores') == 1 || ~strcmp(reportUI,'matlab') && nproc() == 1) 0326 fprintf('Setting ''parallel_flag'' to 0 since only 1 core detected on this machine.\n') 0327 options.parallel_flag = 0; 0328 end 0329 0330 if options.compile_flag && ~options.reduce_function_calls_flag 0331 fprintf('Setting ''reduce_function_calls_flag'' to 1 for compatibility with ''compile_flag=1'' (coder does not support anonymous functions).\n'); 0332 options.reduce_function_calls_flag=1; 0333 end 0334 0335 % Make sure that data is either saved to disk, saved to variable, or plotted 0336 if (nargout==0 || options.cluster_flag) && ~options.save_data_flag && ~options.save_results_flag && isempty(options.plot_functions) 0337 if ~isempty(options.analysis_functions) 0338 fprintf('Setting ''save_results_flag'' to 1 since output from dsSimulate is not stored in a variable and analysis functions specified.\n') 0339 options.save_results_flag = 1; 0340 else 0341 fprintf('Setting ''save_data_flag'' to 1 since output from dsSimulate is not stored in a variable and no plot or analysis functions specified.\n') 0342 options.save_data_flag = 1; 0343 end 0344 end 0345 0346 if options.cluster_flag && ~options.save_data_flag 0347 options.save_results_flag=1; 0348 if options.verbose_flag 0349 fprintf('Setting ''save_results_flag'' to 1 for storing results of batch jobs for later access.\n'); 0350 end 0351 end 0352 0353 if any(strcmp(options.solver, {'ode23','ode45','ode113','ode15s','ode23s','ode23t','ode23tb'})) 0354 matlabSolverBool = 1; 0355 else 0356 matlabSolverBool = 0; 0357 end 0358 0359 if options.disk_flag && matlabSolverBool 0360 fprintf('Since using built-in solver, setting options.disk_flag=1.\n'); 0361 options.disk_flag = 1; 0362 end 0363 0364 % if ischar(options.study_dir) && options.save_data_flag==0 0365 % options.save_data_flag=1; 0366 % if options.verbose_flag 0367 % fprintf('Setting ''save_data_flag'' to 1 for storing results in study_dir: %s.\n',options.study_dir); 0368 % end 0369 % end 0370 0371 % Convert matlab solver options from key/value to struct using odeset if 0372 % necessary. 0373 if iscell(options.matlab_solver_options) 0374 options.matlab_solver_options = odeset(options.matlab_solver_options{:}); 0375 end 0376 0377 %% 0.4 Non-Batch checks. 0378 if isempty(options.sim_id) % not in part of a batch sim 0379 if options.optimize_big_vary 0380 options.cluster_flag = 1; 0381 options.qsub_mode = 'array'; 0382 options.compile_flag = 1; 0383 options.downsample_factor = max(1/options.dt, options.downsample_factor); % at most 1000Hz sampling 0384 options.one_solve_file_flag = 1; 0385 options.sims_per_job = 2; 0386 end 0387 0388 % check for one_solve_file_flag 0389 if options.one_solve_file_flag && ~options.cluster_flag 0390 % One file flag only for cluster 0391 fprintf('Since cluster_flag==0, setting options.one_solve_file_flag=0\n') 0392 options.one_solve_file_flag = 0; 0393 % TODO: this is a temp setting until iss_90 is fully implemented 0394 end 0395 0396 if options.one_solve_file_flag && ~options.overwrite_flag 0397 % One file flag will overwrite 0398 fprintf('Since one_solve_file_flag==1, setting options.overwrite_flag=1\n') 0399 options.overwrite_flag = 1; 0400 % TODO: this is a temp setting until iss_90 is fully implemented 0401 end 0402 0403 if options.one_solve_file_flag && ~strcmp(options.qsub_mode, 'array') 0404 % One file flag needs array mode 0405 fprintf('Since one_solve_file_flag==1, setting options.qsub_mode=''array''\n') 0406 options.qsub_mode = 'array'; 0407 % TODO: this is a temp setting until iss_90 is fully implemented 0408 end 0409 0410 if options.one_solve_file_flag && options.parallel_flag 0411 % One file flag can't do parallel_flag 0412 fprintf('Since one_solve_file_flag==1, setting options.parallel_flag=0\n') 0413 options.parallel_flag = 0; 0414 % TODO: this is a temp setting until iss_90 is fully implemented 0415 end 0416 0417 if options.one_solve_file_flag && isa(options.experiment,'function_handle') 0418 error('one_solve_file_flag doesn''t work with experiments.') 0419 end 0420 0421 if options.one_solve_file_flag && ~options.save_parameters_flag 0422 fprintf('Since one_solve_file_flag==1, setting options.save_parameters_flag=1\n') 0423 options.save_parameters_flag = 1; 0424 % TODO: this is a temp setting until iss_90 is fully implemented 0425 end 0426 end % isempty(options.sim_id) 0427 0428 %% 0.5 Prepare analysis functions and options. 0429 if ~isempty(options.analysis_functions) 0430 if ~iscell(options.analysis_functions) 0431 % convert function handle into cell array of function handles 0432 options.analysis_functions={options.analysis_functions}; 0433 end 0434 0435 if any(~cellfun(@(x)isa(x,'function_handle'),options.analysis_functions)) 0436 error('at least one analysis function was not provided as a function handle.'); 0437 end 0438 0439 if isempty(options.analysis_options) 0440 % convert to empty option cell array 0441 options.analysis_options={}; 0442 end 0443 0444 if ~iscell(options.analysis_options) 0445 error('''analysis_options'' must be a cell array of options or option cell arrays'); 0446 end 0447 0448 % force to be a cell array of option cell arrays 0449 if isempty(options.analysis_options) || ischar(options.analysis_options{1}) % first element is an option 0450 options.analysis_options={options.analysis_options}; 0451 end 0452 0453 % make sure there is one option cell array per analysis function 0454 if length(options.analysis_options)==1 && length(options.analysis_functions)>1 0455 % copy options for each analysis function 0456 options.analysis_options=repmat(options.analysis_options,[1 length(options.analysis_functions)]); 0457 elseif length(options.analysis_options) ~= length(options.analysis_functions) 0458 error('there must be one option cell array per analysis function.'); 0459 end 0460 0461 end 0462 0463 %% 0.6 Prepare plot functions and options. 0464 if ~isempty(options.plot_functions) 0465 if ~iscell(options.plot_functions) 0466 % convert function handle into cell array of function handles 0467 options.plot_functions={options.plot_functions}; 0468 end 0469 0470 if any(~cellfun(@(x)isa(x,'function_handle'),options.plot_functions)) 0471 error('at least one plot function was not provided as a function handle.'); 0472 end 0473 0474 if isempty(options.plot_options) 0475 % convert to empty option cell array 0476 options.plot_options={}; 0477 end 0478 0479 if ~iscell(options.plot_options) 0480 error('''plot_options'' must be a cell array of options or option cell arrays'); 0481 end 0482 0483 % force to be a cell array of option cell arrays 0484 if isempty(options.plot_options) || ischar(options.plot_options{1}) % first element is an option 0485 options.plot_options={options.plot_options}; 0486 end 0487 0488 % make sure there is one option cell array per plot function 0489 if length(options.plot_options)==1 && length(options.plot_functions)>1 0490 % copy options for each plot function 0491 options.plot_options=repmat(options.plot_options,[1 length(options.plot_functions)]); 0492 elseif length(options.plot_options) ~= length(options.plot_functions) 0493 error('there must be one option cell array per plot function.'); 0494 end 0495 end 0496 0497 %% 1.0 Prepare model and study structures for simulation. 0498 0499 % Handle special case of input equations with vary() statement. 0500 [model,options] = extract_vary_statement(model,options); 0501 0502 % Check/standardize model. 0503 model = dsCheckModel(model, varargin{:}); % handles conversion when input is a string w/ equations or a DynaSim specification structure 0504 0505 %% 1.1 Apply modifications before simulation and optional further variation across simulations. 0506 if ~isempty(options.modifications) 0507 [model,options.modifications] = dsApplyModifications(model,options.modifications, varargin{:}); 0508 end 0509 0510 %% 1.2 Incorporate user-supplied initial conditions. 0511 if ~isempty(options.ic) 0512 if isstruct(options.ic) 0513 % TODO: create subfunc that converts numeric ICs to strings for 0514 % consistency (call here and in ProcessNumericICs <-- use code already there) 0515 % user provided structure with ICs 0516 warning('off','catstruct:DuplicatesFound'); 0517 model.ICs=catstruct(model.ICs,options.ic); 0518 elseif isnumeric(options.ic) 0519 % user provided numeric array with one value per state variable per 0520 % cell with state variables ordered according to model.state_variables. 0521 model.ICs=ProcessNumericICs; 0522 end 0523 end 0524 0525 % expand set of things to vary across simulations 0526 if ~isempty(options.vary) 0527 modifications_set=dsVary2Modifications(options.vary,model); % Note: options.vary can also be itself a modifications set. 0528 else 0529 modifications_set={[]}; 0530 end 0531 0532 % check for one_solve_file_flag 0533 if options.one_solve_file_flag && is_varied_mech_list() 0534 % Can't vary mechs if using 1 file mode 0535 error('Can''t vary mechanism_list if using one_solve_file_flag') 0536 0537 % TODO: this is a temp setting until iss_90 is fully implemented 0538 end 0539 0540 %% 1.3 Handle parallel simulations. 0541 %% 1.3.1 Manage cluster computing. 0542 % whether to write jobs for distributed processing on cluster 0543 if options.cluster_flag 0544 % add to model any parameters in 'vary' not explicit in current model 0545 % approach: use dsApplyModifications(), it does that automatically 0546 for i = 1:length(modifications_set) 0547 if ~isempty(modifications_set{i}) && ~strcmp(modifications_set{i}{2},'mechanism_list') && ~strcmp(modifications_set{i}{2},'equations') 0548 model = dsApplyModifications(model,modifications_set{i}, varargin{:}); 0549 break 0550 end 0551 end 0552 keyvals = dsOptions2Keyval(options); 0553 studyinfo = dsCreateBatch(model,modifications_set,'simulator_options',options,'process_id',options.sim_id,keyvals{:}); 0554 0555 if options.one_solve_file_flag 0556 % copy params.mat from project_dir to batchdirs 0557 param_file_path = fullfile(options.project_dir, options.study_dir, 'solve','params.mat'); 0558 [~,home]=system('echo $HOME'); 0559 batch_dir = fullfile(strtrim(home),'batchdirs',options.study_dir); 0560 batch_param_file_path = fullfile(batch_dir,'params.mat'); 0561 [success,msg]=copyfile(param_file_path, batch_param_file_path); 0562 end 0563 0564 %if options.overwrite_flag==0 0565 % check status of study 0566 % [~,s]=dsMonitorStudy(studyinfo.study_dir,'verbose_flag',0,'process_id',options.sim_id); 0567 % if s==1 % study finished 0568 % if options.verbose_flag 0569 % fprintf('Study already finished. Importing data...\n'); 0570 % end 0571 % studyinfo=dsCheckStudyinfo(studyinfo.study_dir,'process_id',options.sim_id); 0572 % data=dsImport(studyinfo,'process_id',options.sim_id); 0573 % end 0574 %end 0575 0576 %% auto_gen_test_data_flag argout 0577 if options.auto_gen_test_data_flag 0578 if isfield(data, 'simulator_options') 0579 data= rmfield(data, 'simulator_options'); % specific to this function 0580 end 0581 if ~isempty(studyinfo) 0582 studyinfo = []; % specific to this function 0583 end 0584 0585 argout = {data, studyinfo}; % specific to this function 0586 0587 % file output dir 0588 % if ~isempty(studyinfo) && ~isempty(studyinfo.study_dir) 0589 % dirOut = studyinfo.study_dir; 0590 % else 0591 dirOut = options.study_dir; 0592 % end 0593 0594 removeStudyinfo(); % remove studyinfo file 0595 renameMexFilesForUnitTesting(); % rename mex files for unit testing 0596 0597 dsUnitSaveAutoGenTestDir(argin, argout, [], dirOut); 0598 end 0599 0600 %% unit test 0601 if options.unit_test_flag 0602 % remove fields that cause issues in unit testing 0603 if isfield(data, 'simulator_options') 0604 data= rmfield(data, 'simulator_options'); % specific to this function 0605 end 0606 if ~isempty(studyinfo) 0607 studyinfo = []; 0608 end 0609 0610 removeStudyinfo(); % remove studyinfo file 0611 renameMexFilesForUnitTesting(); % rename mex files for unit testing 0612 end 0613 0614 return; 0615 end 0616 0617 %% 1.3.2 Manage parallel computing on local machine. 0618 % TODO: debug local parallel sims, doesn't seem to be working right... 0619 % (however SCC cluster+parallel works) 0620 if options.parallel_flag 0621 % prepare studyinfo 0622 [studyinfo,options]=dsSetupStudy(model,'simulator_options',options,'modifications_set',modifications_set); 0623 0624 if options.compile_flag 0625 % Ensure mex_dir is absolute 0626 if (ispc && options.mex_dir(2) == ':') || (~ispc && options.mex_dir(1) == filesep) 0627 relMexPath = false; 0628 else 0629 relMexPath = true; 0630 end 0631 0632 if relMexPath % then make absolute by prepending study_dir 0633 options.mex_dir = fullfile(getAbsolutePath(studyinfo.study_dir), options.mex_dir); 0634 end 0635 0636 % Create mex_dir if it does not yet exist 0637 if ~exist(options.mex_dir,'dir') && ~options.cluster_flag && options.mex_dir_flag 0638 mkdir(options.mex_dir); 0639 end 0640 end 0641 0642 % prepare options 0643 options_temp = rmfield(options,{'vary','modifications','solve_file','parallel_flag','studyinfo','in_parfor_loop_flag'}); 0644 keyvals=dsOptions2Keyval(options_temp); 0645 0646 keyvals{find(strcmp(keyvals, 'auto_gen_test_data_flag'))+1} = 0; 0647 if options.auto_gen_test_data_flag || options.unit_test_flag 0648 keyvals{find(strcmp(keyvals, 'unit_test_flag'))+1} = 1; % prevents time-stamped outputs 0649 end 0650 0651 % run embarrassingly-parallel simulations 0652 0653 % Previous parallel code would overwrite the same params.mat file on each 0654 % parallel iteration, resulting in the same parameters being used for all 0655 % simulations. 0656 0657 % % List any core files - these should be deleted, as they are huge (debug) 0658 % system (['ls ' fullfile(options.study_dir,'output*')],'-echo'); 0659 % system('find * -name "core*"','-echo'); 0660 0661 clear data 0662 0663 % note that parfor currently acts just as a regular for in Octave 0664 parfor sim=1:length(modifications_set) 0665 data(sim)=dsSimulate(model, 'modifications', modifications_set{sim}, keyvals{:},... 0666 'studyinfo', studyinfo, 'sim_id',sim, 'in_parfor_loop_flag', 1); % My modification; now specifies a separate study directory for each sim. 0667 %disp(sim); 0668 end 0669 0670 % Clean up files leftover from sim 0671 % Unfortunately we can't remove the folders due to locked .nfs files. 0672 % Need to do this manually later... 0673 0674 if strcmp(reportUI,'matlab') % parfor currently acts just as a regular for in Octave 0675 % Delete any core files in parent directory 0676 delete(fullfile(options.study_dir,'core*')); 0677 0678 % Verify all core files are deleted 0679 [~,result] = system('find * -name "core*"','-echo'); 0680 if ~isempty(result); fprintf(strcat(result,'\n')); warning('Core files found. Consider deleting to free up space'); end 0681 end 0682 0683 % TODO: sort data sets by things varied in modifications_set 0684 % TODO: Figure out how to delete locked .nfs files 0685 0686 if options.verbose_flag 0687 fprintf('\nParallel simulations complete.\n\n') 0688 end 0689 0690 %% auto_gen_test_data_flag argout 0691 if options.auto_gen_test_data_flag 0692 if isfield(data, 'simulator_options') 0693 data= rmfield(data, 'simulator_options'); % specific to this function 0694 end 0695 if ~isempty(studyinfo) 0696 studyinfo = []; % specific to this function 0697 end 0698 0699 argout = {data, studyinfo}; % specific to this function 0700 0701 % file output dir 0702 % if ~isempty(studyinfo) && ~isempty(studyinfo.study_dir) 0703 % dirOut = studyinfo.study_dir; 0704 % else 0705 dirOut = options.study_dir; 0706 % end 0707 0708 removeStudyinfo(); % remove studyinfo file 0709 renameMexFilesForUnitTesting(); % rename mex files for unit testing 0710 0711 dsUnitSaveAutoGenTestDir(argin, argout, [], dirOut); 0712 end 0713 0714 %% unit test 0715 if options.unit_test_flag 0716 % remove fields that cause issues in unit testing 0717 if isfield(data, 'simulator_options') 0718 data= rmfield(data, 'simulator_options'); % specific to this function 0719 end 0720 if ~isempty(studyinfo) 0721 studyinfo = []; 0722 end 0723 0724 removeStudyinfo(); % remove studyinfo file 0725 renameMexFilesForUnitTesting(); % rename mex files for unit testing 0726 end 0727 0728 return 0729 end %parallel_flag 0730 0731 %% 1.4 prepare study_dir and studyinfo if saving data 0732 if isempty(options.studyinfo) 0733 [studyinfo,options] = dsSetupStudy(model,'modifications_set',modifications_set,'simulator_options',options,'process_id',options.sim_id); 0734 else % in parfor loop and/or cluster job 0735 studyinfo = options.studyinfo; 0736 end 0737 0738 if options.compile_flag 0739 % Ensure mex_dir is absolute 0740 if (ispc && options.mex_dir(2) == ':') || (~ispc && options.mex_dir(1) == filesep) 0741 relMexPath = false; 0742 else 0743 relMexPath = true; 0744 end 0745 0746 if relMexPath % then make absolute by prepending study_dir 0747 if ~isempty(studyinfo) && ~isempty(studyinfo.study_dir) 0748 options.mex_dir = fullfile(getAbsolutePath(studyinfo.study_dir), options.mex_dir); 0749 else 0750 0751 options.mex_dir = fullfile(getAbsolutePath(options.study_dir), options.mex_dir); 0752 end 0753 end 0754 0755 % Create mex_dir if it does not yet exist 0756 if ~exist(options.mex_dir,'dir') && ~options.cluster_flag && options.mex_dir_flag 0757 mkdir(options.mex_dir); 0758 end 0759 end 0760 0761 % put solver blocks in try statement to catch and handle errors/cleanup 0762 cwd=pwd; % record current directory 0763 0764 data_index=0; 0765 tmpdata = []; 0766 0767 if ~options.debug_flag 0768 try 0769 tryFn(nargout) 0770 catch err % error handling 0771 if options.compile_flag && ~isempty(options.solve_file) && ~options.one_solve_file_flag 0772 0773 if options.verbose_flag 0774 fprintf('Removing failed compiled solve file: %s\n',options.solve_file); 0775 end 0776 0777 delete([options.solve_file '*']); 0778 end 0779 0780 displayError(err); 0781 0782 % update studyinfo 0783 if options.save_data_flag && ~options.one_solve_file_flag 0784 studyinfo=dsUpdateStudy(studyinfo.study_dir,'process_id',sim_id,'status','failed','verbose_flag',options.verbose_flag); 0785 data=studyinfo; 0786 end 0787 cleanup('error'); 0788 0789 rethrow(err) 0790 end 0791 else 0792 tryFn(nargout) % outside of try block if options.debug_flag 0793 end 0794 0795 if options.verbose_flag 0796 fprintf('\nSimulations complete.\n\n') 0797 end 0798 0799 if ~options.in_parfor_loop_flag % if not inside of parfor loop 0800 %% auto_gen_test_data_flag argout 0801 if options.auto_gen_test_data_flag 0802 if isfield(data, 'simulator_options') 0803 data = rmfield(data, 'simulator_options'); % specific to this function 0804 end 0805 if ~isempty(studyinfo) 0806 studyinfo = []; % specific to this function 0807 end 0808 0809 argout = {data, studyinfo}; % specific to this function 0810 0811 % file output dir 0812 % if ~isempty(studyinfo) && ~isempty(studyinfo.study_dir) 0813 % dirOut = studyinfo.study_dir; 0814 % else 0815 dirOut = options.study_dir; 0816 % end 0817 0818 removeStudyinfo(); % remove studyinfo file 0819 renameMexFilesForUnitTesting(); % rename mex files for unit testing 0820 0821 dsUnitSaveAutoGenTestDir(argin, argout, [], dirOut); 0822 end 0823 0824 %% unit test 0825 if options.unit_test_flag 0826 % remove fields that cause issues in unit testing 0827 if isfield(data, 'simulator_options') 0828 data= rmfield(data, 'simulator_options'); % specific to this function 0829 end 0830 if ~isempty(studyinfo) 0831 studyinfo = []; 0832 end 0833 0834 removeStudyinfo(); % remove studyinfo file 0835 renameMexFilesForUnitTesting(); % rename mex files for unit testing 0836 end 0837 end % in_parfor_loop_flag 0838 0839 % --------------------------------------------- 0840 % TODO: 0841 % - create helper function that handles log files (creation, standardized format,...) 0842 % --------------------------------------------- 0843 0844 0845 0846 %% ------------------------- 0847 % NESTED FUNCTIONS 0848 % ------------------------- 0849 function tryFn(nargoutmain) 0850 % functions: outputs: 0851 % dsWriteDynaSimSolver m-file for DynaSim solver 0852 % dsWriteMatlabSolver m-file for Matlab solver (including @odefun) 0853 % dsPrepareMEX mex-file for m-file 0854 0855 %% 1.5 loop over simulations, possibly varying things 0856 base_model=model; 0857 0858 for sim=1:length(modifications_set) 0859 if ~strcmp(pwd,cwd) % move back to original directory before potentially regenerating to make sure the model files used are the same 0860 cd(cwd); 0861 end 0862 0863 % get index for this simulation 0864 if ~isempty(options.sim_id) 0865 sim_ind=find([studyinfo.simulations.sim_id]==options.sim_id); 0866 sim_id=options.sim_id; 0867 else 0868 sim_ind=sim; 0869 sim_id=sim; 0870 end 0871 0872 if options.save_data_flag 0873 % check if output data already exists. load if so and skip simulation 0874 data_file=studyinfo.simulations(sim_ind).data_file; 0875 if exist(data_file,'file') && ~options.overwrite_flag 0876 if 1%options.verbose_flag 0877 % note: this is important, should always display 0878 fprintf('Loading data from %s\n',data_file); 0879 end 0880 tmpdata=dsImport(data_file,'process_id',sim_id); 0881 update_data; % concatenate data structures across simulations 0882 continue; % skip to next simulation 0883 end 0884 end 0885 0886 % apply modifications for this point in search space 0887 if ~isempty(modifications_set{sim}) 0888 model=dsApplyModifications(base_model,modifications_set{sim}, varargin{:}); 0889 end 0890 0891 % update studyinfo 0892 if options.save_data_flag 0893 %studyinfo=dsUpdateStudy(studyinfo.study_dir,'process_id',sim_id,'status','started','model',model,'simulator_options',options,'verbose_flag',options.verbose_flag); 0894 end 0895 0896 %% Experiment 0897 if isa(options.experiment,'function_handle') 0898 % EXPERIMENT (wrapping around a set of simulations) 0899 if options.cluster_flag && options.compile_flag 0900 warning('compiled solver is not available for experiments on the cluster. Simulation will be run in Matlab.'); 0901 end 0902 0903 % from varargin... 0904 % remove 'experiment', 'modifications', 'vary', 'cluster_flag' to avoid undesired recursive action in experiment function 0905 % remove 'save_data_flag' to prevent individual simulations from being saved during experiment 0906 keyvals=dsRemoveKeyval(varargin,{'experiment','cluster_flag','vary','modifications','save_data_flag'}); 0907 0908 if ~isempty(options.experiment_options) 0909 % user-supplied experiment options override any found in dsSimulate options 0910 keyvals=dsRemoveKeyval(keyvals,options.experiment_options(1:2:end)); 0911 keyvals=cat(2,keyvals,options.experiment_options); 0912 end 0913 0914 tmpdata=feval(options.experiment,model,keyvals{:}); 0915 else 0916 %% NOT AN EXPERIMENT (single simulation) 0917 %% 2.0 prepare solver function (solve_ode.m/mex) 0918 % - Matlab solver: create @odefun with vectorized state variables 0919 % - DynaSim solver: write solve_ode.m and params.mat (based on dnsimulator()) 0920 % check if model solver needs to be created 0921 % (i.e., if is first simulation or a search space varying mechanism list) 0922 0923 if sim==1 || ( ~isempty(modifications_set{1}) && is_varied_mech_list() ) 0924 % prepare file that solves the model system 0925 if isempty(options.solve_file) || (~exist(options.solve_file,'file') &&... 0926 ~exist([options.solve_file '.mexa64'],'file') &&... 0927 ~exist([options.solve_file '.mexa32'],'file') &&... 0928 ~exist([options.solve_file '.mexmaci64'],'file')) 0929 options.solve_file = dsGetSolveFile(model,studyinfo,options); % store name of solver file in options struct 0930 end 0931 0932 % TODO: consider providing better support for studies that produce different m-files per sim (e.g., varying mechanism_list) 0933 0934 if options.verbose_flag 0935 fprintf('\nSIMULATING MODEL:\n'); 0936 fprintf('Solving system using %s\n',options.solve_file); 0937 end 0938 else 0939 % use previous solve_file 0940 end 0941 [fpath,fname,fext]=fileparts(options.solve_file); 0942 0943 %% 3.0 integrate model with solver of choice and prepare output data 0944 % - matlab solver: solve @odefun with feval and solver_options 0945 % - DynaSim solver: run solve_ode.m or create/run MEX 0946 % move to directory with solver file 0947 0948 if options.verbose_flag 0949 fprintf('Changing directory to %s\n',fpath); 0950 end 0951 0952 cd(fpath); 0953 0954 % save parameters there 0955 warning('off','catstruct:DuplicatesFound'); 0956 p = catstruct(dsCheckSolverOptions(options),model.parameters); 0957 0958 if matlabSolverBool 0959 % add IC to p for use in matlab solver 0960 if isempty(options.ic) 0961 [~, p.ic]=dsDynasim2odefun( dsPropagateParameters( dsPropagateFunctions(model, varargin{:}), varargin{:} ), 'odefun_output','func_body', varargin{:} ); 0962 else 0963 p.ic = options.ic; 0964 end 0965 0966 % add matlab_solver_options to p 0967 if ~isempty(options.matlab_solver_options) 0968 p.matlab_solver_options = options.matlab_solver_options; 0969 end 0970 end 0971 0972 param_file = fullfile(fpath,'params.mat'); 0973 if options.verbose_flag 0974 fprintf('Saving model parameters: %s\n',param_file); 0975 end 0976 %pause(.01); 0977 0978 %% Solve System 0979 if options.disk_flag % ### data stored on disk during simulation ### 0980 sim_start_time=tic; 0981 if ~options.one_solve_file_flag 0982 save(param_file,'p','-v7'); % save params immediately before solving 0983 end 0984 csv_data_file=feval(fname); % returns name of file storing the simulated data 0985 duration=toc(sim_start_time); 0986 0987 if nargout>0 || options.save_data_flag 0988 tmpdata=dsImport(csv_data_file,'process_id',sim_id); % eg, data.csv 0989 end 0990 else % ### data stored in memory during simulation ### 0991 % create list of output variables to capture 0992 output_variables=cat(2,'time',model.state_variables); 0993 0994 if ~isempty(model.monitors) 0995 output_variables=cat(2,output_variables,fieldnames(model.monitors)'); 0996 end 0997 0998 if ~isempty(model.fixed_variables) 0999 fields=fieldnames(model.fixed_variables)'; 1000 output_variables=cat(2,output_variables,fields); 1001 num_fixed_variables=length(fields); 1002 else 1003 num_fixed_variables=0; 1004 end 1005 1006 % run simulation 1007 if options.verbose_flag 1008 fprintf('\nRunning simulation %g/%g (solver=''%s'', dt=%g, tspan=[%g %g]) ...\n',sim,length(modifications_set),options.solver,options.dt,options.tspan); 1009 end 1010 sim_start_time=tic; 1011 1012 outputs=cell(1,length(output_variables)); % preallocate for PCT compatibility 1013 1014 if ~options.one_solve_file_flag 1015 save(param_file,'p','-v7'); % save params immediately before solving 1016 end 1017 1018 % feval solve file 1019 if ~options.one_solve_file_flag 1020 [outputs{1:length(output_variables)}]=feval(fname); 1021 else 1022 % pass sim_id for slicing params 1023 [outputs{1:length(output_variables)}]=feval(fname, sim_id); 1024 end 1025 1026 duration=toc(sim_start_time); 1027 1028 1029 % prepare DynaSim data structure 1030 % organize simulated data in data structure (move time to last) 1031 tmpdata.labels=output_variables([2:length(output_variables)-num_fixed_variables 1]); 1032 for i=1:length(output_variables) 1033 if ~isempty(model.fixed_variables) && isfield(model.fixed_variables,output_variables{i}) 1034 % store fixed variables in model substructure 1035 model.fixed_variables.(output_variables{i})=outputs{i}; 1036 else 1037 % store state variables and monitors as data fields 1038 tmpdata.(output_variables{i})=outputs{i}; 1039 end 1040 1041 outputs{i}=[]; % clear assigned outputs from memory 1042 end 1043 end 1044 1045 if options.verbose_flag 1046 fprintf('Elapsed time: %g seconds.\n',duration); 1047 end 1048 1049 % add metadata to tmpdata 1050 tmpdata.simulator_options=options; % store simulator controls 1051 1052 if options.store_model_flag==1 % optionally store the simulated model 1053 tmpdata.model=model; 1054 end 1055 end 1056 1057 tmpdata = dsModifications2Vary(tmpdata,options.modifications,options,modifications_set,sim); 1058 1059 if (options.auto_gen_test_data_flag || options.unit_test_flag) && isfield(tmpdata, 'simulator_options') 1060 tmpdata= rmfield(tmpdata, 'simulator_options'); 1061 end 1062 1063 % save single data set and update studyinfo 1064 if options.save_data_flag 1065 dsExportData(tmpdata,'filename',data_file,'format','mat','matCompatibility_flag',options.matCompatibility_flag,'verbose_flag',options.verbose_flag); 1066 %studyinfo=dsUpdateStudy(studyinfo.study_dir,'process_id',sim_id,'status','finished','duration',duration,'solve_file',options.solve_file,'email',options.email,'verbose_flag',options.verbose_flag,'model',model,'simulator_options',options); 1067 end 1068 1069 % do post-simulation analysis and plotting 1070 if ~isempty(options.analysis_functions) || ~isempty(options.plot_functions) 1071 if options.save_data_flag || options.save_results_flag 1072 % do analysis and plotting while saving results 1073 siminfo=studyinfo.simulations(sim_ind); 1074 for f=1:length(siminfo.result_functions) 1075 tmpresult=dsAnalyze(tmpdata,siminfo.result_functions{f},'result_file',siminfo.result_files{f},'save_data_flag',1,'save_results_flag',1,siminfo.result_options{f}{:},'parallel_flag',options.parallel_flag); 1076 1077 % since the plots are saved, close all generated figures 1078 if all(ishandle(tmpresult)) 1079 close(tmpresult); 1080 end 1081 end 1082 else 1083 % do analysis and plotting without saving results 1084 if ~isempty(options.analysis_functions) && nargoutmain > 2 1085 for f=1:length(options.analysis_functions) 1086 tmpresult=dsAnalyze(tmpdata,options.analysis_functions{f},'result_file',[],'save_data_flag',0,'save_results_flag',options.save_results_flag,options.analysis_options{f}{:},'parallel_flag',options.parallel_flag); 1087 end 1088 end 1089 1090 if ~isempty(options.plot_functions) 1091 for f=1:length(options.plot_functions) 1092 dsAnalyze(tmpdata,options.plot_functions{f},'result_file',[],'save_data_flag',0,'save_results_flag',options.save_results_flag,options.plot_options{f}{:},'parallel_flag',options.parallel_flag); 1093 end 1094 end 1095 end 1096 end 1097 1098 if nargoutmain>0 1099 update_data; % concatenate data structures across simulations 1100 end 1101 1102 if nargoutmain>2 1103 update_result; 1104 end 1105 end % end loop over sims 1106 1107 cleanup('success'); 1108 end %tryfn 1109 1110 function update_data 1111 % store tmpdata 1112 if sim==1 1113 % replicate first data set as preallocation for all 1114 data=repmat(tmpdata,[1 length(modifications_set)]); 1115 data_index=length(tmpdata); 1116 else 1117 inds=data_index+(1:length(tmpdata)); % support multiple data sets returned by experiments 1118 data(inds)=tmpdata; 1119 data_index=inds(end); 1120 end 1121 end 1122 1123 function update_result 1124 % store tmpdata 1125 if sim==1 1126 % replicate first data set as preallocation for all 1127 result=repmat(tmpresult,[1 length(modifications_set)]); 1128 result_index=length(tmpresult); 1129 else 1130 inds=data_index+(1:length(tmpresult)); % support multiple data sets returned by experiments 1131 result(inds)=tmpresult; 1132 result_index=inds(end); 1133 end 1134 end 1135 1136 function cleanup(status) 1137 % remove temporary files and optionally store info for debugging 1138 % ... 1139 % return to original directory 1140 if options.verbose_flag 1141 fprintf('Changing directory to %s\n',cwd); 1142 end 1143 cd(cwd); 1144 switch status 1145 case 'success' 1146 % ... 1147 % TODO: consider removing solve folder if nothing being saved 1148 case 'error' 1149 % ... error logs 1150 end 1151 end 1152 1153 function all_ICs=ProcessNumericICs 1154 % first, figure out how many IC values we need (i.e., how many state 1155 % variables we need across all cells). 1156 var_names=model.state_variables; 1157 [nvals_per_var,monitor_counts]=dsGetOutputCounts(model); 1158 num_state_variables=sum(nvals_per_var); 1159 % check that the correct number of IC values was provided 1160 if length(options.ic)~=num_state_variables 1161 error('incorrect number of initial conditions. %g values are needed for %g state variables across %g cells',num_state_variables,length(model.state_variables),sum(pop_sizes)); 1162 end 1163 % organize user-supplied ICs into array for each state variable (assume 1164 cnt=0; all_ICs=[]; 1165 for i=1:length(var_names) 1166 ICs=options.ic(cnt+(1:nvals_per_var(i))); 1167 % store ICs as string for writing solve_ode and consistent evaluation 1168 all_ICs.(var_names{i})=sprintf('[%s]',num2str(ICs)); 1169 cnt=cnt+nvals_per_var(i); 1170 end 1171 end 1172 1173 function logicalOut = is_varied_mech_list() 1174 if ~isempty(modifications_set{1}) 1175 logicalOut = any(cellfun(@(x) strcmp(x{2},'mechanism_list'),modifications_set)); 1176 else 1177 logicalOut = false; 1178 end 1179 end 1180 1181 function removeStudyinfo() 1182 % Problem: studyinfo file has many timestamps and absolute paths 1183 % Solution: remove studyinfo file 1184 studyinfoFile = fullfile(options.study_dir, 'studyinfo.mat'); 1185 if exist(studyinfoFile, 'file') 1186 delete(studyinfoFile) 1187 end 1188 end 1189 1190 function renameMexFilesForUnitTesting() 1191 % Problem: different systems make different mex files, and mex files are 1192 % different from simulation to simulation on same computer 1193 % Solution: rename extension to a general one, mex4unittest, and skip these 1194 % files when testing 1195 1196 studyDirFiles = rls(options.study_dir); 1197 studyDirFiles = studyDirFiles(~cellfun(@isempty, strfind(studyDirFiles, '.mex'))); 1198 for k = 1:length(studyDirFiles) 1199 thisFile = studyDirFiles{k}; 1200 renamedFile = regexprep(thisFile, '\.mex.+', '.mex4unittest'); 1201 movefile(thisFile, renamedFile); 1202 end 1203 end 1204 1205 end %main fn 1206 1207 1208 %% Subfunctions 1209 function [model,options]=extract_vary_statement(model,options) 1210 % Purpose: extract vary statement, remove from model, and set options.vary 1211 if ischar(model) && any(regexp(model,';\s*vary\(.*\)','once')) 1212 % extract vary statement 1213 str=regexp(model,';\s*(vary\(.*\);?)','tokens','once'); 1214 1215 % remove from model 1216 model=strrep(model,str{1},''); 1217 1218 % set options 1219 var=regexp(str{1},'\((.*)=','tokens','once'); % variable 1220 val=regexp(str{1},'=(.*)\)','tokens','once'); % values 1221 options.vary={'pop1',var{1},eval(val{1})}; 1222 end 1223 end 1224 1225 function options = backward_compatibility(options) 1226 % option_names: (old_name, new_name; ...} 1227 option_names = {... 1228 'override','modifications'; 1229 'timelimits','tspan'; 1230 'time_limits','tspan'; 1231 'IC','ic'; 1232 'verbose','verbose_flag'; 1233 'SOLVER','solver'; 1234 'nofunctions','reduce_function_calls_flag'; 1235 'dsfact','downsample_factor'; 1236 'memlimit','memory_limit'; 1237 }; 1238 1239 if any(ismember(option_names(:,1),options(1:2:end))) 1240 for i=1:size(option_names,1) 1241 % check if any options have this old name 1242 if ismember(option_names{i,1},options(1:2:end)) 1243 ind=find(ismember(options(1:2:end),option_names{i,1})); 1244 1245 % replace old option name by new option name 1246 options{2*ind-1}=option_names{i,2}; 1247 end 1248 end 1249 end 1250 1251 end