dsWriteMatlabSolver

WRITEMATLABSOLVER - write m-file that numerically integrates the model

 Usage:
   filepath = dsWriteMatlabSolver(model,varargin)

 Inputs:
   - model: DynaSim model structure (see dsGenerateModel)
   - options:
     'tspan'         : units must be consistent with dt and equations
                       {[beg,end]} (default: [0 100])
     'ic'            : initial conditions; this overrides definition in model structure
     'solver'        : built-in Matlab solvers
                         {'ode23','ode45','ode113','ode15s','ode23s','ode23t','ode23tb'}
     'matlab_solver_options': options from odeset for use with built-in Matlab solvers
     'dt'            :  time step used for fixed step DSSim solvers (default: 0.01)
     'modifications' : DynaSim modifications structure
     'reduce_function_calls_flag': whether to eliminate internal function
                                   calls {0 or 1} (default: 1)
     'coder_flag'    : whether to compile using coder instead of interpreting
                       Matlab (default: exist('codegen')==6 TODO is this correct?
                       what does this mean?)
     'downsample_factor': downsampling applied during simulation. Only every
                          downsample_factor-time point is stored in memory or
                          written to disk (default: 1)
     'random_seed'   : seed for random number generator (usage:
                       rng(random_seed)) (default: now)

 Outputs:
   - filepath (solve_ode.m)
   - odefun_filepath (solve_ode_odefun.m)

 Dependencies: dsCheckOptions, dsCheckModel

 See also: dsSimulate, dsDynasim2odefun

