Description of dsApplyModifications

0001 function [output,modifications] = dsApplyModifications(model, modifications, varargin)
0002 %APPLYMODIFICATIONS - Apply modifications to DynaSim specification or model structure
0003 %
0004 % dsApplyModifications returns the same kind of structure with
0005 % modifications applied. In all cases it first modifies the specification
0006 % (model.specification if input is a model structure). Then it returns
0007 % the modified specification or regenerates the model using the new
0008 % specification.
0009 %
0010 % Inputs:
0011 %   - model: DynaSim model or specification structure
0012 %     - see dsCheckModel and dsCheckSpecification for details
0013 %   - modifications: modifications to make to specification structure
0014 %       {X,Y,Z; X,Y,Z; ...}
0015 %       X = population name or connection source->target
0016 %       Y = thing to modify ('name', 'size', or parameter name)
0017 %       set Y=Z if Y = name, size, or value
0018 %       Note: (X1,X2) or (Y1,Y2): modify these simultaneously in the same way
0019 %       Note: Z can be a scalar, row vector, column vector, or matrix. Columns
0020 %       of Z are applied across items Y1, Y2, etc (e.g. parameters); rows of
0021 %       Z are applied to X1, X2, etc (e.g. population names). See examples
0022 %       in dsVary2Modifications.
0023 %
0024 % Outputs:
0025 %   - TODO
0026 %
0027 % Examples:
0028 %   - modifying population size and parameters:
0029 %       modifications={'E','size',5; 'E','gNa',120};
0030 %       model=dsApplyModifications(model,modifications);
0031 %
0032 %   - modifying mechanism_list:
0033 %       m=dsApplyModifications('dv/dt=10+@current; {iNa,iK}',...
0034 %                            {'pop1','mechanism_list','-iNa'});
0035 %       m.populations.mechanism_list
0036 %       m=dsApplyModifications('dv/dt=10+@current; {iNa,iK}',...
0037 %                            {'pop1','mechanism_list','+iCa'});
0038 %       m.populations.mechanism_list
0039 %       m=dsApplyModifications('dv/dt=10+@current; {iNa,iK}',...
0040 %                            {'pop1','mechanism_list','+(iCa,iCan,CaBuffer)'});
0041 %       m.populations.mechanism_list
0042 %
0043 %   - modifying equations (using special "cat()" operator or direct substitution)
0044 %       m=dsApplyModifications('dv/dt=10+@current; {iNa,iK}',...
0045 %                            {'pop1','equations','cat(dv/dt,+I)'});
0046 %       m.populations.equations
0047 %       m=dsApplyModifications('dv/dt=I(t)+@current; I(t)=10; {iNa,iK}',...
0048 %                            {'pop1','equations','cat(I(t),+sin(2*pi*t))'});
0049 %       m.populations.equations
0050 %       m=dsApplyModifications('dv/dt=I(t)+@current; I(t)=10; {iNa,iK}',...
0051 %                            {'pop1','equations','dv/dt=10+@current'});
0052 %       m.populations.equations
0053 %       m.populations.mechanism_list
0054 %
0055 %   - modifying equations with reserved keywords "ODEn" and "FUNCTIONn"
0056 %     - Note:
0057 %       'ODEn' = reserved keyword referencing the n-th ODE in equations
0058 %       'ODE' = aliases ODE1
0059 %       similarly: 'FUNCTIONn' and 'FUNCTION'
0060 %
0061 %       m=dsApplyModifications('dv/dt=10+@current; {iNa,iK}',...
0062 %                            {'pop1','equations','cat(ODE,+I)'});
0063 %       m.populations.equations
0064 %       m=dsApplyModifications('dv/dt=10+@current; du/dt=-u; {iNa,iK}',...
0065 %                            {'pop1','equations','cat(ODE2,+I)'});
0066 %       m.populations.equations
0067 %       m=dsApplyModifications('dv/dt=I(t)+@current; I(t)=10; {iNa,iK}',...
0068 %                            {'pop1','equations','cat(FUNCTION,+sin(2*pi*t))'});
0069 %       m.populations.equations
0070 %
0071 % See also: dsGenerateModel, dsSimulate, dsVary2Modifications
0072 %
0073 % Author: Jason Sherfey, PhD <jssherfey@gmail.com>
0074 % Copyright (C) 2016 Jason Sherfey, Boston University, USA
0075 
0076 %% localfn output
0077 if ~nargin
0078   output = localfunctions; % output var name specific to this fn
0079   return
0080 end
0081 
0082 %% auto_gen_test_data_flag argin
0083 options = dsCheckOptions(varargin,{'auto_gen_test_data_flag',0,{0,1}},false);
0084 if options.auto_gen_test_data_flag
0085   varargs = varargin;
0086   varargs{find(strcmp(varargs, 'auto_gen_test_data_flag'))+1} = 0;
0087   varargs(end+1:end+2) = {'unit_test_flag',1};
0088   argin = [{model},{modifications}, varargs]; % specific to this function
0089 end
0090 
0091 %%
0092 % check for modifications
0093 if isempty(modifications)
0094   % nothing to do
0095   output = model;
0096   return;
0097 end
0098 
0099 % check input type
0100 if ~isfield(model,'state_variables')
0101   ismodel = 0;
0102 else
0103   ismodel = 1;
0104 end
0105 
0106 % check specification
0107 if ismodel
0108   specification=dsCheckSpecification(model.specification, varargin{:});
0109 else
0110   specification=dsCheckSpecification(model, varargin{:});
0111 end
0112 
0113 % update specification with whatever is in modifications
0114 modifications = standardize_modifications(modifications,specification,varargin{:});
0115 % if ismodel % TODO: test this
0116   specification = modify_specification(specification,modifications,varargin{:});
0117 % end
0118 
0119 % update model if input was a model structure
0120 if ismodel
0121   output=dsGenerateModel(specification);
0122 else
0123   output = specification;
0124 end
0125 
0126 
0127 %% auto_gen_test_data_flag argout
0128 if options.auto_gen_test_data_flag
0129   argout = {output, modifications}; % specific to this function
0130 
0131   dsUnitSaveAutoGenTestData(argin, argout);
0132 end
0133 
0134 end
0135 
0136 function modifications = standardize_modifications(modifications,specification, varargin)
0137 % convert all modifications into 3-column cell matrix format
0138 % (namespace,variable,value)
0139 
0140 % 1. modifications structure
0141 % 2. partial specification structure
0142 
0143 %% auto_gen_test_data_flag argin
0144 options = dsCheckOptions(varargin,{'auto_gen_test_data_flag',0,{0,1}},false);
0145 if options.auto_gen_test_data_flag
0146   varargs = varargin;
0147   varargs{find(strcmp(varargs, 'auto_gen_test_data_flag'))+1} = 0;
0148   varargs(end+1:end+2) = {'unit_test_flag',1};
0149   argin = [{modifications},{specification}, varargs]; % specific to this function
0150 end
0151 
0152 if isstruct(modifications)
0153   % TODO
0154   % convert structure to cell matrix
0155   % ...
0156 end
0157 
0158 % check backward compatibility
0159 modifications=backward_compatibility(modifications);
0160 
0161 % check for empty object specification
0162 missing_objects=find(cellfun(@isempty,modifications(:,1)));
0163 if any(missing_objects)
0164   % set to default population name
0165   for i=1:length(missing_objects)
0166     modifications{missing_objects(i),1}=specification.populations(1).name;%'pop1';
0167   end
0168 end
0169 
0170 % support modifying multiple elements (of namespace or variable) simultaneously
0171 % approach -- add extra entry for each thing to modify
0172 % eg) expand {'(E,I)','Cm',2} to {'E','Cm',2; 'I','Cm',2}
0173 % note: should be able to support {'(E,I)','(EK,EK2)',-80}
0174 if any(~cellfun(@isempty,regexp(modifications(:,1),'^\(.*\)$'))) || ...
0175    any(~cellfun(@isempty,regexp(modifications(:,2),'^\(.*\)$')))
0176 
0177   % loop over modifications
0178   modifications_={};
0179   for i=1:size(modifications,1)
0180     % check namespace for ()
0181     namespaces=regexp(modifications{i,1},'[\w\.\-<>]+','match');
0182 
0183     % check variable for ()
0184     variables=regexp(modifications{i,2},'[\w\.-]+','match');
0185 
0186     if ischar(modifications{i,3})
0187 
0188         % expand list of modifications
0189         for j=1:length(namespaces)
0190             for k=1:length(variables)
0191                 modifications_(end+1,1:3)={namespaces{j},variables{k},modifications{i,3}};
0192             end
0193         end
0194 
0195     elseif isnumeric(modifications{i,3})
0196 
0197         % check size of values matches number of namespaces, variables
0198         if isscalar(modifications{i,3}) % in case number of values is one
0199             modifications{i,3} = repmat(modifications{i,3},length(variables),length(namespaces));
0200         else
0201             if size(modifications{i,3},1) ~= length(variables) || size(modifications{i,3},2) ~= length(namespaces)
0202                 % in case values is number of variables x 1
0203                 if size(modifications{i,3},1) == length(variables) && size(modifications{i,3},2) == 1
0204                     modifications{i,3} = repmat(modifications{i,3},1,length(namespaces));
0205                     % in case values is 1 x number of namespaces
0206                 elseif size(modifications{i,3},2) == length(namespaces) && size(modifications{i,3},1) == 1
0207                     modifications{i,3} = repmat(modifications{i,3},length(variables),1);
0208                     % TODO: char inputs
0209                     % elseif ischar(modifications{i,3})
0210                     % string input
0211                 else % if ~ischar(modifications{i,3})
0212                     error(['Numerical values varied over must be in array format,',...
0213                         'where dimensions 1, 2, and 3 correspond to mechanisms, values, and populations varied over.'])
0214                 end
0215             end
0216         end
0217 
0218         % expand list of modifications
0219         for j=1:length(namespaces)
0220             for k=1:length(variables)
0221                 modifications_(end+1,1:3)={namespaces{j},variables{k},modifications{i,3}(k,j)};
0222             end
0223         end
0224 
0225     end
0226 
0227   end
0228 
0229   modifications=modifications_;
0230 
0231 end
0232 
0233 %% auto_gen_test_data_flag argout
0234 if options.auto_gen_test_data_flag
0235   argout = {modifications}; % specific to this function
0236 
0237   dsUnitSaveAutoGenTestDataLocalFn(argin, argout); % localfn
0238 end
0239 
0240 end
0241 
0242 
0243 function spec = modify_specification(spec,mods, varargin)
0244 %% auto_gen_test_data_flag argin
0245 options = dsCheckOptions(varargin,{'auto_gen_test_data_flag',0,{0,1}},false);
0246 if options.auto_gen_test_data_flag
0247   varargs = varargin;
0248   varargs{find(strcmp(varargs, 'auto_gen_test_data_flag'))+1} = 0;
0249   varargs(end+1:end+2) = {'unit_test_flag',1};
0250   argin = [{spec},{mods}, varargs]; % specific to this function
0251 end
0252 
0253 precision=8; % number of digits allowed for user-supplied values
0254 predefined_variables={'name','size','parameters','mechanism_list','equations'};
0255 pop_names={spec.populations.name};
0256 if ~isempty(spec.connections)
0257   con_names=arrayfun(@(x)[x.source '->' x.target],spec.connections,'uni',0);
0258 else
0259   con_names=[];
0260 end
0261 
0262 % loop over modifications to apply
0263 for i=1:size(mods,1)
0264   obj=mods{i,1}; % population name or connection source-target
0265   fld=mods{i,2}; % population name, size, or parameter name
0266 
0267   if ~ischar(obj)
0268     error('modification must be applied to population name or connection source-target');
0269   end
0270 
0271   if ~ischar(fld) %|| ~ismember(fld,{'name','size','parameters','mechanism_list','equations'})
0272     error('modification must be applied to population ''name'',''size'',or a parameter referenced by its name');
0273   end
0274 
0275   % standardize connection object: convert target<-source to source->target
0276   if any(strfind(obj,'<-'))
0277     ind=strfind(obj,'<-');
0278     obj=[obj(ind(1)+2:end) '->' obj(1:ind(1)-1)];
0279   end
0280 
0281   % check and adjust for mechanism-specific parameter identifier
0282 %   MECH='';
0283   if any(obj=='.') % OBJECT.MECH
0284     tmp=regexp(obj,'\.','split');
0285     obj=tmp{1};
0286     MECH=tmp{2};
0287     fld=[MECH '.' fld];
0288   end
0289 %   if any(fld=='.') % MECH.PARAM
0290 %     tmp=regexp(fld,'\.','split');
0291 %     MECH=tmp{1};
0292 %     fld=tmp{2};
0293 %   end
0294 %   if ~isempty(MECH)
0295 %     fld=[MECH '.' fld];
0296 %   end
0297 
0298   val=mods{i,3}; % value for population name or size, or parameter to modify
0299   if ismember(obj,pop_names)
0300     type='populations';
0301     names=pop_names;
0302   elseif ismember(obj,con_names)
0303     type='connections';
0304     names=con_names;
0305   else
0306     warning('name of object to modify not found in populations or connections.');
0307     continue
0308   end
0309 
0310   index=ismember(names,obj);
0311 
0312   if strcmp(fld,'mechanism_list')
0313     % support --
0314     % 'E'    'mechanism_list'    '(iNa,iK)'
0315     % 'E'    'mechanism_list'    '-(iNa,iK)'
0316     % 'E'    'mechanism_list'    '+(iK)'
0317     % 'E'    'mechanism_list'    '+iK'
0318     % 'E'    'mechanism_list'    'iK'
0319 
0320     elems=regexp(val,'[\w@]+','match');
0321 
0322     if strcmp(val(1),'+')
0323       % add mechanisms to existing list
0324       spec.(type)(index).mechanism_list=unique_wrapper(cat(2,spec.(type)(index).mechanism_list,elems),'stable');
0325     elseif strcmp(val(1),'-')
0326       % remove mechanisms from existing list
0327       spec.(type)(index).mechanism_list=setdiff(spec.(type)(index).mechanism_list,elems,'stable');
0328     else
0329       % redefine mechanism list
0330       spec.(type)(index).mechanism_list=elems;
0331     end
0332   elseif strcmp(fld,'equations') && ~isempty(regexp(val,'^cat(','once'))
0333     % process special case: cat(TARGET,EXPRESSION)
0334     eqns=spec.(type)(index).equations;
0335     args=regexp(val,'cat\((.+)\)','tokens','once');
0336     ind=find(args{1}==',',1,'first');
0337     target=args{1}(1:ind-1);
0338     expression=args{1}(ind+1:end);
0339 %     args=regexp(args{1},',','split');
0340 %     target=args{1};
0341 %     expression=args{2};
0342     % support keywords: ODEn, ODE=ODE1, FUNCTIONn, FUNCTION=FUNCTION1
0343     if strncmp('ODE',target,3)
0344       % support target = ODEn and ODE=ODE1
0345       % replace target by n-th ODE LHS
0346       lines=cellfun(@strtrim,strsplit(eqns,';'),'uni',0); % split equations into statements
0347       lines=lines(~cellfun(@isempty,lines)); % eliminate empty elements
0348       pattern='^((\w+'')|(d\w+/dt))\s*='; % pattern for ODEs
0349       inds=regexp(lines,pattern,'once'); % indices to ODE statements
0350       inds=find(~cellfun(@isempty,inds));
0351       % get index to the ODE statement to modify
0352       if strcmp(target,'ODE')
0353         ind=inds(1);
0354       else
0355         n=str2num(target(4:end));
0356         ind=inds(n);
0357       end
0358       % get LHS of ODE
0359       %LHS=regexp(lines{ind},'^(.+)=','tokens','once');
0360       LHS=regexp(lines{ind},'^([^=]+)=','tokens','once');
0361       target=LHS{1};
0362     elseif strncmp('FUNCTION',target,8)
0363       % support target = FUNCTIONn and FUNCTION=FUNCTION1
0364       % replace target by n-th FUNCTION LHS
0365       lines=cellfun(@strtrim,strsplit(eqns,';'),'uni',0); % split equations into statements
0366       lines=lines(~cellfun(@isempty,lines)); % eliminate empty elements
0367       pattern='^\w+\([a-zA-Z][\w,]*\)\s*='; % pattern for functions
0368       inds=regexp(lines,pattern,'once'); % indices to function statements
0369       inds=find(~cellfun(@isempty,inds));
0370 
0371       % get index to the function statement to modify
0372       if strcmp(target,'FUNCTION')
0373         ind=inds(1);
0374       else
0375         n=str2num(target(9:end));
0376         ind=inds(n);
0377       end
0378       % get LHS of ODE
0379       %LHS=regexp(lines{ind},'^(.+)=','tokens','once');
0380       LHS=regexp(lines{ind},'^([^=]+)=','tokens','once');
0381       target=LHS{1};
0382     end
0383     % add escape character for using regular expression to match function statements
0384     target=strrep(target,'(','\(');
0385     target=strrep(target,')','\)');
0386 
0387     % modify equations
0388     old=regexp(eqns,[target '\s*=[^;]+'],'match');
0389     if ~isempty(old)
0390       eqns=strrep(eqns,old{1},[old{1} expression]);
0391       spec.(type)(index).equations=eqns;
0392     end
0393   elseif ismember(fld,predefined_variables) % not a single parameter to modify
0394     spec.(type)(index).(fld)=val;
0395     if strcmp(type,'populations') && strcmp(fld,'name')
0396       % update name in connections
0397       for j=1:length(spec.connections)
0398         if strcmp(pop_names{index},spec.connections(j).source)
0399           spec.connections(j).source=val;
0400         end
0401 
0402         if strcmp(pop_names{index},spec.connections(j).target)
0403           spec.connections(j).target=val;
0404         end
0405       end
0406     end
0407   else % modify a single parameter in the populations.parameters cell array
0408     param_names=spec.(type)(index).parameters(1:2:end);
0409     if isempty(spec.(type)(index).parameters)
0410       % no parameters set; start cell array of parameters
0411       spec.(type)(index).parameters={fld,val};%toString(val2,precision)};
0412     elseif ismember(fld,param_names)
0413       % parameter found in list; update its value
0414       val_pos=2*find(ismember(param_names,fld));
0415       spec.(type)(index).parameters{val_pos}=val;%toString(val2,precision);
0416     else
0417       % parameter not in existing list; append to end
0418       spec.(type)(index).parameters{end+1}=fld;
0419       spec.(type)(index).parameters{end+1}=val;%toString(val2,precision);
0420     end
0421   end
0422 end
0423 
0424 % todo: split X.MECH, set (type,index) from X, store {MECH.fld,val} in .parameters
0425 
0426 %% auto_gen_test_data_flag argout
0427 if options.auto_gen_test_data_flag
0428   argout = {spec}; % specific to this function
0429 
0430   dsUnitSaveAutoGenTestDataLocalFn(argin, argout); % localfn
0431 end
0432 
0433 end
0434 
0435 
0436 function modifications=backward_compatibility(modifications)
0437 % convert 2-column specification to 3-column specification with empty object name
0438 if size(modifications,2)==2
0439   tmp={};
0440   for i=1:size(modifications,1)
0441     tmp{i,1}='';
0442     tmp{i,2}=modifications{i,1};
0443     tmp{i,3}=modifications{i,2};
0444   end
0445   modifications=tmp;
0446 end
0447 
0448 % convert 4-column specification to 3-column
0449 if size(modifications,2)==4
0450   for i=1:size(modifications,1)
0451     if strcmp(modifications{i,2},'parameters')
0452       % shift parameter name to 2nd column
0453       modifications{i,2}=modifications{i,3};
0454 
0455       % shift parameter value to 3rd column
0456       modifications{i,3}=modifications{i,4};
0457     end
0458   end
0459 
0460   % remove fourth column
0461   modifications=modifications(:,1:3);
0462 end
0463 
0464 % convert connection reference source-target to source->target
0465 if any(~cellfun(@isempty,regexp(modifications(:,1),'\w-\w')))
0466   % find elements to adjust
0467   inds=find(~cellfun(@isempty,regexp(modifications(:,1),'\w-\w')));
0468   for i=1:length(inds)
0469     modifications{inds(i),1}=strrep(modifications{inds(i),1},'-','->');
0470   end
0471 end
0472 
0473 end
dsApplyModifications

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
