dsWriteDynaSimSolver

WRITEDYNASIMSOLVER - write m-file that numerically inteegrates the model

 Usage:
   solver_file=dsWriteDynaSimSolver(model,varargin)

 Inputs:
   - model: DynaSim model structure (see dsGenerateModel)
   - options:
     'solver'     : solver for numerical integration (see dsGetSolveFile)
                    {'euler'/'rk1', 'rk2'/'modified_euler', 'rungekutta'/'rk'/'rk4'} (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)
                          - Note: 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. This overrides definition in model structure (default:
                    all zeros)
     'random_seed': seed for random number generator (default: 'shuffle', set
                    randomly) (usage: rng(options.random_seed))
     'disk_flag'  : whether to write to disk during simulation instead of
                    storing in memory {0 or 1} (default: 0)
     'sparse_flag' : whether to convert numeric fixed variables to sparse matrices {0 or 1} (default: 0)

 Outputs:
   - solver_file (e.g., solve_ode.m): function that numerically integrates the model

 Examples:
   - Example 1: test solver production, display function in standard output
       model=dsGenerateModel; % test model
       without writing anything to disk:
       dsWriteDynaSimSolver(model,'fileID',1,'save_parameters_flag',0,'solver','rk4');
       dsWriteDynaSimSolver(model,'fileID',1,'save_parameters_flag',1,'solver','rk4');
       model=dsPropagateFunctions(model);
       dsWriteDynaSimSolver(model,'fileID',1,'save_parameters_flag',0,'solver','rk4');
       dsWriteDynaSimSolver(model,'fileID',1,'save_parameters_flag',1,'solver','rk4');

   - Example 2: real-time downsampling
       dsWriteDynaSimSolver(model,'downsample_factor',10,'fileID',1,'solver','rk4');

   - Example 3: real-time writing data to disk (reduce memory demand)
       dsWriteDynaSimSolver(model,'disk_flag',1,'fileID',1,'solver','rk4');

   - Example 4: maximally conserve memory: downsample and write to disk
       dsWriteDynaSimSolver(model,'disk_flag',1,'downsample_factor',10,'fileID',1,'solver','rk4');

 More Examples:
     dsWriteDynaSimSolver(model,'solver','euler');
     dsWriteDynaSimSolver(model,'solver','rk2');
     dsWriteDynaSimSolver(model,'solver','rk4');

     model=dsGenerateModel; % test model
     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',0,'reduce_function_calls_flag',0,'solver','rk4','filename','solve_ode.m'); v=solve_ode; plot(v); toc
     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',1,'reduce_function_calls_flag',0,'solver','rk4','filename','solve_ode.m'); v=solve_ode; plot(v); toc
     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',0,'reduce_function_calls_flag',1,'solver','rk4','filename','solve_ode.m'); v=solve_ode; plot(v); toc
     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',1,'reduce_function_calls_flag',1,'solver','rk4','filename','solve_ode.m'); v=solve_ode; plot(v); toc
     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',1,'reduce_function_calls_flag',1,'solver','rk2','filename','solve_ode.m'); v=solve_ode; plot(v); toc
     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',1,'reduce_function_calls_flag',1,'solver','rk1','filename','solve_ode.m'); v=solve_ode; plot(v); toc
     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',1,'reduce_function_calls_flag',1,'solver','rk1','filename','solve_ode.m','downsample_factor',10); v=solve_ode; plot(v); toc
     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',1,'reduce_function_calls_flag',1,'solver','rk2','filename','solve_ode.m','dt',.001,'downsample_factor',10); v=solve_ode; plot(v); toc
     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',1,'reduce_function_calls_flag',1,'solver','rk1','filename','solve_ode.m','dt',.001,'downsample_factor',10); v=solve_ode; plot(v); toc
     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',1,'reduce_function_calls_flag',1,'solver','rk1','filename','solve_ode.m','dt',.005,'tspan',[0 200],'downsample_factor',10); v=solve_ode; plot(v); toc

 Dependencies: dsCheckOptions, dsCheckModel

 See also: dsGetSolveFile, dsSimulate, dsWriteMatlabSolver

 Author: Jason Sherfey, PhD <jssherfey@gmail.com>
 Copyright (C) 2016 Jason Sherfey, Boston University, USA