0001 function solve_ode_filepath = dsWriteMatlabSolver(model,varargin)
0002 %WRITEMATLABSOLVER - write m-file that numerically integrates the model
0003 %
0004 % Usage:
0005 %   filepath = dsWriteMatlabSolver(model,varargin)
0006 %
0007 % Inputs:
0008 %   - model: DynaSim model structure (see dsGenerateModel)
0009 %   - options:
0010 %     'tspan'         : units must be consistent with dt and equations
0011 %                       {[beg,end]} (default: [0 100])
0012 %     'ic'            : initial conditions; this overrides definition in model structure
0013 %     'solver'        : built-in Matlab solvers
0014 %                         {'ode23','ode45','ode113','ode15s','ode23s','ode23t','ode23tb'}
0015 %     'matlab_solver_options': options from odeset for use with built-in Matlab solvers
0016 %     'dt'            :  time step used for fixed step DSSim solvers (default: 0.01)
0017 %     'modifications' : DynaSim modifications structure
0018 %     'reduce_function_calls_flag': whether to eliminate internal function
0019 %                                   calls {0 or 1} (default: 1)
0020 %     'coder_flag'    : whether to compile using coder instead of interpreting
0021 %                       Matlab (default: exist('codegen')==6 TODO is this correct?
0022 %                       what does this mean?)
0023 %     'downsample_factor': downsampling applied during simulation. Only every
0024 %                          downsample_factor-time point is stored in memory or
0025 %                          written to disk (default: 1)
0026 %     'random_seed'   : seed for random number generator (usage:
0027 %                       rng(random_seed)) (default: now)
0028 %
0029 % Outputs:
0030 %   - filepath (solve_ode.m)
0031 %   - odefun_filepath (solve_ode_odefun.m)
0032 %
0033 % Dependencies: dsCheckOptions, dsCheckModel
0034 %
0035 % See also: dsSimulate, dsDynasim2odefun
0036 
0037 % Check inputs
0038 options=dsCheckOptions(varargin,{...
0039   'ic',[],[],...                  % initial conditions (overrides definition in model structure)
0040   'tspan',[0 100],[],...          % [beg,end] (units must be consistent with dt and equations)
0041   'dt',.01,[],...                 % time step used for fixed step DynaSim solvers
0042   'downsample_factor',1,[],...    % downsampling applied after simulation (only every downsample_factor-time point is returned)
0043   'random_seed','shuffle',[],...        % seed for random number generator (usage: rng(random_seed))
0044   'solver','ode45',{'ode23','ode45','ode113','ode15s','ode23s','ode23t','ode23tb'},... % built-in Matlab solvers
0045   'solver_type','matlab',{'matlab', 'matlab_no_mex'},... % if compile_flag==1, will decide whether to mex solve_file or odefun_file
0046   'matlab_solver_options',[],[],... % options from odeset for use with built-in Matlab solvers
0047   'reduce_function_calls_flag',1,{0,1},...   % whether to eliminate internal (anonymous) function calls
0048   'save_parameters_flag',1,{0,1},...
0049   'filename',[],[],...         % name of solver file that integrates model
0050   'fileID',1,[],...
0051   'compile_flag',0,{0,1},... % whether to prepare script for being compiled using coder instead of interpreting Matlab
0052   'verbose_flag',1,{0,1},...
0053   'one_solve_file_flag',0,{0,1},... % use only 1 solve file of each type, but can't vary mechs yet
0054   'benchmark_flag',0,{0,1},...
0055   },false);
0056 
0057 % Check inputs
0058 model=dsCheckModel(model, varargin{:});
0059 
0060 % convert matlab solver options from key/value to struct using odeset if necessary
0061 if iscell(options.matlab_solver_options) && ~isempty(options.matlab_solver_options)
0062   options.matlab_solver_options = odeset(options.matlab_solver_options{:});
0063 end
0064 
0065 %% 1.0 Get ode_fun
0066 
0067 % create function that calls feval(@solver,...) and has subfunction
0068 % defining odefun (including optional conditionals)...
0069 
0070 propagatedModel = dsPropagateParameters( dsPropagateFunctions(model, varargin{:}), varargin{:} );
0071 propagatedModel = dsPropagateParameters(propagatedModel, 'param_type', 'fixed_variables', varargin{:});
0072 [odefun,IC,elem_names] = dsDynasim2odefun(propagatedModel, 'odefun_output','func_body', varargin{:});
0073 
0074 
0075 %% 2.0 prepare model info
0076 parameter_prefix='p.';%'pset.p.';
0077 % state_variables=model.state_variables;
0078 
0079 % 1.1 eliminate internal (anonymous) function calls from model equations
0080 % if options.reduce_function_calls_flag==1
0081   model=dsPropagateFunctions(model, varargin{:});
0082 % end
0083 
0084 % 1.1 prepare parameters
0085 if options.save_parameters_flag
0086   % add parameter struct prefix to parameters in model equations
0087   model=dsPropagateParameters(model,'action','prepend','prop_prefix',parameter_prefix, varargin{:});
0088 
0089   % set and capture numeric seed value
0090   if options.compile_flag==1
0091     % todo: make seed string (eg, 'shuffle') from param struct work with coder (options.compile_flag=1)
0092     % (currently raises error: "String input must be constant")
0093     % workaround: (shuffle here and get numeric seed for MEX-compatible params.mat)
0094     rng_wrapper(options.random_seed);
0095     options.random_seed=getfield(rng_wrapper,'Seed');  % <-- current active seed
0096     rng_function = 'rng';
0097   else
0098     rng_function = 'rng_wrapper';
0099   end
0100 
0101   % set parameter file name (save with m-file)
0102   [fpath,fname,fext]=fileparts2(options.filename);
0103   odefun_filename = [fname '_odefun'];
0104   param_file_name = fullfile(fpath,'params.mat');
0105 
0106   % save parameters to disk
0107   warning('off','catstruct:DuplicatesFound');
0108 
0109   % make p struct
0110   p=catstruct(dsCheckSolverOptions(options),model.parameters);
0111 
0112   % add IC to p
0113   %   NOTE: will get done again in simulateModel
0114   if isempty(options.ic)
0115     p.ic = IC;
0116   else %if overridden from options
0117     p.ic = options.ic;
0118   end
0119 
0120   % add matlab_solver_options to p
0121   if ~isempty(options.matlab_solver_options)
0122     p.matlab_solver_options = options.matlab_solver_options;
0123   end
0124 
0125   if options.one_solve_file_flag
0126     % fill p flds that were varied with vectors of length = nSims
0127 
0128     vary=dsCheckOptions(varargin,{'vary',[],[],},false);
0129     vary = vary.vary;
0130 
0131     mod_set = dsVary2Modifications(vary);
0132     % The first 2 cols of modifications_set are idenitical to vary, it just
0133     % has the last column distributed out to the number of sims
0134 
0135 
0136     % Get param names
0137     iMod = 1;
0138     % Split extra entries in first 2 cols of mods, so each row is a single pop and param
0139     [~, first_mod_set] = dsApplyModifications([],mod_set{iMod}, varargin{:});
0140 
0141     % replace '->' with '_'
0142     first_mod_set(:,1) = strrep(first_mod_set(:,1), '->', '_');
0143 
0144     % add col of underscores
0145     first_mod_set = cat(2,first_mod_set(:,1), repmat({'_'},size(first_mod_set,1), 1), first_mod_set(:,2:end));
0146     nParamMods = size(first_mod_set, 1);
0147 
0148     % get param names
0149     mod_params = cell(nParamMods,1);
0150     for iRow = 1:nParamMods
0151       mod_params{iRow} = [first_mod_set{iRow,1:3}];
0152 
0153       %check if variable in namespace
0154       if ~any(strcmp(model.namespaces(:,2), mod_params{iRow}))
0155         % find correct entry based on param and pop
0156         nsInd = logical(~cellfun(@isempty, strfind(model.namespaces(:,2), [first_mod_set{iRow,1} '_'])) .* ...
0157           ~cellfun(@isempty, strfind(model.namespaces(:,2), first_mod_set{iRow,3})));
0158 
0159         assert(sum(nsInd) == 1)
0160 
0161         % add mech names using namespace
0162         mod_params{iRow} = model.namespaces{nsInd,2};
0163       end
0164     end
0165 
0166     % Get param values for each sim
0167     param_values = nan(nParamMods, length(mod_set));
0168     for iMod = 1:length(mod_set)
0169       % Split extra entries in first 2 cols of mods, so each row is a single pop and param
0170       [~, mod_set{iMod}] = dsApplyModifications([],mod_set{iMod}, varargin{:});
0171 
0172       % Get scalar values as vector
0173       param_values(:, iMod) = [mod_set{iMod}{:,3}];
0174     end
0175 
0176     % Assign value vectors to params
0177     for iParam = 1:nParamMods
0178       p.(mod_params{iParam}) = param_values(iParam,:);
0179     end
0180   end % one_solve_file_flag
0181 
0182   if options.verbose_flag
0183     fprintf('saving params.mat\n');
0184   end
0185   save(param_file_name,'p','-v7');
0186 else
0187   % insert parameter values into model expressions
0188   model=dsPropagateParameters(model,'action','substitute', varargin{:});
0189 end
0190 
0191 % 1.2 prepare list of outputs (state variables and monitors)
0192 tmp=cellfun(@(x)[x ','],model.state_variables,'uni',0);
0193 tmp=[tmp{:}];
0194 output_string=tmp(1:end-1);
0195 
0196 if ~isempty(model.monitors)
0197   tmp=cellfun(@(x)[x ','],fieldnames(model.monitors),'uni',0);
0198   tmp=[tmp{:}];
0199   output_string=[output_string ',' tmp(1:end-1)];
0200 end
0201 
0202 if ~isempty(model.fixed_variables)
0203   tmp=cellfun(@(x)[x ','],fieldnames(model.fixed_variables),'uni',0);
0204   tmp=[tmp{:}];
0205   output_string=[output_string ',' tmp(1:end-1)];
0206 end
0207 
0208 output_string=['[T,' output_string ']']; % state vars, monitors, time vector
0209 
0210 % HACK to get IC to work
0211 if options.downsample_factor == 1
0212   for fieldNameC = fieldnames(model.ICs)'
0213     model.ICs.(fieldNameC{1}) = regexprep(model.ICs.(fieldNameC{1}), '_t0', '(1,:)');
0214   end
0215 end
0216 
0217 
0218 %% 3.0 write m-file solver
0219 % 2.1 create mfile
0220 if ~isempty(options.filename)
0221   if options.verbose_flag
0222     fprintf('Creating solver file: %s\n',options.filename);
0223   end
0224 
0225   fid=fopen(options.filename,'wt');
0226 else
0227   fid=options.fileID;
0228 end
0229 
0230 % get abs file path
0231 solve_ode_filepath = fopen(fid);
0232 
0233 if ~options.one_solve_file_flag
0234   fprintf(fid,'function %s=solve_ode\n',output_string);
0235 else
0236   fprintf(fid,'function %s=solve_ode(simID)\n',output_string);
0237 end
0238 
0239 % Benchmark tic
0240 if options.benchmark_flag
0241   fprintf(fid, 'tic;');
0242 end
0243 
0244 % 2.3 load parameters
0245 if options.save_parameters_flag
0246   fprintf(fid,'%% ------------------------------------------------------------\n');
0247   fprintf(fid,'%% Parameters:\n');
0248   fprintf(fid,'%% ------------------------------------------------------------\n');
0249   fprintf(fid,'params = load(''params.mat'',''p'');\n');
0250 
0251   if options.one_solve_file_flag && options.compile_flag
0252     fprintf(fid,'pVecs = params.p;\n');
0253   else
0254      fprintf(fid,'p = params.p;\n');
0255   end
0256 end
0257 
0258 if options.one_solve_file_flag
0259   % loop through p and for any vector, take simID index of it (ignores tspan)
0260   if ~options.compile_flag
0261     fprintf(fid,'\n%% For vector params, select index for this simID\n');
0262     fprintf(fid,'flds = fields(rmfield(p,''tspan''));\n'); % remove tspan
0263     fprintf(fid,'for fld = flds''\n');
0264     fprintf(fid,'  fld = fld{1};\n');
0265     fprintf(fid,'  if isnumeric(p.(fld)) && length(p.(fld)) > 1\n');
0266     fprintf(fid,'    p.(fld) = p.(fld)(simID);\n');
0267     fprintf(fid,'  end\n');
0268     fprintf(fid,'end\n\n');
0269   else %compile_flag
0270     % slice scalar from vector params
0271     for iParam = 1:nParamMods
0272       fprintf(fid,'p.%s = pVecs.%s(simID);\n', mod_params{iParam}, mod_params{iParam});
0273     end
0274 
0275     % take scalar from scalar params
0276     [~,sharedFlds] = intersect(fields(p), mod_params);
0277     scalar_params = fields(p);
0278     scalar_params(sharedFlds) = [];
0279     nScalarParams = length(scalar_params);
0280     for iParam = 1:nScalarParams
0281       fprintf(fid,'p.%s = pVecs.%s;\n', scalar_params{iParam}, scalar_params{iParam});
0282     end
0283   end
0284 end
0285 
0286 % write tspan, dt, and downsample_factor
0287 if options.save_parameters_flag
0288   fprintf(fid,'downsample_factor = %sdownsample_factor;\n',parameter_prefix);
0289   fprintf(fid,'dt = %sdt;\n',parameter_prefix);
0290   fprintf(fid,'T = (%stspan(1):downsample_factor*dt:%stspan(2))'';\n',parameter_prefix,parameter_prefix);
0291 else
0292   fprintf(fid,'tspan=[%g %g];\ndt = %g;\ndownsample_factor = %g;\n',options.tspan,options.dt,options.downsample_factor);
0293   fprintf(fid,'T = (tspan(1):downsample_factor*dt:tspan(2))'';\n');
0294 end
0295   % NOTE: T is different here since we take into account downsampling
0296 
0297 % write calculation of time vector and derived parameters: ntime, nsamp
0298 % fprintf(fid,'ntime = length(T);\nnsamp = length(1:downsample_factor*dt:ntime);\n');
0299 
0300 % 2.4 evaluate fixed variables
0301 if ~isempty(model.fixed_variables)
0302   fprintf(fid,'\n%% ------------------------------------------------------------\n');
0303   fprintf(fid,'%% Fixed variables:\n');
0304   fprintf(fid,'%% ------------------------------------------------------------\n');
0305   
0306   % 2.2 set random seed
0307   fprintf(fid,'%% seed the random number generator\n');
0308   fprintf(fid,'%% seed the random number generator\n');
0309   if options.save_parameters_flag
0310     fprintf(fid,'%s(%srandom_seed);\n',rng_function,parameter_prefix);
0311   else
0312     if ischar(options.random_seed)
0313       fprintf(fid,'%s(''%s'');\n',rng_function,options.random_seed);
0314     elseif isnumeric(options.random_seed)
0315       fprintf(fid,'%s(%g);\n',rng_function,options.random_seed);
0316     end
0317   end
0318   
0319   names=fieldnames(model.fixed_variables);
0320   expressions=struct2cell(model.fixed_variables);
0321   for i=1:length(names)
0322     fprintf(fid,'%s = %s;\n',names{i},expressions{i});
0323   end
0324 end
0325 
0326 % 2.5 evaluate function handles
0327 if ~isempty(model.functions) && options.reduce_function_calls_flag==0
0328   fprintf(fid,'\n%% ------------------------------------------------------------\n');
0329   fprintf(fid,'%% Functions:\n');
0330   fprintf(fid,'%% ------------------------------------------------------------\n');
0331   names=fieldnames(model.functions);
0332   expressions=struct2cell(model.functions);
0333   for i=1:length(names)
0334     fprintf(fid,'%s = %s;\n',names{i},expressions{i});
0335   end
0336 end
0337 
0338 % 2.6 prepare storage
0339 fprintf(fid,'\n%% ------------------------------------------------------------\n');
0340 fprintf(fid,'%% Initial conditions:\n');
0341 fprintf(fid,'%% ------------------------------------------------------------\n');
0342 
0343 % 2.2 set random seed
0344 fprintf(fid,'%% seed the random number generator\n');
0345 fprintf(fid,'%% seed the random number generator\n');
0346 if options.save_parameters_flag
0347   fprintf(fid,'%s(%srandom_seed);\n',rng_function,parameter_prefix);
0348 else
0349   if ischar(options.random_seed)
0350     fprintf(fid,'%s(''%s'');\n',rng_function,options.random_seed);
0351   elseif isnumeric(options.random_seed)
0352     fprintf(fid,'%s(%g);\n',rng_function,options.random_seed);
0353   end
0354 end
0355 
0356 %% Numerical integration
0357 % write code to do numerical integration
0358 fprintf(fid,'%% ###########################################################\n');
0359 fprintf(fid,'%% Numerical integration:\n');
0360 fprintf(fid,'%% ###########################################################\n');
0361 
0362 if options.compile_flag && strcmp(options.solver_type,'matlab_no_mex')
0363   odefun_str_name = odefun_filename;
0364 
0365   if options.compile_flag
0366     odefun_str_name = [odefun_str_name '_mex']; % switch to mex version
0367   end
0368 else
0369   odefun_str_name = 'odefun';
0370 end
0371 
0372 if ~isempty(options.matlab_solver_options)
0373   fprintf(fid,'[T,data] = %s(@%s, T, p.ic, p.matlab_solver_options);\n', options.solver, odefun_str_name);
0374 else
0375   fprintf(fid,'[T,data] = %s(@%s, T, p.ic);\n', options.solver, odefun_str_name);
0376 end
0377 
0378 %% Get vars from odefun output
0379 fprintf(fid,'\n%%Extract linear data into original state variables\n');
0380 
0381 % evaluate ICs to get (# elems) per state var and set up generic state var X
0382 num_vars=length(model.state_variables);
0383 state_var_index=0;
0384 for i=1:num_vars
0385   var=model.state_variables{i};
0386 
0387   % check ICs for use of inital state_var value and put in proper starting value
0388   if regexp(model.ICs.(var), '_last')
0389     stateVars = regexp(model.ICs.(var), '([\w_]+)_last', 'tokens');
0390     model.ICs.(var) = regexprep(model.ICs.(var), '_last', '');
0391 
0392     for iSVar = 1:length(stateVars)
0393       thisSvar = stateVars{iSVar}{1};
0394       model.ICs.(var) = regexprep(model.ICs.(var), thisSvar, model.ICs.(thisSvar));
0395     end
0396   end
0397 
0398   % evaluate ICs to get (# elems) per state var
0399   num_elems=length(eval([model.ICs.(var) ';']));
0400 
0401   % set state var indices a variables for generic state vector X
0402   data_inds = state_var_index + [1,num_elems];
0403 
0404   assert(strcmp(elem_names{data_inds(1)}, var)) %current elem should be same as var
0405 
0406   fprintf(fid,'%s = data(:, %i:%i);\n', var, data_inds(1), data_inds(end));
0407 
0408   state_var_index = state_var_index + num_elems;
0409 end
0410 
0411 %% Monitors
0412 if ~isempty(model.monitors)
0413   fprintf(fid,'\n%%Calculate monitors from params and state vars\n');
0414   monitor_names = fields(model.monitors);
0415   for iMon = 1:length(monitor_names)
0416     thisMonName = monitor_names{iMon};
0417     thisMonFcn = regexp(model.functions.(thisMonName),'@\([a-zA-Z][\w,]*\)\s*(.*)','tokens','once');
0418     thisMonFcn = thisMonFcn{1};
0419     fprintf(fid,'%s = %s;\n', thisMonName, thisMonFcn);
0420   end
0421 end
0422 
0423 %% Benchmark toc
0424 if options.benchmark_flag
0425   fprintf(fid, 'fprintf(''Sim Time: %%g seconds\\n'', toc);');
0426 end
0427 
0428 %% end solve function
0429 fprintf(fid,'\nend\n\n');
0430 
0431 %% ODEFUN
0432 if options.compile_flag && strcmp(options.solver_type,'matlab_no_mex') % save ode function as separate m-file for mex compilation
0433   %open file
0434   odefun_filepath = fullfile(fpath, [odefun_filename fext]);
0435   odefun_fid = fopen(odefun_filepath,'wt');
0436 
0437   %write to file
0438   fprintf(odefun_fid,'function dydt = %s(t,X)\n', odefun_filename);
0439   fprintf(odefun_fid,['dydt = [\n\n' odefun '\n]'';\n']); % make row into col vector
0440   fprintf(odefun_fid,'end\n');
0441 
0442   %close file
0443   fclose(odefun_fid);
0444 
0445   %% mex compile odefun
0446   options.codegen_args = {0,IC};
0447   dsPrepareMEX(odefun_filepath, options);
0448 
0449 else % use subfunction
0450   fprintf(fid,'\n%% ###########################################################\n');
0451   fprintf(fid,'%% SUBFUNCTIONS\n');
0452   fprintf(fid,'%% ###########################################################\n\n');
0453 
0454   % make sub function (no shared variables with main function workspace for max performance)
0455   fprintf(fid,'function dydt = odefun(t,X)\n');
0456   fprintf(fid,['dydt = [\n\n' odefun '\n]'';\n']); % make row into col vector
0457   fprintf(fid,'end\n');
0458 end
0459 
0460 if ~strcmp(solve_ode_filepath,'"stdout"')
0461   fclose(fid);
0462   % wait for file before continuing to simulation
0463   while ~exist(solve_ode_filepath,'file')
0464     pause(.01);
0465   end
0466 end
0467 
0468 end %function
0469 %% END MAIN FUNC