Description of dsGenerateModel

0001 function [model,name_map] = dsGenerateModel(specification, varargin)
0002 %GENERATEMODEL - Parse DynaSim specification and organize model data in DynaSim model structure
0003 %
0004 % Usage:
0005 %   [model,name_map]=dsGenerateModel(specification,'option',value,...)
0006 %
0007 % Inputs:
0008 %   - specification: one of:
0009 %     - DynaSim specification structure (see below and dsCheckSpecification for more details)
0010 %     - string with name of MAT-file containing DynaSim specification structure
0011 %     - string with equations
0012 %     - string with name of file containing equations (.eqns)
0013 %       note: .eqns files can also be converted into model structure using LoadModel()
0014 %   - options (with defaults): 'option1',value1,'option2',value2,...
0015 %     'modifications'  : specify modifications to apply to specification
0016 %                        before generating the model, see dsApplyModifications
0017 %                        for more details (default?: []).
0018 %     'open_link_flag' : whether to leave linker identifiers in place (default: 0)
0019 %     'auto_gen_test_data_flag': whether to save model for unit testing (default: 0)
0020 %
0021 % Outputs:
0022 %   - model: DynaSim model structure (see dsCheckModel for more details):
0023 %     .parameters      : substructure with model parameters
0024 %     .fixed_variables : substructure with fixed variable definitions
0025 %     .functions       : substructure with function definitions
0026 %     .monitors        : substructure with monitor definitions
0027 %     .state_variables : cell array listing state variables
0028 %     .ODEs            : substructure with one ordinary differential
0029 %                             equation (ODE) per state variable
0030 %     .ICs             : substructure with initial conditions (ICs) for
0031 %                             each state variable
0032 %     .conditionals(i) : structure array with each element indicating
0033 %                             conditional actions specified in subfields
0034 %                             "condition","action","else" (see NOTE 1 in dsCheckModel)
0035 %     .linkers(i)      : structure array with each element indicating
0036 %                             an "expression" that should be inserted
0037 %                             (according to "operation") into any equations
0038 %                             where the "target" appears. (see NOTE 2 in dsCheckModel)
0039 %       .target    : string giving the target where expression should be inserted
0040 %       .expression: string giving the expression to insert
0041 %       .operation : string giving the operation to use to insert expression
0042 %     .comments{i}     : cell array of comments found in model files
0043 %     .specification   : specification used to generate the model
0044 %     .namespaces      : (see NOTE 3 in dsCheckModel)
0045 %   - name_map: cell matrix mapping parameter, variable, and function names
0046 %       between the user-created specification (population equations and mechanism
0047 %       files) and the full model with automatically generated namespaces. It
0048 %       has four columns with: user-specified name, name with namespace prefix,
0049 %       namespace, and type ('parameters', 'fixed_variables', 'state_variables',
0050 %       'functions', or 'monitors') indicating the category to which the named
0051 %       element belongs.
0052 %
0053 % - DynaSim specification structure (see dsCheckSpecification for more details)
0054 %   .populations(i) (required): contains info for defining independent population models
0055 %       .name (default: 'pop1')      : name of population
0056 %       .size (default: 1)           : number of elements in population (i.e., # cells)
0057 %       .equations (required)        : string listing equations (see NOTE 1 in dsCheckSpecification)
0058 %       .mechanism_list (default: []): cell array listing mechanisms (see NOTE 2
0059 %                                      in dsCheckSpecification)
0060 %       .parameters (default: [])    : parameters to assign across all equations in
0061 %         the population. provide as cell array list of key/value pairs
0062 %         {'param1',value1,'param2',value2,...}
0063 %       .model (default: [])   : optional DynaSim model structure
0064 %   .connections(i) (default: []): contains info for linking population models
0065 %       .source (required if >1 pops): name of OUT population
0066 %       .target (required if >1 pops): name of target population
0067 %       .mechanism_list (required)   : list of mechanisms that link two populations
0068 %       .parameters (default: [])    : parameters to assign across all equations in
0069 %         mechanisms of this connection's mechanism_list.
0070 %
0071 % Examples:
0072 %   - Example 0:
0073 %     model=dsGenerateModel('db/dt=3')
0074 %
0075 %   - Example 1: Lorenz equations
0076 %     eqns={
0077 %       's=10; r=27; b=2.666';
0078 %       'dx/dt=s*(y-x)';
0079 %       'dy/dt=r*x-y-x*z';
0080 %       'dz/dt=-b*z+x*y';
0081 %     };
0082 %     model=dsGenerateModel(eqns)
0083 %
0084 %   - Example 2: Leaky integrate-and-fire neuron
0085 %     model=dsGenerateModel('tau=10; R=10; E=-70; dV/dt=(E-V+R*1.55)/tau; if(V>-55)(V=-75)')
0086 %
0087 %   - Example 3: Hodgkin-Huxley neuron with sinusoidal drive
0088 %     model=dsGenerateModel('dv/dt=current+sin(2*pi*t); {iNa,iK}')
0089 %
0090 %   - Example 4: HH with self inhibition and sinusoidal drive
0091 %     specification.populations(1).equations='dv/dt=current+sin(2*pi); v(0)=-65';
0092 %     specification.populations(1).mechanism_list={'iNa','iK'};
0093 %     specification.connections(1).mechanism_list={'iGABAa'};
0094 %     specification.connections(1).parameters={'tauDx',15};
0095 %     model=dsGenerateModel(specification)
0096 %
0097 %   - Example 5: using custom mechanism alias in equations (for modularization)
0098 %     specification.populations(1).equations='dv/dt=@M+sin(2*pi); v(0)=-65';
0099 %     specification.populations(1).mechanism_list={'iNa@M','iK@M'};
0100 %     model=dsGenerateModel(specification)
0101 %
0102 %     or:
0103 %
0104 %     specification.populations(1).equations='dv/dt=@M+sin(2*pi); {iNa,iK}@M; v(0)=-65';
0105 %     model=dsGenerateModel(specification)
0106 %
0107 %   - Example 6: directly incorporating mechanism models from online repositories:
0108 %
0109 %     model=dsGenerateModel('dv/dt=@M; {ib:57,iK}@M')
0110 %
0111 %     where "ib" is a known alias for the infinitebrain.org repository,
0112 %     and "57" is the Na+ current at http://infinitebrain.org/models/57.
0113 %     note: currently not supported on *most* machines...
0114 %
0115 % See also: dsCheckSpecification, dsCheckModel, dsParseModelEquations, dsSimulate
0116 %
0117 % Author: Jason Sherfey, PhD <jssherfey@gmail.com>
0118 % Copyright (C) 2016 Jason Sherfey, Boston University, USA
0119 
0120 % Check inputs
0121 % ------------------------------------------------------
0122 if nargin==0
0123   % use default model
0124   specification=[];
0125   specification.populations(1).equations='dv/dt=10+@current/Cm; Cm=1; v(0)=-65';
0126   specification.populations(1).mechanism_list={'iNa','iK'};
0127   specification.populations(1).parameters={'Cm',1};
0128   specification.connections(1).mechanism_list={'iGABAa'};
0129   varargin={'modifications',[]};
0130 end
0131 % ------------------------------------------------------
0132 
0133 options=dsCheckOptions(varargin,{...
0134   'modifications',[],[],...
0135   'open_link_flag',0,{0,1},...
0136   'auto_gen_test_data_flag',0,{0,1},...
0137   },false);
0138 
0139 if options.auto_gen_test_data_flag
0140   varargs = varargin;
0141   varargs{find(strcmp(varargs, 'auto_gen_test_data_flag'))+1} = 0;
0142   varargs(end+1:end+2) = {'unit_test_flag',1};
0143   argin = [{specification}, varargs]; % specific to this function
0144 end
0145 
0146 % check if a model
0147 if isfield(specification,'state_variables')
0148   % do nothing
0149   model=specification;
0150   return;
0151 %   TODO: consider the following --
0152 %   if isfield(specification,'specification')
0153 %     % regenerate from specification
0154 %     specification=specification.specification;
0155 %   end
0156 end
0157 % standardize specification
0158 specification=dsCheckSpecification(specification, varargin{:}); % standardize & auto-populate as needed
0159 
0160 % Apply modifications to specification before generating model
0161 if ~isempty(options.modifications)
0162   specification=dsApplyModifications(specification,options.modifications, varargin{:});
0163 end
0164 
0165 % specification metadata:
0166 npops=length(specification.populations); % number of populations
0167 ncons=length(specification.connections); % number of connections
0168 
0169 %{
0170 % Dev notes on improving implementation:
0171 % Ideally (1.0)-(3.0) could be packaged into external functions and run as:
0172 % -------------------------------------------------------------------------
0173 %% 1.0 load sub-models, assign namespaces, and combine across all equations and mechanisms in specification
0174 % [model,name_map]=LoadModelSet(specification) % bug: disrupted subsequent namespace propagation (without raising an error)
0175 %% 2.0 propagate namespaces through variable and function names
0176 % model = dsPropagateNamespaces(model,name_map); % this works
0177 %% 3.0 expand population equations according to mechanism linkers
0178 % model = LinkMechanisms(model,name_map);      % problem: unable to identify linker population from model.linkers; see notes below (3.0) for more details
0179 % -------------------------------------------------------------------------
0180 %}
0181 
0182 % support full modularization of mechanisms
0183 % (eg, dv/dt=@M; {Na,K}@M w/ Na.mech: @current += I(IN,m,h)).
0184 %     approach taken below:
0185 %     - add support for dv/dt=@M; {Na@M,K@M}
0186 %       have dsGenerateModel split mech_name on '@' and replace first
0187 %       linker in mech (e.g., @current) by what follows '@' (e.g., @M)
0188 %     - then have dsCheckSpecification convert {Na,K}@M into {Na@M,K@M}
0189 
0190 %% 1.0 load sub-models, assign namespaces, and combine across all equations and mechanisms in specification
0191 % use empty struct for Octave compatibility
0192 model.parameters=struct('');
0193 model.fixed_variables=struct('');
0194 model.functions=struct('');
0195 model.monitors=struct('');
0196 model.state_variables={};
0197 model.ODEs=struct('');
0198 model.ICs=struct('');
0199 model.conditionals=struct('');
0200 model.linkers=struct('');
0201 model.comments={};
0202 name_map={}; % {name, namespace_name, namespace, type}, used for namespacing
0203 linker_pops={}; % list of populations associated with mechanism linkers
0204 
0205 % 1.1 load and combine population sub-models from population equations and mechanisms
0206 for i=1:npops
0207   % does the population model already exist?
0208   if ~isempty(specification.populations(i).model)
0209     tmpmodel=specification.populations(i).model; % get model structure
0210     tmpname=tmpmodel.specification.populations.name; % assumes one population sub-model
0211 
0212     % adjust the name if necessary
0213     if ~strcmp(specification.populations(i).name,tmpname)
0214       % use the name in the specification
0215       tmpmodel=dsApplyModifications(tmpmodel,{tmpname,'name',specification.populations(i).name}, varargin{:});
0216     elseif strcmp(tmpname,'pop1') % if default name
0217       % use default name for this population index
0218       tmpmodel=dsApplyModifications(tmpmodel,{tmpname,'name',sprintf('pop%g',i)}, varargin{:});
0219     end
0220 
0221     tmpmodel.linkers=[]; % remove old linkers from original model construction
0222     model=dsCombineModels(model,tmpmodel, varargin{:});
0223     name_map=cat(1,name_map,tmpmodel.namespaces);
0224     continue;
0225   end
0226   % construct new population model
0227   PopScope=specification.populations(i).name;
0228     % NOTE: dsParseModelEquations adds a '_' suffix to the namespace; therefore,
0229     % a '_' suffix is added to PopScope when used below for consistency of
0230     % namespaces/namespaces. (this could be cleaned up by adding '_' to PopScope
0231     % here, removing it below, and removing the additional '_' from
0232     % dsParseModelEquations).
0233 
0234   % 1.1.1 parse population equations
0235   equations=specification.populations(i).equations;
0236   parameters=specification.populations(i).parameters;
0237   nmechs=length(specification.populations(i).mechanism_list);
0238 
0239   % parse population equations
0240   if ~isempty(equations)
0241     [tmpmodel,tmpmap]=dsImportModel(equations,'namespace',PopScope,'ic_pop',specification.populations(i).name,'user_parameters',parameters);
0242     model=dsCombineModels(model,tmpmodel, varargin{:});
0243     name_map=cat(1,name_map,tmpmap);
0244   end
0245 
0246   % 1.1.2 parse population mechanisms
0247   for j=1:nmechs
0248     mechanism_=specification.populations(i).mechanism_list{j};
0249 
0250     % support separation of linker names in pop equations vs mechanisms
0251     mechanism_=regexp(mechanism_,'@','split');
0252     mechanism=mechanism_{1};
0253 
0254     if numel(mechanism_)>1, new_linker=mechanism_{2}; else new_linker=[]; end
0255 
0256     % set mechanism namespace
0257     [~,MechID]=fileparts2(mechanism);
0258     if any(MechID==':')
0259       % exclude host name from namespace
0260       tmp=regexp(MechID,':','split');
0261       MechScope=[specification.populations(i).name '_' tmp{2}];
0262     else
0263       % extract mechanism file name without path
0264       MechScope=[specification.populations(i).name '_' MechID];
0265     end
0266     % use mechanism equations in specification if present
0267     if isfield(specification.populations,'mechanisms') && ~isempty(specification.populations(i).mechanisms)
0268       if ismember(MechID,{specification.populations(i).mechanisms.name})
0269         idx=ismember({specification.populations(i).mechanisms.name},MechID);
0270         mechanism=specification.populations(i).mechanisms(idx).equations;
0271       end
0272       % parse mechanism equations
0273       if ~isempty(mechanism)
0274         [tmpmodel,tmpmap]=dsImportModel(mechanism,'namespace',MechScope,'ic_pop',specification.populations(i).name,'user_parameters',parameters);
0275         % replace 1st linker name by the one in specification
0276         if ~isempty(new_linker) && ~isempty(tmpmodel.linkers)
0277           % first try to find 1st linker target starting with @
0278           links_at=find(~cellfun(@isempty,regexp({tmpmodel.linkers.target},'^@','once')));
0279           if ~isempty(links_at)
0280             % use first link with target prepended by '@'
0281             link_ind=links_at(1);
0282           else
0283             % use first link
0284             link_ind=1;
0285           end
0286           tmpmodel.linkers(link_ind).target=['@' new_linker];
0287         end
0288         % combine sub-model with other sub-models
0289         model=dsCombineModels(model,tmpmodel, varargin{:});
0290         name_map=cat(1,name_map,tmpmap);
0291         linker_pops=cat(2,linker_pops,repmat({specification.populations(i).name},[1 length(tmpmodel.linkers)]));
0292       end
0293     end
0294   end
0295   pop=specification.populations(i).name;
0296 
0297   % add reserved keywords (parameters and state variables) to name_map
0298   add_keywords(pop,pop,[PopScope '_']);
0299   %model.parameters.([pop '_Npop'])=num2str(specification.populations(i).size);
0300   model.parameters(1).([pop '_Npop'])=toString(specification.populations(i).size,0);
0301 end
0302 
0303 % 1.2 load and combine sub-models from connection mechanisms
0304 for i=1:ncons
0305   % parse connection mechanisms
0306   source=specification.connections(i).source;
0307   target=specification.connections(i).target;
0308   parameters=specification.connections(i).parameters;
0309   ConScope=[target '_' source '_'];
0310     % NOTE: in contrast to PopScope above, ConScope is never passed to
0311     % dsParseModelEquations; thus the '_' should be added here for consistency
0312     % with mechanism namespaces (which are modified by dsParseModelEquations).
0313   for j=1:length(specification.connections(i).mechanism_list)
0314     mechanism_=specification.connections(i).mechanism_list{j};
0315 
0316     % support separation of linker names in pop equations vs mechanisms
0317     mechanism_=regexp(mechanism_,'@','split');
0318     mechanism=mechanism_{1};
0319 
0320     if numel(mechanism_)>1, new_linker=mechanism_{2}; else new_linker=[]; end
0321 
0322     % extract mechanism file name without path
0323     [~,MechID]=fileparts2(mechanism);
0324     MechScope=[target '_' source '_' MechID];
0325         % NOTE: must use target_source_mechanism for connection mechanisms
0326         % to distinguish their parent namespaces from those of population mechanisms
0327         % see: dsGetParentNamespace
0328 
0329     % use mechanism equations in specification if present
0330     if ~isempty(specification.connections(i).mechanisms) && ismember(MechID,{specification.connections(i).mechanisms.name})
0331       idx=ismember({specification.connections(i).mechanisms.name},MechID);
0332       mechanism=specification.connections(i).mechanisms(idx).equations;
0333     end
0334 
0335     % parse model equations
0336     [tmpmodel,tmpmap]=dsImportModel(mechanism,'namespace',MechScope,'ic_pop',source,'user_parameters',parameters);
0337 
0338     % replace 1st linker name by the one in specification
0339     if ~isempty(new_linker) && ~isempty(tmpmodel.linkers)
0340       tmpmodel.linkers(1).target=['@' new_linker];
0341     end
0342 
0343     model=dsCombineModels(model,tmpmodel, varargin{:});
0344     name_map=cat(1,name_map,tmpmap);
0345 
0346     % link this mechanism to the target population
0347     linker_pops=cat(2,linker_pops,repmat(target,[1 length(tmpmodel.linkers)]));
0348 
0349     % edit names of connection monitors specified in population equations
0350     % TODO: consider design changes to avoid specifying connection monitors
0351     %   in population equations; this is an undesirable hack:
0352     %     eg, convert E_iGABAa_functions -> I_E_iGABAa_functions
0353     if ~isempty(model.monitors)
0354       % get indices to all model.monitors that have incorrect connection namespace
0355       con_mon_to_update=find(~cellfun(@isempty,regexp(fieldnames(model.monitors),['^' target '_' mechanism])));
0356       if any(con_mon_to_update)
0357         % get list of current model.monitors
0358         monitor_names=fieldnames(model.monitors);
0359         for m=1:length(con_mon_to_update)
0360           % get name of monitor with incorrect connection namespace
0361           old=monitor_names{con_mon_to_update(m)};
0362 
0363           % get name of monitor with correct connection namespace
0364           new=strrep(old,[target '_' mechanism '_'],[MechScope '_']);
0365 
0366           % add new monitor with correct namespace
0367           model.monitors.(new)=model.monitors.(old);
0368 
0369           % remove old monitor with incorrect namespace
0370           model.monitors=rmfield(model.monitors,old);
0371         end
0372       end
0373     end
0374   end
0375   % add reserved keywords (parameters and state variables) to name_map
0376   add_keywords(source,target,ConScope);
0377 end
0378 
0379 % check for monitoring functions (e.g., 'monitor functions' or 'monitor Na.functions')
0380 if ~isempty(model.monitors)
0381   % get list of functions
0382   if ~isempty(model.functions)
0383     function_names=fieldnames(model.functions);
0384   else
0385     function_names={};
0386   end
0387   % get list of monitor names
0388   monitor_names=fieldnames(model.monitors);
0389 
0390   % get indices to monitors with names ending in _functions
0391   function_monitor_index=find(~cellfun(@isempty,regexp(monitor_names,'_functions$','once')));
0392 
0393   % create list of functions with namespaces matching monitors ending in _functions
0394   functions_to_monitor={};
0395 
0396   for i=1:length(function_monitor_index)
0397     % get namespace of functions to monitor
0398     monitor_name=monitor_names{function_monitor_index(i)};
0399     monitor_namespace=regexp(monitor_name,'(.*)_functions$','tokens','once');
0400     monitor_namespace=monitor_namespace{1};
0401 
0402     % get list of functions with matching namespace
0403     function_index=find(~cellfun(@isempty,regexp(function_names,['^' monitor_namespace],'once')));
0404 
0405     % add functions to list
0406     functions_to_monitor=cat(1,functions_to_monitor,function_names(function_index));
0407 
0408     % remove "function" monitor name
0409     model.monitors=rmfield(model.monitors,monitor_name);
0410   end
0411 
0412   % eliminate duplicate function names
0413   functions_to_monitor=unique_wrapper(functions_to_monitor);
0414 
0415   % add functions to monitor list
0416   for i=1:length(functions_to_monitor)
0417     model.monitors.(functions_to_monitor{i})=[];
0418   end
0419 end
0420 
0421   % ----------------------------------
0422   % NESTED FUNCTIONS
0423   % ----------------------------------
0424   function add_keywords(src,dst,namespace)
0425     % NOTE: this needs to be coordinated with update_keywords() in dsSimulate()
0426     %   for parameters
0427     Nsrc=[src '_Npop'];
0428     Ndst=[dst '_Npop'];
0429 
0430     old={'Npre','N[1]','N_pre','Npost','N_post','N[0]','Npop','N_pop','tspike_pre','tspike_post','tspike'};
0431     new={Nsrc,Nsrc,Nsrc,Ndst,Ndst,Ndst,Ndst,Ndst,[src '_tspike'],[dst '_tspike'],[dst '_tspike']};
0432     for p=1:length(old)
0433       name_map(end+1,:)={old{p},new{p},namespace,'parameters'};
0434     end
0435 
0436     % for state variables
0437     new={};
0438     old={};
0439     src_excluded=~cellfun(@isempty,regexp(name_map(:,1),['pre' '$']));
0440     dst_excluded=~cellfun(@isempty,regexp(name_map(:,1),['post' '$']));
0441     excluded=src_excluded|dst_excluded;
0442 
0443     PopScope=[src '_'];
0444     var_idx=strcmp(PopScope,name_map(:,3)) & strcmp('state_variables',name_map(:,4)) & ~excluded;
0445     if any(var_idx)
0446       Xsrc_old_vars=name_map(var_idx,1);
0447       Xsrc_new_vars=name_map(var_idx,2);
0448       % default for IN is first Xsrc state var
0449       Xsrc=Xsrc_new_vars{1};
0450       old=cat(2,old,{'IN','Xpre','X_pre'});
0451       new=cat(2,new,{Xsrc,Xsrc,Xsrc});
0452     else
0453       Xsrc_old_vars=[];
0454       Xsrc_new_vars=[];
0455       Xsrc=[];
0456     end
0457 
0458     PopScope=[dst '_'];
0459     var_idx=strcmp(PopScope,name_map(:,3)) & strcmp('state_variables',name_map(:,4)) & ~excluded;
0460     if any(var_idx)
0461       Xdst_old_vars=name_map(var_idx,1);
0462       Xdst_new_vars=name_map(var_idx,2);
0463       % default for OUT and X is first Xdst state var
0464       Xdst=Xdst_new_vars{1};
0465       old=cat(2,old,{'OUT','X','Xpost','X_post'});
0466       new=cat(2,new,{Xdst,Xdst,Xdst,Xdst});
0467     else
0468       Xdst_old_vars=[];
0469       Xdst_new_vars=[];
0470       Xdst=[];
0471     end
0472 
0473     % add variants [var_pre,var_post,varpre,varpost]
0474     if ~isempty(Xsrc_old_vars)
0475       [Xsrc_old_vars,IA]=setdiff(Xsrc_old_vars,old);
0476       Xsrc_new_vars=Xsrc_new_vars(IA);
0477     end
0478 
0479     if ~isempty(Xdst_old_vars)
0480       [Xdst_old_vars,IA]=setdiff(Xdst_old_vars,old);
0481       Xdst_new_vars=Xdst_new_vars(IA);
0482     end
0483 
0484     for p=1:length(Xsrc_old_vars)
0485       old{end+1}=[Xsrc_old_vars{p} '_pre'];
0486       new{end+1}=Xsrc_new_vars{p};
0487       old{end+1}=[Xsrc_old_vars{p} 'pre'];
0488       new{end+1}=Xsrc_new_vars{p};
0489     end
0490 
0491     for p=1:length(Xdst_old_vars)
0492       old{end+1}=[Xdst_old_vars{p} '_post'];
0493       new{end+1}=Xdst_new_vars{p};
0494       old{end+1}=[Xdst_old_vars{p} 'post'];
0495       new{end+1}=Xdst_new_vars{p};
0496     end
0497 
0498     for p=1:length(old)
0499       name_map(end+1,:)={old{p},new{p},namespace,'state_variables'};
0500     end
0501   end
0502 
0503 %% 2.0 propagate namespaces through variable and function names
0504 %      i.e., to establish uniqueness of names by adding namespace/namespace prefixes)
0505 model.specification=specification;
0506 model = dsPropagateNamespaces(model,name_map, varargin{:});
0507 
0508 %% 3.0 expand population equations according to mechanism linkers
0509 % purpose: expand population equations according to linkers
0510 % - link populations.equations to mechanism sub-models
0511 % - link mechanism functions and state variables across mechanisms in a given population
0512 
0513 % store indices to all expressions and conditionals that are linked (this
0514 % is necessary for efficiently removing linker targets from expressions after linking)
0515 all_expression_inds=[];
0516 all_expression_targets={};
0517 all_conditionals_inds=[];
0518 all_conditionals_targets={};
0519 
0520 % add variables to linked expression if its a function without ()
0521 if ~isempty(model.functions) && ~isempty(model.linkers)
0522   function_names=fieldnames(model.functions);
0523   expressions={model.linkers.expression};
0524   [~,I,J]=intersect(function_names,expressions);
0525   for i=1:length(I)
0526     e=model.functions.(function_names{I(i)}); % function expression (eg,'@(x,y,z)x-(y-z)')
0527     v=regexp(e,'@(\([\w,]+\))','tokens','once'); % function input list (eg, '(x,y,z)')
0528     if ~isempty(v)
0529       model.linkers(J(i)).expression=[model.linkers(J(i)).expression v{1}];
0530     end
0531   end
0532 end
0533 
0534 % loop over linkers
0535 for i=1:length(model.linkers)
0536   % determine how to link
0537   operation=model.linkers(i).operation;
0538   oldstr=model.linkers(i).target;
0539   newstr=model.linkers(i).expression;
0540   switch operation % see dsClassifyEquation and dsParseModelEquations   % ('((\+=)|(-=)|(\*=)|(/=)|(=>))')
0541     case '+='
0542       operator='+';
0543     case '-='
0544       operator='-';
0545     case '*='
0546       operator='.*';
0547     case '/='
0548       operator='./';
0549     otherwise
0550       operator='+';
0551   end
0552   % determine what to link (ie, link across everything belonging to the linker population)
0553   % explicitly constrain to linker population
0554   expressions_in_pop=~cellfun(@isempty,regexp(name_map(:,3),['^' linker_pops{i}]));
0555 
0556   if ~isempty(model.conditionals)
0557     conditionals_in_pop=~cellfun(@isempty,regexp({model.conditionals.namespace},['^' linker_pops{i}]));
0558   end
0559 
0560   if ~isempty(model.linkers)
0561     linkers_in_pop=~cellfun(@isempty,regexp({model.linkers.namespace},['^' linker_pops{i}]));
0562   end
0563 
0564   % constrain to namespace
0565   names_in_namespace=cellfun(@(x,y)strncmp(y,x,length(y)),name_map(:,2),name_map(:,3));
0566 
0567   % get list of (functions,monitors,ODEs) belonging to the linker population
0568   eqn_types={'ODEs','monitors','functions'};%{'monitors','ODEs'};
0569   search_types={'state_variables','monitors','functions'};%{'monitors','state_variables'};
0570 
0571   % indices to expressions in the linker population with the correct search_types and namespace
0572   inds=find(expressions_in_pop&names_in_namespace&ismember(name_map(:,4),search_types));
0573 
0574   % eliminate duplicates (e.g., state_variables replacing OUT and X)
0575   [jnk,ia,ib]=unique_wrapper(name_map(inds,2),'stable');
0576   inds=inds(ia);
0577   all_expression_inds=[all_expression_inds inds'];
0578   all_expression_targets=cat(2,all_expression_targets,repmat({oldstr},[1 length(inds)]));
0579 
0580   % substitute link
0581   for j=1:length(inds)
0582     name=name_map{inds(j),2}; % name of variable as stored in model structure
0583     type=name_map{inds(j),4}; % search_types
0584     eqn_type=eqn_types{strcmp(type,search_types)}; % corresponding equation type
0585 
0586     % update expression with the current link
0587     if isfield(model.(eqn_type),name)
0588       % note: name will not be a field of eqn_type for special monitors
0589       % (e.g., monitor functions)
0590       model.(eqn_type).(name)=linker_strrep(model.(eqn_type).(name),oldstr,newstr,operator);
0591     end
0592   end
0593 
0594   if ~isempty(model.conditionals)
0595     fields={'condition','action','else'};
0596 
0597     % get list of conditionals belonging to the linker population
0598     inds=find(conditionals_in_pop);
0599     all_conditionals_inds=[all_conditionals_inds inds];
0600     all_conditionals_targets=cat(2,all_conditionals_targets,repmat({oldstr},[1 length(inds)]));
0601 
0602     % substitute link
0603     for j=1:length(inds)
0604       for field_index=1:length(fields)
0605         field=fields{field_index};
0606         model.conditionals(inds(j)).(field)=linker_strrep(model.conditionals(inds(j)).(field),oldstr,newstr,operator);
0607       end
0608     end
0609   end
0610 
0611   if ~isempty(model.linkers)
0612     inds=find(linkers_in_pop);
0613     for j=1:length(inds)
0614       model.linkers(inds(j)).expression=linker_strrep(model.linkers(inds(j)).expression,oldstr,newstr,operator);
0615     end
0616   end
0617 end
0618 
0619 if options.open_link_flag==0
0620   % remove target placeholders from expressions and conditionals
0621   for i=1:length(all_expression_inds)
0622     oldstr=all_expression_targets{i};
0623     newstr='';
0624     name=name_map{all_expression_inds(i),2};
0625     type=name_map{all_expression_inds(i),4};
0626     eqn_type=eqn_types{strcmp(type,search_types)};
0627     pattern = ['\)\.?[-\+\*/]' oldstr '\)']; % pattern accounts for all possible newstr defined for linking
0628     replace = [newstr '))'];
0629     if isfield(model.(eqn_type),name) && ischar(model.(eqn_type).(name))
0630         % NOTE: name will not be a field of eqn_type for special monitors
0631         % (e.g., monitor functions)
0632       model.(eqn_type).(name)=regexprep(model.(eqn_type).(name),pattern,replace);
0633     end
0634   end
0635 end
0636 
0637 if ~isempty(model.conditionals)
0638   for i=1:length(all_conditionals_inds)
0639     oldstr=all_conditionals_targets{i};
0640     newstr='';
0641     pattern = ['\)\.?[-\+\*/]' oldstr '\)']; % pattern accounts for all possible newstr defined for linking
0642     replace = [newstr '))'];
0643     for field_index=1:length(fields)
0644       field=fields{field_index};
0645       if model.conditionals(all_conditionals_inds(i)).(field)
0646         model.conditionals(all_conditionals_inds(i)).(field)=regexprep(model.conditionals(all_conditionals_inds(i)).(field),pattern,replace);
0647       end
0648     end
0649   end
0650 end
0651 
0652 % ------------------------------------------
0653 % NOTE on non-ideal implementation of 3.0: model.linkers does not contain enough
0654 % information to determine the population to which it belongs in all cases
0655 % (due to namespace format differences for population vs connection mechanisms & models
0656 % with one vs more populations). consequently, had to perform linking in this
0657 % function using info stored above while parsing the model; ideally, the
0658 % linking could occur independently of this function, informed by info in
0659 % model.linkers, and be packaged in its own external function LinkMechanisms().
0660 % ------------------------------------------
0661 
0662 %% 4.0 finalize
0663 
0664 % 4.1 sort .ODEs and .ICs wrt .state_variables
0665 if ~isempty(model.ODEs)
0666   model.ODEs = orderfields(model.ODEs,model.state_variables);
0667   model.ICs = orderfields(model.ICs,model.state_variables);
0668 end
0669 
0670 % 4.2 convert to numeric parameters
0671 c = struct2cell(model.parameters);
0672 
0673 % get index of strings
0674 idx1=find(cellfun(@ischar,c));
0675 
0676 % which strings contain numeric values?
0677 idx2=find(cellfun(@isempty,regexp(c(idx1),'[a-z_A-Z]')) | ~cellfun(@isempty,regexp(c(idx1),'^\s*\[*\s*\+?inf\s*\]*\s*$','ignorecase')));
0678 
0679 % convert those strings which contain numeric values
0680 c(idx1(idx2)) = cellfun(@eval,c(idx1(idx2)),'uni',0);
0681 f = fieldnames(model.parameters);
0682 model.parameters = cell2struct(c,f,1);
0683 
0684 % 4.3 store original specification
0685 model.specification = specification; % store specification to enable modifications to be applied later
0686 model.namespaces = name_map; % store name_map for transparency
0687 
0688 model=dsCheckModel(model, varargin{:});
0689 
0690 %% auto_gen_test_data_flag argout
0691 if options.auto_gen_test_data_flag
0692   argout = {model, name_map}; % specific to this function
0693 
0694   dsUnitSaveAutoGenTestData(argin, argout);
0695 end
0696 
0697 end % main function
0698 
0699 
0700 %% SUBFUNCTIONS
0701 function str=linker_strrep(str,oldstr,newstr,operator)
0702   if isempty(str)
0703     return;
0704   end
0705   % if inserting one word (e.g., a state variable), just replace it
0706   % WARNING: could cause problems in future if there is an additive
0707   % substitution of different state variables into the same place followed
0708   % by non-additive operations (e.g., @cai+=cai1 and @cai+=cai2 into
0709   % v'=f(v)*cai where cai1 & cai2 are defined in mechanisms for the same v;
0710   % workaround: insert into v'=f(v)*(cai)).
0711 
0712   % check if anything besides a single variable:
0713   if isempty(regexp(newstr,'[^a-z_A-Z\d]+','once'))
0714     str=dsStrrep(str,oldstr,newstr);
0715   else
0716     % otherwise do substitution with operator and parenthesis
0717     pat=['([^\w]+)' oldstr '([^\w]+)']; % in the middle
0718     rep=['$1((' newstr ')' operator oldstr ')$2'];
0719     str=regexprep(str,pat,rep);
0720     pat=['([^\w]+)' oldstr '$'];        % at the end
0721     rep=['$1((' newstr ')' operator oldstr ')'];
0722     str=regexprep(str,pat,rep);
0723     pat=['^' oldstr '([^\w]+)'];        % at the beginning
0724     rep=['((' newstr ')' operator oldstr ')$1'];
0725     str=regexprep(str,pat,rep);
0726     pat=['^' oldstr '$'];               % all there is
0727     rep=['((' newstr ')' operator oldstr ')'];
0728     str=regexprep(str,pat,rep);
0729   end
0730 end
dsGenerateModel

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