0001 function [outfile,options] = dsWriteDynaSimSolver(model,varargin)
0002 %WRITEDYNASIMSOLVER - write m-file that numerically inteegrates the model
0003 %
0004 % Usage:
0005 %   solver_file=dsWriteDynaSimSolver(model,varargin)
0006 %
0007 % Inputs:
0008 %   - model: DynaSim model structure (see dsGenerateModel)
0009 %   - options:
0010 %     'solver'     : solver for numerical integration (see dsGetSolveFile)
0011 %                    {'euler'/'rk1', 'rk2'/'modified_euler', 'rungekutta'/'rk'/'rk4'} (default: 'rk4')
0012 %     'tspan'      : time limits of simulation [begin,end] (default: [0 100]) [ms]
0013 %                    - Note: units must be consistent with dt and model equations
0014 %     'dt'         : time step used for DynaSim solvers (default: .01) [ms]
0015 %     'downsample_factor': downsampling applied during simulation (default: 1, no downsampling)
0016 %                          - Note: only every downsample_factor-time point is
0017 %                                  stored in memory and/or written to disk
0018 %     'ic'         : numeric array of initial conditions, one value per state
0019 %                    variable. This overrides definition in model structure (default:
0020 %                    all zeros)
0021 %     'random_seed': seed for random number generator (default: 'shuffle', set
0022 %                    randomly) (usage: rng(options.random_seed))
0023 %     'disk_flag'  : whether to write to disk during simulation instead of
0024 %                    storing in memory {0 or 1} (default: 0)
0025 %     'sparse_flag' : whether to convert numeric fixed variables to sparse matrices {0 or 1} (default: 0)
0026 %
0027 % Outputs:
0028 %   - solver_file (e.g., solve_ode.m): function that numerically integrates the model
0029 %
0030 % Examples:
0031 %   - Example 1: test solver production, display function in standard output
0032 %       model=dsGenerateModel; % test model
0033 %       without writing anything to disk:
0034 %       dsWriteDynaSimSolver(model,'fileID',1,'save_parameters_flag',0,'solver','rk4');
0035 %       dsWriteDynaSimSolver(model,'fileID',1,'save_parameters_flag',1,'solver','rk4');
0036 %       model=dsPropagateFunctions(model);
0037 %       dsWriteDynaSimSolver(model,'fileID',1,'save_parameters_flag',0,'solver','rk4');
0038 %       dsWriteDynaSimSolver(model,'fileID',1,'save_parameters_flag',1,'solver','rk4');
0039 %
0040 %   - Example 2: real-time downsampling
0041 %       dsWriteDynaSimSolver(model,'downsample_factor',10,'fileID',1,'solver','rk4');
0042 %
0043 %   - Example 3: real-time writing data to disk (reduce memory demand)
0044 %       dsWriteDynaSimSolver(model,'disk_flag',1,'fileID',1,'solver','rk4');
0045 %
0046 %   - Example 4: maximally conserve memory: downsample and write to disk
0047 %       dsWriteDynaSimSolver(model,'disk_flag',1,'downsample_factor',10,'fileID',1,'solver','rk4');
0048 %
0049 % More Examples:
0050 %     dsWriteDynaSimSolver(model,'solver','euler');
0051 %     dsWriteDynaSimSolver(model,'solver','rk2');
0052 %     dsWriteDynaSimSolver(model,'solver','rk4');
0053 %
0054 %     model=dsGenerateModel; % test model
0055 %     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',0,'reduce_function_calls_flag',0,'solver','rk4','filename','solve_ode.m'); v=solve_ode; plot(v); toc
0056 %     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',1,'reduce_function_calls_flag',0,'solver','rk4','filename','solve_ode.m'); v=solve_ode; plot(v); toc
0057 %     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',0,'reduce_function_calls_flag',1,'solver','rk4','filename','solve_ode.m'); v=solve_ode; plot(v); toc
0058 %     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',1,'reduce_function_calls_flag',1,'solver','rk4','filename','solve_ode.m'); v=solve_ode; plot(v); toc
0059 %     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',1,'reduce_function_calls_flag',1,'solver','rk2','filename','solve_ode.m'); v=solve_ode; plot(v); toc
0060 %     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',1,'reduce_function_calls_flag',1,'solver','rk1','filename','solve_ode.m'); v=solve_ode; plot(v); toc
0061 %     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',1,'reduce_function_calls_flag',1,'solver','rk1','filename','solve_ode.m','downsample_factor',10); v=solve_ode; plot(v); toc
0062 %     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',1,'reduce_function_calls_flag',1,'solver','rk2','filename','solve_ode.m','dt',.001,'downsample_factor',10); v=solve_ode; plot(v); toc
0063 %     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',1,'reduce_function_calls_flag',1,'solver','rk1','filename','solve_ode.m','dt',.001,'downsample_factor',10); v=solve_ode; plot(v); toc
0064 %     tic; dsWriteDynaSimSolver(model,'save_parameters_flag',1,'reduce_function_calls_flag',1,'solver','rk1','filename','solve_ode.m','dt',.005,'tspan',[0 200],'downsample_factor',10); v=solve_ode; plot(v); toc
0065 %
0066 % Dependencies: dsCheckOptions, dsCheckModel
0067 %
0068 % See also: dsGetSolveFile, dsSimulate, dsWriteMatlabSolver
0069 %
0070 % Author: Jason Sherfey, PhD <jssherfey@gmail.com>
0071 % Copyright (C) 2016 Jason Sherfey, Boston University, USA
0072 
0073 % Check inputs
0074 options=dsCheckOptions(varargin,{...
0075   'tspan',[0 100],[],...          % [beg,end] (units must be consistent with dt and equations)
0076   'downsample_factor',1,[],...    % downsampling applied during simulation (only every downsample_factor-time point is stored in memory or written to disk)
0077   'dt',.01,[],...                 % time step used for fixed step DynaSim solvers
0078   'random_seed','shuffle',[],...        % seed for random number generator (usage: rng(random_seed))
0079   'solver','rk4',{'euler','rk1','rk2','rk4','modified_euler','rungekutta','rk'},... % DynaSim solvers
0080   'disk_flag',0,{0,1},...            % whether to write to disk during simulation instead of storing in memory
0081   'reduce_function_calls_flag',1,{0,1},...   % whether to eliminate internal (anonymous) function calls
0082   'save_parameters_flag',1,{0,1},...
0083   'filename',[],[],...         % name of solver file that integrates model
0084   'data_file','data.csv',[],... % name of data file if disk_flag=1
0085   'fileID',1,[],...
0086   'compile_flag',0,{0,1},... % whether to prepare script for being compiled using coder instead of interpreting Matlab
0087   'verbose_flag',1,{0,1},...
0088   'sparse_flag',0,{0,1},...
0089   'one_solve_file_flag',0,{0,1},... % use only 1 solve file of each type, but can't vary mechs yet
0090   'benchmark_flag',0,{0,1},...
0091   },false);
0092 model=dsCheckModel(model, varargin{:});
0093 separator=','; % ',', '\\t'
0094 
0095 %% 1.0 prepare model info
0096 parameter_prefix='p.';
0097 state_variables=model.state_variables;
0098 
0099 % 1.1 eliminate internal (anonymous) function calls from model equations
0100 if options.reduce_function_calls_flag==1
0101   model=dsPropagateFunctions(model, varargin{:});
0102 end
0103 
0104 % 1.1 prepare parameters
0105 if options.save_parameters_flag
0106   % add parameter struct prefix to parameters in model equations
0107   model=dsPropagateParameters(model,'action','prepend', 'prop_prefix',parameter_prefix, varargin{:});
0108 
0109   % set and capture numeric seed value
0110   if options.compile_flag==1
0111     % todo: make seed string (eg, 'shuffle') from param struct work with coder (options.compile_flag=1)
0112     % (currently raises error: "String input must be constant")
0113     % workaround: (shuffle here and get numeric seed for MEX-compatible params.mat)
0114     rng_wrapper(options.random_seed);
0115     options.random_seed=getfield(rng_wrapper,'Seed');  % <-- current active seed
0116     rng_function = 'rng';
0117   else
0118     rng_function = 'rng_wrapper';
0119   end
0120 
0121   % set parameter file name (save with m-file)
0122   [fpath,fname,fext]=fileparts2(options.filename);
0123   param_file_path = fullfile(fpath,'params.mat');
0124 
0125   % save parameters to disk
0126   warning('off','catstruct:DuplicatesFound');
0127   p = catstruct(dsCheckSolverOptions(options),model.parameters);
0128 
0129   if options.one_solve_file_flag
0130     % fill p flds that were varied with vectors of length = nSims
0131 
0132     vary=dsCheckOptions(varargin,{'vary',[],[],},false);
0133     vary = vary.vary;
0134 
0135     mod_set = dsVary2Modifications(vary);
0136     % The first 2 cols of modifications_set are idenitical to vary, it just
0137     % has the last column distributed out to the number of sims
0138 
0139 
0140     % Get param names
0141     iMod = 1;
0142     % Split extra entries in first 2 cols of mods, so each row is a single pop and param
0143     [~, first_mod_set] = dsApplyModifications([],mod_set{iMod}, varargin{:});
0144 
0145     % replace '->' with '_'
0146     first_mod_set(:,1) = strrep(first_mod_set(:,1), '->', '_');
0147 
0148     % add col of underscores
0149     first_mod_set = cat(2,first_mod_set(:,1), repmat({'_'},size(first_mod_set,1), 1), first_mod_set(:,2:end));
0150     nParamMods = size(first_mod_set, 1);
0151 
0152     % get param names
0153     mod_params = cell(nParamMods,1);
0154     for iRow = 1:nParamMods
0155       mod_params{iRow} = [first_mod_set{iRow,1:3}];
0156 
0157       %check if variable in namespace
0158       if ~any(strcmp(model.namespaces(:,2), mod_params{iRow}))
0159         % find correct entry based on param and pop
0160         nsInd = logical(~cellfun(@isempty, strfind(model.namespaces(:,2), [first_mod_set{iRow,1} '_'])) .* ...
0161           ~cellfun(@isempty, strfind(model.namespaces(:,2), first_mod_set{iRow,3})));
0162 
0163         assert(sum(nsInd) == 1)
0164 
0165         % add mech names using namespace
0166         mod_params{iRow} = model.namespaces{nsInd,2};
0167       end
0168     end
0169 
0170     % Get param values for each sim
0171     param_values = nan(nParamMods, length(mod_set));
0172     for iMod = 1:length(mod_set)
0173       % Split extra entries in first 2 cols of mods, so each row is a single pop and param
0174       [~, mod_set{iMod}] = dsApplyModifications([],mod_set{iMod}, varargin{:});
0175 
0176       % Get scalar values as vector
0177       param_values(:, iMod) = [mod_set{iMod}{:,3}];
0178     end
0179 
0180     % Assign value vectors to params
0181     for iParam = 1:nParamMods
0182       p.(mod_params{iParam}) = param_values(iParam,:);
0183     end
0184   end % one_solve_file_flag
0185 
0186   if options.verbose_flag
0187     fprintf('Saving params.mat\n');
0188   end
0189   save(param_file_path,'p','-v7');
0190 else
0191   % insert parameter values into model expressions
0192   model=dsPropagateParameters(model,'action','substitute', varargin{:});
0193 end
0194 
0195 % 1.2 prepare list of outputs (state variables and monitors)
0196 tmp=cellfun(@(x)[x ','],model.state_variables,'uni',0);
0197 tmp=[tmp{:}];
0198 output_string=tmp(1:end-1);
0199 
0200 if ~isempty(model.monitors)
0201   tmp=cellfun(@(x)[x ','],fieldnames(model.monitors),'uni',0);
0202   tmp=[tmp{:}];
0203   output_string=[output_string ',' tmp(1:end-1)];
0204 end
0205 
0206 if ~isempty(model.fixed_variables)
0207   tmp=cellfun(@(x)[x ','],fieldnames(model.fixed_variables),'uni',0);
0208   tmp=[tmp{:}];
0209   output_string=[output_string ',' tmp(1:end-1)];
0210 end
0211 
0212 output_string=['[T,' output_string ']']; % state vars, monitors, time vector
0213 
0214 % HACK to get IC to work
0215 if options.downsample_factor == 1
0216   for fieldNameC = fieldnames(model.ICs)'
0217     model.ICs.(fieldNameC{1}) = regexprep(model.ICs.(fieldNameC{1}), '_last', '(1,:)');
0218   end
0219 end
0220 
0221 %% 2.0 write m-file solver
0222 % 2.1 create mfile
0223 if ~isempty(options.filename)
0224   if options.verbose_flag
0225     fprintf('Creating solver file: %s\n',options.filename);
0226   end
0227 
0228   fid=fopen(options.filename,'wt');
0229 else
0230   fid=options.fileID;
0231 end
0232 
0233 outfile=fopen(fid);
0234 
0235 if options.disk_flag==1
0236   if ~options.one_solve_file_flag
0237     fprintf(fid,'function data_file=solve_ode\n');
0238   else
0239     fprintf(fid,'function data_file=solve_ode(simID)\n');
0240     if options.compile_flag
0241       fprintf(fid, 'assert(isa(simID, ''double''));\n');
0242     end
0243   end
0244 
0245   % create output data file
0246   fprintf(fid,'%% ------------------------------------------------------------\n');
0247   fprintf(fid,'%% Open output data file:\n');
0248   fprintf(fid,'data_file=''%s'';\n',options.data_file);
0249 
0250   %fprintf(fid,'fileID=fopen(data_file,''wt'');\n'); % <-- 'wt' does not
0251     % compile in linux. 't' may be necessary on PC. may need to look into this
0252   fprintf(fid,'fileID=fopen(data_file,''w'');\n');
0253 
0254   % write headers
0255   [state_var_counts,monitor_counts]=dsGetOutputCounts(model);
0256   fprintf(fid,'fprintf(fileID,''time%s'');\n',separator);
0257 
0258   if ~isempty(model.state_variables)
0259     for i=1:length(model.state_variables)
0260       fprintf(fid,'for i=1:%g, fprintf(fileID,''%s%s''); end\n',state_var_counts(i),model.state_variables{i},separator);
0261     end
0262   end
0263 
0264   if ~isempty(model.monitors)
0265     monitor_names=fieldnames(model.monitors);
0266     for i=1:length(monitor_names)
0267       fprintf(fid,'for i=1:%g, fprintf(fileID,''%s%s''); end\n',monitor_counts(i),monitor_names{i},separator);
0268     end
0269   end
0270 
0271   fprintf(fid,'fprintf(fileID,''\\n'');\n');
0272 else %options.disk_flag==0
0273   if ~options.one_solve_file_flag
0274     fprintf(fid,'function %s=solve_ode\n',output_string);
0275   else
0276     fprintf(fid,'function %s=solve_ode(simID)\n',output_string);
0277     if options.compile_flag
0278       fprintf(fid, 'assert(isa(simID, ''double''));\n');
0279     end
0280   end
0281 end
0282 
0283 % Benchmark tic
0284 if options.benchmark_flag
0285   fprintf(fid, 'tic;');
0286 end
0287 
0288 % 2.3 load parameters
0289 if options.save_parameters_flag
0290   fprintf(fid,'%% ------------------------------------------------------------\n');
0291   fprintf(fid,'%% Parameters:\n');
0292   fprintf(fid,'%% ------------------------------------------------------------\n');
0293   fprintf(fid,'params = load(''params.mat'',''p'');\n');
0294 
0295   if options.one_solve_file_flag && options.compile_flag
0296     fprintf(fid,'pVecs = params.p;\n');
0297   else
0298      fprintf(fid,'p = params.p;\n');
0299   end
0300 end
0301 
0302 if options.one_solve_file_flag
0303   % loop through p and for any vector, take simID index of it (ignores tspan)
0304   if ~options.compile_flag
0305     fprintf(fid,'\n%% For vector params, select index for this simID\n');
0306     fprintf(fid,'flds = fields(rmfield(p,''tspan''));\n'); % remove tspan
0307     fprintf(fid,'for fld = flds''\n');
0308     fprintf(fid,'  fld = fld{1};\n');
0309     fprintf(fid,'  if isnumeric(p.(fld)) && length(p.(fld)) > 1\n');
0310     fprintf(fid,'    p.(fld) = p.(fld)(simID);\n');
0311     fprintf(fid,'  end\n');
0312     fprintf(fid,'end\n\n');
0313   else %compile_flag
0314     % slice scalar from vector params
0315     for iParam = 1:nParamMods
0316       fprintf(fid,'p.%s = pVecs.%s(simID);\n', mod_params{iParam}, mod_params{iParam});
0317     end
0318 
0319     % take scalar from scalar params
0320     [~,sharedFlds] = intersect(fields(p), mod_params);
0321     scalar_params = fields(p);
0322     scalar_params(sharedFlds) = [];
0323     nScalarParams = length(scalar_params);
0324     for iParam = 1:nScalarParams
0325       fprintf(fid,'p.%s = pVecs.%s;\n', scalar_params{iParam}, scalar_params{iParam});
0326     end
0327   end
0328 end
0329 
0330 % write tspan, dt, and downsample_factor
0331 if options.save_parameters_flag
0332   fprintf(fid,'downsample_factor=%sdownsample_factor;\n',parameter_prefix);
0333   fprintf(fid,'dt=%sdt;\n',parameter_prefix);
0334   fprintf(fid,'T=(%stspan(1):dt:%stspan(2))'';\n',parameter_prefix,parameter_prefix);
0335 else
0336   fprintf(fid,'tspan=[%g %g];\ndt=%g;\ndownsample_factor=%g;\n',options.tspan,options.dt,options.downsample_factor);
0337   fprintf(fid,'T=(tspan(1):dt:tspan(2))'';\n');
0338 end
0339 
0340 % write calculation of time vector and derived parameters: ntime, nsamp
0341 fprintf(fid,'ntime=length(T);\nnsamp=length(1:downsample_factor:ntime);\n');
0342 
0343 % 2.4 evaluate fixed variables
0344 if ~isempty(model.fixed_variables)
0345   fprintf(fid,'%% ------------------------------------------------------------\n');
0346   fprintf(fid,'%% Fixed variables:\n');
0347   fprintf(fid,'%% ------------------------------------------------------------\n');
0348   % 2.2 set random seed
0349   fprintf(fid,'%% seed the random number generator\n');
0350   if options.save_parameters_flag
0351     fprintf(fid,'%s(%srandom_seed);\n',rng_function,parameter_prefix);
0352   else
0353     if ischar(options.random_seed)
0354       fprintf(fid,'%s(''%s'');\n',rng_function,options.random_seed);
0355     elseif isnumeric(options.random_seed)
0356       fprintf(fid,'%s(%g);\n',rng_function,options.random_seed);
0357     end
0358   end
0359   
0360   names=fieldnames(model.fixed_variables);
0361   expressions=struct2cell(model.fixed_variables);
0362   for i=1:length(names)
0363     fprintf(fid,'%s = %s;\n',names{i},expressions{i});
0364     if options.sparse_flag
0365       % create sparse matrix
0366       fprintf(fid,'%s = sparse(%s);\n',names{i},names{i});
0367     end
0368   end
0369 end
0370 
0371 % 2.5 evaluate function handles
0372 if ~isempty(model.functions) && options.reduce_function_calls_flag==0
0373   fprintf(fid,'%% ------------------------------------------------------------\n');
0374   fprintf(fid,'%% Functions:\n');
0375   fprintf(fid,'%% ------------------------------------------------------------\n');
0376   names=fieldnames(model.functions);
0377   expressions=struct2cell(model.functions);
0378   for i=1:length(names)
0379     fprintf(fid,'%s = %s;\n',names{i},expressions{i});
0380   end
0381 end
0382 
0383 % 2.6 prepare storage
0384 fprintf(fid,'%% ------------------------------------------------------------\n');
0385 fprintf(fid,'%% Initial conditions:\n');
0386 fprintf(fid,'%% ------------------------------------------------------------\n');
0387 
0388 % 2.2 set random seed (do this a 2nd time, so earlier functions don't mess
0389 % up the random seed)
0390 fprintf(fid,'%% seed the random number generator\n');
0391 if options.save_parameters_flag
0392   fprintf(fid,'%s(%srandom_seed);\n',rng_function,parameter_prefix);
0393 else
0394   if ischar(options.random_seed)
0395     fprintf(fid,'%s(''%s'');\n',rng_function,options.random_seed);
0396   elseif isnumeric(options.random_seed)
0397     fprintf(fid,'%s(%g);\n',rng_function,options.random_seed);
0398   end
0399 end
0400 
0401 % initialize time
0402 fprintf(fid,'t=0; k=1;\n');
0403 
0404 % todo: get coder varsize working with new format:
0405 
0406 % prepare for compilation
0407 % if options.compile_flag==1
0408 %   for i = 1:length(state_variables)
0409 %     %fprintf(fid,'coder.varsize(''%s'',[1e4,1],[true,false]);\n',state_variables{i}); % population size up to 1e4 for variable initial conditions
0410 %     fprintf(fid,'coder.varsize(''%s'',[1e8,1e4],[true,true]);\n',state_variables{i}); % population size up to 1e4 for variable initial conditions
0411 %   end
0412 % end
0413 
0414 % preallocate and initialize matrices to store results
0415 if options.disk_flag==1
0416   % add time zero to output data file
0417   fprintf(fid,'fprintf(fileID,''0%s'');\n',separator);
0418 end
0419 
0420 %% STATE_VARIABLES
0421 fprintf(fid,'\n%% STATE_VARIABLES:\n');
0422 nvals_per_var=zeros(1,length(state_variables)); % number of elements
0423 ndims_per_var=zeros(1,length(state_variables)); % number of dimensions
0424 IC_expressions=struct2cell(model.ICs);
0425 for i=1:length(state_variables)
0426   % initialize var_last
0427   if options.downsample_factor>1 || options.disk_flag==1
0428     % set var_last=IC;
0429     fprintf(fid,'%s_last = %s;\n',state_variables{i},IC_expressions{i});
0430   end
0431 
0432   if options.disk_flag==1
0433     % print var_last
0434     var_last=sprintf('%s_last',state_variables{i});
0435     fprintf(fid,'for i=1:numel(%s), fprintf(fileID,''%%g%s'',%s(i)); end\n',var_last,separator,var_last);
0436   else
0437     % preallocate state variables
0438     [pop_size,pop_name]=dsGetPopSizeFromName(model,state_variables{i});
0439     ndims_per_var(i)=length(pop_size);
0440     nvals_per_var(i)=prod(model.parameters.([pop_name '_Npop']));
0441     if options.save_parameters_flag
0442       % use pop size in saved params structure
0443       if ndims_per_var(i)==1
0444         % 1D population (time index is first dimension)
0445         fprintf(fid,'%s = zeros(nsamp,%s%s_Npop);\n',state_variables{i},parameter_prefix,pop_name);
0446       else
0447         % 2D population (time index is final dimension; will be shifted after simulation)
0448         % note: time index is last to avoid needing to squeeze() the matrix
0449         fprintf(fid,'%s = zeros([%s%s_Npop,nsamp]);\n',state_variables{i},parameter_prefix,pop_name);
0450       end
0451     else
0452       % hard-code the pop size
0453       if ndims_per_var(i)==1
0454         % 1D population (time index is first dimension)
0455         fprintf(fid,'%s = zeros(nsamp,%g);\n',state_variables{i},model.parameters.([pop_name '_Npop']));
0456       else
0457         % 2D population (time index is final dimension; will be shifted after simulation)
0458         % note: time index is last to avoid needing to squeeze() the matrix
0459         fprintf(fid,'%s = zeros([[%s],nsamp]);\n',state_variables{i},num2str(model.parameters.([pop_name '_Npop'])));
0460       end
0461     end
0462 
0463     % initialize state variables
0464     if options.downsample_factor==1
0465       if ndims_per_var(i)==1
0466         % set var(1,:)=IC;
0467         fprintf(fid,'  %s(1,:) = %s;\n',state_variables{i},IC_expressions{i});
0468       elseif ndims_per_var(i)==2
0469         % set var(:,:,1)=IC;
0470         fprintf(fid,'  %s(:,:,1) = %s;\n',state_variables{i},IC_expressions{i});
0471       else
0472         error('only 1D and 2D populations are supported a this time.');
0473       end
0474     else
0475       if ndims_per_var(i)==1
0476         % set var(1,:)=var_last;
0477         fprintf(fid,'  %s(1,:) = %s_last;\n',state_variables{i},state_variables{i});
0478       elseif ndims_per_var(i)==2
0479         % set var(:,:,1)=var_last;
0480         fprintf(fid,'  %s(:,:,1) = %s_last;\n',state_variables{i},state_variables{i});
0481       else
0482         error('only 1D and 2D populations are supported a this time.');
0483       end
0484     end
0485   end %disk_flag
0486 end %state_variables
0487 
0488 %% MONITORS
0489 if ~isempty(model.monitors)
0490   if strcmp(reportUI,'matlab') || options.disk_flag==1
0491     fprintf(fid,'\n%% MONITORS:\n');
0492   end
0493 
0494   monitor_names=fieldnames(model.monitors);
0495   monitor_expression=struct2cell(model.monitors);
0496 
0497   for i=1:length(monitor_names)
0498     if ~isempty(regexp(monitor_names{i},'_spikes$','once'))
0499     % set expression if monitoring spikes
0500       if options.disk_flag==1
0501         error('spike monitoring is not supported for writing data to disk at this time.');
0502         % todo: add support for spike monitoring with data written to
0503         % disk. approach: load data at end of simulation and do post-hoc
0504         % spike finding. if saving *.mat, use matfile() to load only the
0505         % variables in which to search. add spikes to output data file.
0506       end
0507 
0508       % number of spike times to store for each cell
0509       spike_buffer_size=2;%5;%100;
0510       % TODO: add support for: monitor VAR.spikes(thresh,buffer_size) (edit next if-then statements)
0511 
0512       if isempty(monitor_expression{i})
0513         spike_threshold=0;
0514       elseif isempty(regexp(monitor_expression{i},'[^\d]','once'))
0515         spike_threshold=str2num(monitor_expression{i});
0516         monitor_expression{i}=[];
0517       else
0518         spike_threshold=monitor_expression{i};
0519         monitor_expression{i}=[];
0520       end
0521 
0522       % approach: add conditional check for upward threshold crossing
0523       pop_name=regexp(monitor_names{i},'_','split');
0524       pop_name=pop_name{1};
0525       var_spikes=regexp(monitor_names{i},'(.*)_spikes$','tokens','once');
0526       var_spikes=var_spikes{1}; % variable to monitor
0527       var_tspikes=[pop_name '_tspike']; % only allow one event type to be tracked per population (i.e., it is ok to use pop_name, like 'E', as namespace instead of pop_var, like 'E_v')
0528       var_buffer_index=[pop_name '_buffer_index'];
0529 
0530       if ismember(var_spikes,model.state_variables)
0531         model.conditionals(end+1).namespace='spike_monitor';
0532 
0533         if (options.downsample_factor>1 || options.disk_flag==1)
0534           index_curr='_last';
0535         else
0536           index_curr='(n,:)';
0537         end
0538         index_last='(n-1,:)';
0539 
0540         if isnumeric(spike_threshold)
0541           model.conditionals(end).condition=...
0542             sprintf('%s%s>=%g&%s%s<%g',var_spikes,index_curr,spike_threshold,var_spikes,index_last,spike_threshold);
0543         else
0544           model.conditionals(end).condition=...
0545             sprintf('%s%s>=%s&%s%s<%s',var_spikes,index_curr,spike_threshold,var_spikes,index_last,spike_threshold);
0546         end
0547 
0548         action1=sprintf('%s(n,conditional_test)=1',monitor_names{i});
0549         action2=sprintf('inds=find(conditional_test); for j=1:length(inds), i=inds(j); %s(%s(i),i)=t; %s(i)=mod(-1+(%s(i)+1),%g)+1; end',var_tspikes,var_buffer_index,var_buffer_index,var_buffer_index,spike_buffer_size);
0550         model.conditionals(end).action=sprintf('%s;%s',action1,action2);
0551         model.conditionals(end).else=[];
0552         % move spike monitor to first position (ie.., to evaluate before other conditionals)
0553         model.conditionals=model.conditionals([length(model.conditionals) 1:length(model.conditionals)-1]);
0554         % remove from monitor list
0555         model.monitors=rmfield(model.monitors,monitor_names{i});
0556       end
0557 
0558       % initialize spike buffer and buffer index
0559       if options.save_parameters_flag
0560         % tspike = -inf(buffer_size,npop):
0561         fprintf(fid,'%s = -1e6*ones(%g,%s%s_Npop);\n',var_tspikes,spike_buffer_size,parameter_prefix,pop_name);
0562         fprintf(fid,'%s = ones(1,%s%s_Npop);\n',var_buffer_index,parameter_prefix,pop_name);
0563       else
0564         fprintf(fid,'%s = -1e6*ones(%g,%g);\n',var_tspikes,spike_buffer_size,model.parameters.([pop_name '_Npop']));
0565         fprintf(fid,'%s = ones(1,%g);\n',var_buffer_index,model.parameters.([pop_name '_Npop']));
0566       end
0567 
0568     elseif isempty(monitor_expression{i}) && isfield(model.functions,monitor_names{i})
0569     % set expression if monitoring function referenced by name
0570       tmp=regexp(model.functions.(monitor_names{i}),'@\([a-zA-Z][\w,]*\)\s*(.*)','tokens','once');
0571       monitor_expression{i}=tmp{1};
0572       model.monitors.(monitor_names{i})=tmp{1};
0573     end
0574 
0575     % initialize mon_last if not storing every time point and this is not a
0576     % spike monitor
0577     if (options.downsample_factor>1 || options.disk_flag==1) && isempty(regexp(monitor_names{i},'_spikes$','once'))
0578       % set mon_last=f(IC);
0579       tmp=cell2struct({monitor_expression{i}},{monitor_names{i}},1);
0580       print_monitor_update(fid,tmp,'_last',state_variables,'_last', varargin{:});
0581     end
0582 
0583     if options.disk_flag==1
0584       % print mon_last
0585       mon_last=sprintf('%s_last',monitor_names{i});
0586       fprintf(fid,'for i=1:numel(%s), fprintf(fileID,''%%g%s'',%s(i)); end\n',mon_last,separator,mon_last);
0587     elseif strcmp(reportUI,'matlab')
0588       % preallocate monitors
0589       [~,~,pop_name] = dsGetPopSizeFromName(model,monitor_names{i});
0590       if options.save_parameters_flag
0591         fprintf(fid,'%s = zeros(nsamp,%s%s_Npop);\n',monitor_names{i},parameter_prefix,pop_name);
0592       else
0593         fprintf(fid,'%s = zeros(nsamp,%g);\n',monitor_names{i},model.parameters.([pop_name '_Npop']));
0594       end
0595 %       parts=regexp(monitor_names{i},'_','split');
0596 %       if options.save_parameters_flag
0597 %         fprintf(fid,'%s = zeros(nsamp,%s%s_Npop);\n',monitor_names{i},parameter_prefix,parts{1});
0598 %       else
0599 %         fprintf(fid,'%s = zeros(nsamp,%g);\n',monitor_names{i},model.parameters.([parts{1} '_Npop']));
0600 %       end
0601 
0602       if isempty(monitor_expression{i})
0603         continue;
0604       end
0605 
0606       % initialize monitors
0607       if options.downsample_factor==1
0608         % set mon(1,:)=f(IC);
0609         tmp=cell2struct({monitor_expression{i}},{monitor_names{i}},1);
0610         print_monitor_update(fid,tmp,'(1,:)',state_variables,'(1,:)', varargin{:});
0611       else
0612         % set mon(1,:)=mon_last;
0613         tmp=cell2struct({monitor_expression{i}},{monitor_names{i}},1);
0614         print_monitor_update(fid,tmp,'(1,:)',state_variables,'_last', varargin{:});
0615       end
0616     end %disk_flag
0617   end %monitor_names
0618 end %monitors
0619 
0620 if options.disk_flag==1
0621   % go to new line for next time point
0622   fprintf(fid,'fprintf(fileID,''\\n'');\n');
0623 end
0624 
0625 % determine how to index each state variable based on how often state
0626 % variables are stored, whether they are written to disk or stored in
0627 % memory, and whether the state variable matrix is one- or two-dimensional
0628 index_lasts=cell(1,length(state_variables));
0629 index_nexts=cell(1,length(state_variables));
0630 index_temps=repmat({'_last'},[1 length(state_variables)]);
0631 for i=1:length(state_variables)
0632   if options.downsample_factor==1 && options.disk_flag==0
0633     % store state directly into state variables on each integration step
0634     if nvals_per_var(i)>1 % use full 2D matrix indexing
0635       if ndims_per_var(i)==1 % 1D population
0636         index_lasts{i}='(n-1,:)';
0637         index_nexts{i}='(n,:)';
0638       elseif ndims_per_var(i)==2 % 2D population
0639         index_lasts{i}='(:,:,n-1)';
0640         index_nexts{i}='(:,:,n)';
0641       end
0642     else % use more concise 1D indexing because it is much faster for some Matlab-specific reason...
0643       index_lasts{i}='(n-1)';
0644       index_nexts{i}='(n)';
0645     end
0646   elseif options.downsample_factor>1 && options.disk_flag==0
0647     % store state in var_last then update state variables on each downsample_factor integration step
0648     index_lasts{i}='_last';
0649     if nvals_per_var(i)>1
0650       if ndims_per_var(i)==1 % 1D population
0651         index_nexts{i}='(n,:)';
0652       elseif ndims_per_var(i)==2 % 2D population
0653         index_nexts{i}='(:,:,n)';
0654       end
0655     else
0656       index_nexts{i}='(n)';
0657     end
0658   elseif options.disk_flag==1
0659     % always store state in var_last and write on each downsample_factor integration step
0660       index_lasts{i}='_last';
0661       index_nexts{i}='_last';
0662   end
0663 end
0664 
0665 % add index to state variables in ODEs and look for delay differential equations
0666 delayinfo=[];
0667 odes = struct2cell(model.ODEs);
0668 for i=1:length(odes)
0669   for j=1:length(state_variables)
0670     odes{i}=dsStrrep(odes{i}, state_variables{j}, [state_variables{j} index_lasts{j}], '', '', varargin{:});
0671     % #####################################################################
0672     % COLLECT DELAY DIFFERENTIAL EQUATION INFO
0673     % search for delays: [state_variables{j} index_lasts{j} '(t-']
0674     % note: this only works for 1D populations
0675     tmp=strrep(strrep([state_variables{j} index_lasts{j}],'(','\('),')','\)');
0676     if ~isempty(regexp(odes{i},[tmp '\(t'],'once'))
0677       matches=regexp(odes{i},[tmp '\(t-[\w\.,:]+\)'],'match');
0678       for k=1:length(matches)
0679         % determine amount of delay for each occurrence of a delay to this state variable
0680         %     note: account for user-specified X(t-tau) and X(t-tau,:)
0681         %     note: account for tau as variable defined elsewhere or numeric
0682         % look for: X(t-#)
0683         delay=cellstr2num(regexp(matches{k},'\(t-([\.\d]+)\)','tokens','once'));
0684         if isempty(delay)
0685           % look for: X(t-#,:)
0686           delay=cellstr2num(regexp(matches{k},'\(t-([\.\d]+),:\)','tokens','once'));
0687         end
0688         if isempty(delay)
0689           % look for: X(t-param)
0690           delay=regexp(matches{k},'\(t-([\w\.]+)\)','tokens','once');
0691         end
0692         if isempty(delay)
0693           % look for: X(t-param,:)
0694           delay=regexp(matches{k},'\(t-([\w\.]+),:\)','tokens','once');
0695         end
0696         if iscell(delay) && ischar(delay{1})
0697           delay=strrep(delay{1},parameter_prefix,''); % remove parameter prefix
0698           delay=strrep(delay,',:',''); % remove population dimension from index to delay matrix
0699           % look for parameter with delay length
0700           if isfield(model.parameters,delay)
0701             delay=model.parameters.(delay);
0702           else
0703             error('delay parameter ''%s'' not found.',delay);
0704           end
0705         end
0706         if ~isempty(delay) && isnumeric(delay)
0707           delay_samp = ceil(delay/options.dt);
0708           delayinfo(end+1).variable=state_variables{j};
0709           delayinfo(end).strmatch=matches{k};
0710           delayinfo(end).delay_samp=delay_samp;
0711           delayinfo(end).ode_index=i;
0712         end
0713       end
0714     end
0715     % #####################################################################
0716   end
0717 end
0718 % #####################################################################
0719 % PROCESS DELAY DIFFERENTIAL EQUATION INFO
0720 % determine max delay for each state variable
0721 if ~isempty(delayinfo)
0722   delay_vars=unique({delayinfo.variable});
0723   delay_maxi=zeros(size(delay_vars));
0724   for i=1:length(delay_vars)
0725     if i==1
0726       fprintf(fid,'%% DELAY MATRICES:\n');
0727     end
0728     % find max delay for this state variable
0729     idx=ismember({delayinfo.variable},delay_vars{i});
0730     Dmax=max([delayinfo(idx).delay_samp]);
0731     delay_maxi(i)=Dmax;
0732     % convert delay indices into delay vector indices based on max delay
0733     tmps=num2cell(Dmax-[delayinfo(idx).delay_samp]);
0734     [delayinfo(idx).delay_index]=deal(tmps{:});
0735     % initialize delay matrix with max delay ICs and all time points
0736     fprintf(fid,'%s_delay = zeros(nsamp+%g,size(%s,2));\n',delayinfo(i).variable,Dmax,delayinfo(i).variable);
0737     fprintf(fid,'  %s_delay(1:%g,:) = repmat(%s(1,:),[%g 1]);\n',delayinfo(i).variable,Dmax,delayinfo(i).variable,Dmax);
0738   end
0739   % replace delays in ODEs with indices to delay matrices
0740   for i=1:length(delayinfo)
0741     rep=sprintf('%s_delay(k+%g,:)',delayinfo(i).variable,delayinfo(i).delay_index);
0742     odes{delayinfo(i).ode_index}=strrep(odes{delayinfo(i).ode_index},delayinfo(i).strmatch,rep);
0743   end
0744 end
0745 % #####################################################################
0746 % remove unused @linkers from ODEs
0747 for i=1:length(odes)
0748   if any(odes{i}=='@')
0749     tmp=regexp(odes{i},'@([\w_]+)','tokens');
0750     if ~isempty(tmp)
0751       tmp=[tmp{:}];
0752       for j=1:length(tmp)
0753         odes{i}=strrep(odes{i},['@' tmp{j}],'0');
0754       end
0755     end
0756   end
0757 end
0758 % #####################################################################
0759 
0760 %% Numerical integration
0761 % write code to do numerical integration
0762 fprintf(fid,'%% ###########################################################\n');
0763 fprintf(fid,'%% Numerical integration:\n');
0764 fprintf(fid,'%% ###########################################################\n');
0765 fprintf(fid,'n=2;\n');
0766 fprintf(fid,'for k=2:ntime\n'); % time index
0767 fprintf(fid,'  t=T(k-1);\n');
0768 if options.downsample_factor==1 && options.disk_flag==0 % store every time point, do not use var_last
0769   % update_vars;      % var(:,k-1)->var(k,:) or var(k-1)->var(k)
0770   update_vars(index_nexts, varargin{:});
0771 
0772   % conditionals;     % var(k,:)->var(k,:) or var(k)->var(k)
0773   print_conditional_update(fid,model.conditionals,index_nexts,state_variables)
0774 
0775   if strcmp(reportUI,'matlab') % update_monitors;  % mon(:,k-1)->mon(k,:) or mon(k-1)->mon(k)
0776     print_monitor_update(fid,model.monitors,index_nexts,state_variables, [], varargin{:});
0777   end
0778 
0779   fprintf(fid,'  n=n+1;\n');
0780 else % store every downsample_factor time point in memory or on disk
0781   % update_vars;      % var_last->var_last
0782   update_vars(index_temps, varargin{:});
0783 
0784   % conditionals;     % var_last->var_last
0785   print_conditional_update(fid,model.conditionals,index_temps,state_variables, varargin{:});
0786 
0787   % check if it is time to store data
0788   fprintf(fid,'if mod(k,downsample_factor)==0 %% store this time point\n');
0789 
0790   if options.disk_flag==1 % write to disk
0791     fprintf(fid,'  %% ------------------------------------------------------------\n');
0792     fprintf(fid,'  %% Write state variables and monitors to disk:\n');
0793     fprintf(fid,'  %% ------------------------------------------------------------\n');
0794 
0795     % print current time point to data file
0796     fprintf(fid,'  fprintf(fileID,''%%g%s'',T(k));\n',separator);
0797 
0798     % print state variables
0799     for i=1:length(state_variables)
0800       var_last=sprintf('%s_last',state_variables{i});
0801       fprintf(fid,'  for i=1:numel(%s), fprintf(fileID,''%%g%s'',%s(i)); end\n',var_last,separator,var_last);
0802     end
0803 
0804     % print monitors
0805     if ~isempty(model.monitors)
0806       for i=1:length(monitor_names)
0807           mon_last=sprintf('%s_last',monitor_names{i});
0808           fprintf(fid,'  for i=1:numel(%s), fprintf(fileID,''%%g%s'',%s(i)); end\n',mon_last,separator,mon_last);
0809       end
0810     end
0811 
0812     fprintf(fid,'  fprintf(fileID,''\\n'');\n');
0813   else            % store in memory
0814     % update_vars;    % var_last -> var(n,:) or var(n)
0815     print_var_update_last(fid,index_nexts,state_variables)
0816     if strcmp(reportUI,'matlab') % update_monitors;% f(var_last) -> mon(n,:) or mon(n)
0817       print_monitor_update(fid,model.monitors,index_nexts,state_variables, [], varargin{:});
0818     end
0819   end %disk_flag
0820 
0821   fprintf(fid,'  n=n+1;\n');
0822   fprintf(fid,'end\n');
0823 end
0824 % update delay matrices
0825 if ~isempty(delayinfo)
0826   fprintf(fid,'  %% ------------------------------------------------------------\n');
0827   fprintf(fid,'  %% Update delay matrices:\n');
0828   fprintf(fid,'  %% ------------------------------------------------------------\n');
0829   for i=1:length(delay_vars)
0830     if options.downsample_factor==1 && options.disk_flag==0
0831       fprintf(fid,'  %s_delay(k+%g,:)=%s(k,:);\n',delay_vars{i},delay_maxi(i),delay_vars{i});
0832     else
0833       fprintf(fid,'  %s_delay(k+%g,:)=%s_last;\n',delay_vars{i},delay_maxi(i),delay_vars{i});
0834     end
0835   end
0836 end
0837 
0838 fprintf(fid,'end\n');
0839 if ~isempty(model.monitors) && ~strcmp(reportUI,'matlab') && options.disk_flag==0 % computing monitors outside the solver loop
0840   fprintf(fid,'%% ------------------------------------------------------------\n');
0841   fprintf(fid,'%% Compute monitors:\n');
0842   fprintf(fid,'%% ------------------------------------------------------------\n');
0843 
0844   monitor_name=fieldnames(model.monitors);
0845   monitor_expression=struct2cell(model.monitors);
0846 
0847   % add indexes to state variables in monitors
0848   for i=1:length(monitor_name)
0849     % hacks to vectorize expression in Octave:
0850     monitor_vectorized_expression = dsStrrep(monitor_expression{i},'t','T');
0851     monitor_vectorized_expression = strrep(monitor_vectorized_expression,'1,p.pop','length(T),p.pop');
0852     monitor_vectorized_expression = strrep(monitor_vectorized_expression,'(k,:)','');
0853     % write monitors to solver function
0854     fprintf(fid,'%s=%s;\n',monitor_name{i},monitor_vectorized_expression);
0855   end
0856 end
0857 fprintf(fid,'\nT=T(1:downsample_factor:ntime);\n');
0858 
0859 if any(ndims_per_var>1) % is there at least one 2D population?
0860   for i=1:length(state_variables)
0861     if ndims_per_var(i)>1
0862       % move time index to first dimension
0863       fprintf(fid,'%s=shiftdim(%s,%g);\n',state_variables{i},state_variables{i},ndims_per_var(i));
0864     end
0865   end
0866 end
0867 
0868 % cleanup
0869 if options.disk_flag==1
0870   fprintf(fid,'%% ------------------------------------------------------------\n');
0871   fprintf(fid,'%% Close output data file:\n');
0872   fprintf(fid,'fclose(fileID);\n');
0873   fprintf(fid,'%% ------------------------------------------------------------\n');
0874 end
0875 
0876 %% Benchmark toc
0877 if options.benchmark_flag
0878   fprintf(fid, 'fprintf(''Sim Time: %%g seconds\\n'', toc);');
0879 end
0880 
0881 %% end solve function
0882 fprintf(fid,'\nend\n\n');
0883 
0884 if ~strcmp(outfile,'"stdout"')
0885   fclose(fid);
0886   % wait for file before continuing to simulation
0887   while ~exist(outfile,'file')
0888     pause(.01);
0889   end
0890 end
0891 
0892   % ########################################
0893   % NESTED FUNCTIONS
0894   % ########################################
0895   function update_vars(index_nexts_, varargin)
0896     switch options.solver
0897       case {'euler','rk1'}
0898         print_k(fid,odes,'_k1',state_variables);                              % write k1 using model.ODEs
0899 
0900         % write update for state variables
0901         print_var_update(fid,index_nexts_,index_lasts,...
0902           'dt*%s_k1',state_variables);
0903       case {'rk2','modified_euler'}
0904         print_k(fid,odes,'_k1',state_variables);                              % write k1 using model.ODEs
0905 
0906         odes_k2=update_odes(odes,'_k1','.5*dt',state_variables,index_lasts, varargin{:});   % F(*,yn+dt*k1/2)
0907         fprintf(fid,'  t=t+.5*dt;\n');                                          % update t before writing k2
0908         print_k(fid,odes_k2,'_k2',state_variables);                           % write k2 using odes_k2
0909 
0910         % write update for state variables
0911         print_var_update(fid,index_nexts_,index_lasts,...
0912           'dt*%s_k2',state_variables);
0913       case {'rk4','rungekutta','rk'}
0914         print_k(fid,odes,'_k1',state_variables);                              % write k1 using model.ODEs
0915 
0916         odes_k2=update_odes(odes,'_k1','.5*dt',state_variables,index_lasts, varargin{:});   % F(*,yn+dt*k1/2)
0917         fprintf(fid,'  t=t+.5*dt;\n');                                          % update t before writing k2
0918         print_k(fid,odes_k2,'_k2',state_variables);                           % write k2 using odes_k2
0919 
0920         odes_k3=update_odes(odes,'_k2','.5*dt',state_variables,index_lasts, varargin{:});   % F(*,yn+dt*k2/2)
0921         print_k(fid,odes_k3,'_k3',state_variables);                           % write k3 using odes_k3
0922 
0923         odes_k4=update_odes(odes,'_k3','dt',state_variables,index_lasts, varargin{:});      % F(*,yn+dt*k3)
0924         fprintf(fid,'  t=t+.5*dt;\n');                                          % update t before writing k4
0925         print_k(fid,odes_k4,'_k4',state_variables);                           % write k4 using odes_k4
0926 
0927         % write update for state variables
0928         print_var_update(fid,index_nexts_,index_lasts,...
0929           '(dt/6)*(%s_k1+2*(%s_k2+%s_k3)+%s_k4)',state_variables);
0930     end
0931   end
0932 
0933 end %main
0934 
0935 
0936 % ########################################
0937 %% SUBFUNCTIONS
0938 % ########################################
0939 
0940 function print_k(fid,odes_k,suffix_k,state_variables,nvals_per_var)
0941   % purpose: write auxiliary calculations (k1-k4) for runge-kutta
0942   for i=1:length(odes_k)
0943     fprintf(fid,'  %s%s=%s;\n',state_variables{i},suffix_k,odes_k{i});
0944   end
0945 end
0946 
0947 function odes_out=update_odes(odes,suffix_k,increment,state_variables,index_lasts, varargin)
0948   % purpose: update expressions for axiliary calculations (k1-k4)
0949   odes_out=odes;
0950   for i=1:length(odes)
0951     for j=1:length(odes)
0952       tmp=[state_variables{j} index_lasts{j}]; % original state variable as appears in ODEs
0953       tmp=strrep(strrep(tmp,')','\)'),'(','\('); % escape parentheses for substitution
0954       odes_out{i}=dsStrrep(odes_out{i}, tmp, sprintf('(%s%s+%s*%s%s)', state_variables{j}, index_lasts{j}, increment, state_variables{j}, suffix_k), '(',')', varargin{:});
0955     end
0956   end
0957 end
0958 
0959 function print_var_update(fid,index_nexts,index_lasts,update_term,state_variables)
0960   % purpose: write statements to update state variables according to their dynamics
0961   % example:
0962   %   update_term='(dt/6)*(%s_k1+2*(%s_k2+%s_k3)+%s_k4)';
0963   %   state_variable='A_v';
0964 
0965   if ~isempty(state_variables)
0966     fprintf(fid,'  %% ------------------------------------------------------------\n');
0967     fprintf(fid,'  %% Update state variables:\n');
0968     fprintf(fid,'  %% ------------------------------------------------------------\n');
0969   else
0970     return;
0971   end
0972 
0973   for i=1:length(state_variables)
0974     this_update_term=strrep(update_term,'%s',state_variables{i});
0975     this_update_expression=sprintf('%s%s=%s%s+%s;',state_variables{i},index_nexts{i},state_variables{i},index_lasts{i},this_update_term);
0976     fprintf(fid,'  %s\n',this_update_expression);
0977   end
0978 end
0979 
0980 function print_var_update_last(fid,index_nexts,state_variables)
0981   if ~isempty(state_variables)
0982     fprintf(fid,'  %% ------------------------------------------------------------\n');
0983     fprintf(fid,'  %% Store state variables:\n');
0984     fprintf(fid,'  %% ------------------------------------------------------------\n');
0985   else
0986     return;
0987   end
0988 
0989   for i=1:length(state_variables)
0990     this_update_expression=sprintf('%s%s=%s_last;\n',state_variables{i},index_nexts{i},state_variables{i});
0991     fprintf(fid,'  %s',this_update_expression);
0992   end
0993 end
0994 
0995 function print_conditional_update(fid,conditionals,index_nexts,state_variables, varargin)
0996   % purpose: write statements to perform conditional actions (that may
0997   %   include updating state variables).
0998   if ~isempty(conditionals)
0999     fprintf(fid,'  %% ------------------------------------------------------------\n');
1000     fprintf(fid,'  %% Conditional actions:\n');
1001     fprintf(fid,'  %% ------------------------------------------------------------\n');
1002   else
1003     return;
1004   end
1005 
1006   switch index_nexts{1}(1)
1007     case '('
1008       action_index='(n,conditional_test)';
1009     case '_'
1010       action_index='_last(conditional_test)';
1011   end
1012 
1013   for i=1:length(conditionals)
1014     condition=conditionals(i).condition;
1015     action=conditionals(i).action;
1016     elseaction=conditionals(i).else;
1017 
1018     % add indexes to state variables in conditional actions
1019     for j=1:length(state_variables)
1020       if strcmp('spike_monitor',conditionals(i).namespace)
1021         % do nothing if spike_monitor
1022       else
1023         condition=dsStrrep(condition, state_variables{j}, [state_variables{j} index_nexts{j}], '', '', varargin{:});
1024       end
1025 
1026       action=dsStrrep(action, state_variables{j}, [state_variables{j} action_index], '', '', varargin{:});
1027 
1028       if ~isempty(elseaction)
1029         elseaction=dsStrrep(elseaction, state_variables{j}, [state_variables{j} action_index], '', '', varargin{:});
1030       end
1031     end
1032 
1033     % write conditional to solver function
1034     fprintf(fid,'  conditional_test=(%s);\n',condition);
1035     action=dsStrrep(action, '\(n,:', '(n,conditional_test', '', '', varargin{:});
1036     fprintf(fid,'  if any(conditional_test), %s; ',action);
1037 
1038     if ~isempty(elseaction)
1039       elseaction=dsStrrep(elseaction, '(n,:', '(n,conditional_test', '', '', varargin{:});
1040       fprintf('else %s; ',elseaction);
1041     end
1042 
1043     fprintf(fid,'end\n');
1044   end
1045 end
1046 
1047 function print_monitor_update(fid,monitors,index_nexts,state_variables,index_lasts, varargin)
1048   if ~isempty(monitors) && iscell(index_nexts) % being called from within the integrator loop
1049     fprintf(fid,'  %% ------------------------------------------------------------\n');
1050     fprintf(fid,'  %% Update monitors:\n');
1051     fprintf(fid,'  %% ------------------------------------------------------------\n');
1052   end
1053 
1054   if nargin<5 || isempty(index_lasts)
1055     index_lasts=index_nexts;
1056   end
1057 
1058   % account for inputs from monitor initialization
1059 
1060   if ~iscell(index_nexts), index_nexts={index_nexts}; end
1061 
1062   if ~iscell(index_lasts), index_lasts={index_lasts}; end
1063 
1064   if length(index_lasts)~=length(state_variables)
1065     index_lasts=repmat(index_lasts,[1 length(state_variables)]);
1066   end
1067 
1068   if isequal(index_nexts{1},'(1,:)')
1069     monitor_index='(1,:)';
1070   else
1071     switch index_nexts{1}(1)
1072       case '('
1073         monitor_index='(n,:)';
1074       case '_'
1075         monitor_index='_last';
1076     end
1077   end
1078 
1079   % purpose: write statements to update monitors given current state variables
1080   %   note: only run this every time a point is recorded (not necessarily on
1081   %   every step of the integration).
1082   if isempty(monitors)
1083     return;
1084   end
1085 
1086   monitor_name=fieldnames(monitors);
1087   monitor_expression=struct2cell(monitors);
1088 
1089   % add indexes to state variables in monitors
1090   for i=1:length(monitor_name)
1091     for j=1:length(state_variables)
1092       %monitor_expression{i}=dsStrrep(monitor_expression{i}, state_variables{j}, [state_variables{j} index_lasts{j}], '', '', varargin{:});
1093       monitor_expression{i}=dsStrrep(monitor_expression{i}, state_variables{j}, [state_variables{j} monitor_index], '', '', varargin{:});
1094     end
1095 
1096     % write monitors to solver function
1097     fprintf(fid,'  %s%s=%s;\n',monitor_name{i},monitor_index,monitor_expression{i});
1098   end
1099 end
1100 
1101 %{
1102 PSEUDOCODE:
1103 % --------------------------------------------------------------------------------
1104 % setup (parameters, fixed_variables, functions)
1105 % preallocation and initial conditions
1106 if downsample_factor>1 || disk_flag==1
1107   % var_last=IC;
1108   % mon_last=f(IC);
1109 end
1110 if disk_flag==1
1111   fid=fopen(options.filename,'wt');
1112   % print var_last
1113   % print mon_last
1114 else
1115   % var=zeros(npop,nsamp);
1116   % mon=zeros(npop,nsamp);
1117   if downsample_factor==1
1118     % var(1,:)=IC;
1119     % mon(1,:)=f(IC);
1120   else
1121     % var(1,:)=var_last;
1122     % mon(1,:)=mon_last;
1123   end
1124 end
1125 
1126 % numerical integration
1127 % n=1; % storage index
1128 % for k=2:ntime % time index
1129   if downsample_factor==1 && disk_flag==0 % do not use var_last
1130     % update_vars;      % var(:,k-1)->var(k,:) or var(k-1)->var(k)
1131     % conditionals;     % var(k,:)->var(k,:) or var(k)->var(k)
1132     % update_monitors;  % mon(:,k-1)->mon(k,:) or mon(k-1)->mon(k)
1133   else % downsample_factor>1 and/or disk_flag==1
1134     % update_vars;      % var_last->var_last
1135     % conditionals;     % var_last->var_last
1136     if mod(t,downsample_factor)==0 % store this time point
1137       if disk_flag==1 % write to disk
1138         % write_vars;     % print: var_last
1139         % write_monitors; % print: mon=f(var_last)
1140       else            % store in memory
1141         % update_vars;    % var_last -> var(n,:) or var(n)
1142         % update_monitors;% f(var_last) -> mon(n,:) or mon(n)
1143       end
1144       % n=n+1;
1145     end
1146   end
1147 % end
1148 % --------------------------------------------------------------------------------
1149 %}