From 7526e4f9c1a49a95621ee37e8d68c59b3f414e67 Mon Sep 17 00:00:00 2001 From: Niklas Neubrand <38315848+niklasneubrand@users.noreply.github.com> Date: Fri, 13 Sep 2024 22:37:47 +0200 Subject: [PATCH] consider cyclic replacements in SBML Export Update to improve compatibility of SBML as a export format: * add check if CONDITIONS can be mapped to initAssignments in SBML * add fnction to help fix dependent conditions * improve encoding of conditions and initial values in SBML --- .../ImportExport/arCheckSBMLCompatibilty.m | 83 ++++++ .../ImportExport/arExportSBML_FullModel.m | 242 ++++++++-------- arFramework3/ImportExport/arParseSBML.m | 20 +- .../ImportExport/arRenameModelCondPars.m | 261 ++++++++++++++++++ 4 files changed, 472 insertions(+), 134 deletions(-) create mode 100644 arFramework3/ImportExport/arCheckSBMLCompatibilty.m create mode 100644 arFramework3/ImportExport/arRenameModelCondPars.m diff --git a/arFramework3/ImportExport/arCheckSBMLCompatibilty.m b/arFramework3/ImportExport/arCheckSBMLCompatibilty.m new file mode 100644 index 00000000..9d6f43db --- /dev/null +++ b/arFramework3/ImportExport/arCheckSBMLCompatibilty.m @@ -0,0 +1,83 @@ +function [qSBMLCompatible, pProblemCond] = arCheckSBMLCompatibilty(m) +% ARCHECHECKSBMLCOMPATIBILTY checks if CONDITIONS in model.def are SBML compatible +% +% BACKGROUND: +% CONDITIONS in model.def correspond to replacements of the original model parameters +% with numerical values of mathematical expressions. In SBML this is best represented as +% "initialAssignment" rules. However, there is an important difference how CONDITIONS +% and initialAssignment rules are applied: +% - CONDITIONS are applied exactly once and independently of each other +% - initialAssignment rules are combined to form a set of equations that are solved +% simultaneously +% +% This means that conditions should be encoded as initialAssignment rules*, but cannot +% be translated directly in all cases. Specifically, there must not be any dependencies +% between the conditions. This means, a parameter that is repalced must not appear in +% any mathematical expression on the RHS of the replacements. +% +% * See arExportSBML_fullmodel.m +% +% EXAMPLE OF DIFFERENT BEHAVIOR: +% consider the expression expr = "p1*x^p2" +% +% d2d: +% CONDITIONS +% p1 "2" +% p2 "p1*2" +% +% The expression after replacements becomes: "2*x^(2*p1)" +% +% SBML: +% +% +% +% 2 +% +% +% +% +% +% +% p1 +% 2 +% +% +% +% +% +% The expression after replacements becomes: "2*x^(4)" +% This is because the assignment p1<-2 is also applied within the assignment p2<-p1*2 + +arguments + m (1,1) double {mustBeInteger, mustBePositive} = 1 +end + + +global ar + +% identify parameters that are replaced in the conditions +modelP = ar.model(m).p'; +modelFPsym = arSym(ar.model(m).fp); +modelFP = cellstr(string(arSym(ar.model(m).fp))); +qReplaced = ~strcmp(modelP, modelFP); + +% identify which of the replaced parameters also appear in modelFP expressions +fpParams = unique(cellstr(string(symvar(modelFPsym(qReplaced))))); +qReplacedAndInFP = ismember(modelP, fpParams) & qReplaced; + +% exclude init parameters that do not appear in the model equations +qStateInit = ismember(modelP, ar.model(m).px0); +qInModelEquations = ismember(modelP, union(ar.model(m).pu, ar.model(m).pv)); +qPureStateInit = qStateInit & ~qInModelEquations; + +% exclude compartment size parameters +% Compartment size parameters do not pose a problem, because they are not used in SBML. +% Instead compartment valume is assigned directly to the compartment ID in the SBML file. +qCompSizeParam = ismember(ar.model(m).p', ar.model(m).pc); + +% final array of parameters with problematic replacements +qProblemCond = qReplacedAndInFP & ~qCompSizeParam & ~qPureStateInit; +pProblemCond = modelP(qProblemCond); + +% Flag that indicates if the model is SBML compatible +qSBMLCompatible = ~any(qProblemCond); \ No newline at end of file diff --git a/arFramework3/ImportExport/arExportSBML_FullModel.m b/arFramework3/ImportExport/arExportSBML_FullModel.m index 134d8403..1f9d218d 100644 --- a/arFramework3/ImportExport/arExportSBML_FullModel.m +++ b/arFramework3/ImportExport/arExportSBML_FullModel.m @@ -11,7 +11,11 @@ function arExportSBML_FullModel(m,name) global ar % simulate once for initial values - arSimu(0,1,0) + try + arSimu(0,1,0) + catch + arSimu(1,1,0) + end copasi = true; @@ -26,12 +30,18 @@ function arExportSBML_FullModel(m,name) % steadystate = false; % end % end + + qCondsSBMLConform = arCheckSBMLCompatibilty(m); + if ~qCondsSBMLConform + warning(['Model conditions are not independet and will be represented incorrectly in SBML. ', ... + 'Consider using "arRenameModelCondPars" to get independent model parameters and conditions.']); + end try M = TranslateSBML(which('empty.xml')); F = TranslateSBML(which('filled.xml')); catch - warning('error in libSBML. Probably backwards compatibility issues with old MATLAB version. Should work with 2019a') + warning('error in libSBML. Probably backwards compatibility issues with old MATLAB version. Should work with 2019a.') end M.id = ar.model(m).name; @@ -180,17 +190,6 @@ function arExportSBML_FullModel(m,name) M.initialAssignment(idxIA).annotation = ''; M.initialAssignment(idxIA).sboTerm = -1; M.initialAssignment(idxIA).symbol = M.compartment(jc).id; - - % Translate the compartment size parameter to a MathML expression - % if isempty(str2num(cSizeReplace)) - % mathType = 'ci'; - % else - % mathType = 'cn'; - % end - % mathExpr = sprintf('<%s>%s', mathType, cSizeReplace, mathType); - % M.initialAssignment(idxIA).math = ... - % sprintf('%s', mathExpr); - M.initialAssignment(idxIA).math = char(cSizeReplace); M.initialAssignment(idxIA).level = 2; M.initialAssignment(idxIA).version = 4; @@ -264,9 +263,7 @@ function arExportSBML_FullModel(m,name) M.species(jx).isSetCharge = 0; M.species(jx).level = 2; M.species(jx).version = 4; - - qp = ismember(ar.pLabel, ar.model(m).px0{jx}); %R2013a compatible - + simulated_ss = 0; % if ( steadystate ) % if (isfield(ar.model(m).condition(c),'ssLink') && ~isempty( ar.model(m).condition(c).ssLink ) ) @@ -275,10 +272,10 @@ function arExportSBML_FullModel(m,name) % simulated_ss = 1; % end % end - qp = ismember(ar.pLabel, ar.model(m).px0{jx}); %R2013a compatible if ( ~simulated_ss ) % check if init parameter still exists in condition parameters + qp = ismember(ar.pLabel, ar.model(m).px0{jx}); %R2013a compatible is_set = sum(ismember(ar.model(m).fp, ar.model(m).px0{jx}))==0; if(sum(qp)==1 && is_set==0) M.species(jx).initialConcentration = 1; @@ -380,121 +377,118 @@ function arExportSBML_FullModel(m,name) function [M] = GetParameters(M,m) global ar - - - %% first: check the model parameters in ar.model(m).p - % allPars = model parameters that appear in at least one condition - % or one input. Compartment izes are excluded. - - allPars = zeros(1, length(ar.model(m).p)); - constPars = ones(1, length(ar.model(m).p)); - for condId = 1:length(ar.model(m).condition) - condPars = ismember(ar.model(m).p,ar.model(m).condition(condId).p)'; - allPars = condPars | allPars; - end - % loop necessary? index unused - for inpId = 1:length(ar.model(m).u) - inpPars = ismember(ar.model(m).p, ar.model(m).pu); - allPars = inpPars | allPars; - end - % if parameter is an initial condition + + %% first: collect available numerical values for model parameters (ar.model.p) + + % logical flags for model parameters isInit = cellfun(@(x) any(strcmp(x, ar.model(m).px0)), ar.model(m).p)'; - % if model parameter is replaced by other parameter - isReplaced = zeros(1, length(ar.model(m).p)); - for jp=1:length(ar.model(m).p) - if isempty(regexp(ar.model(m).p{jp}, ar.model(m).fp{jp}, "once")) - isReplaced(jp) = true; - end - end - isReplaced = logical(isReplaced); - % if parameter is compartment size - isCompSize = cellfun(@(x) any(strcmp(x, ar.model(m).pc)), ar.model(m).p)'; - allPars = allPars & ~isCompSize; % exclude compartment sizes - + isReplaced = ~strcmp(ar.model.p', string(arSym(ar.model.fp))); + isCompSize = cellfun(@(x) any(strcmp(x, string(arSym(ar.model(m).pc)))), ar.model(m).p)'; + % isConst = (isReplaced & ~isInit) | (isInit & ~isReplaced); + isConst = true(1, length(ar.model(m).p)); % all model parameters are constant (i.e. not time-dependent) + for jp = 1:length(ar.model(m).p) - if isReplaced(jp) && ~isInit(jp) && ~isCompSize(jp) - - allPars(jp) = true; - constPars(jp) = false; + % All model parameters should be defined in SBML as parameters + % irrespective of initAssigns or undefined values. + + % compound parameters are handeled separately in "GetCompartments" + if isCompSize(jp) + continue + end + + % Is there a numeric value for the model parameter in ar.p? + qp = strcmp(ar.model(m).p(jp), ar.pLabel); %R2013a compatible + if any(qp) + % get parameter value from ar.p + pvalue = ar.p(qp); + if(ar.qLog10(qp) == 1) + pvalue = 10^pvalue; + end + isSetValue = 1; + constant = double(isConst(qp)); + else + % no numeric value in ar.p found + pvalue = NaN; + isSetValue = 0; + constant = 1; + end + + id_tmp = length(M.parameter) + 1; + M.parameter(id_tmp).typecode = 'SBML_PARAMETER'; + M.parameter(id_tmp).metaid = ''; + M.parameter(id_tmp).notes = ''; + M.parameter(id_tmp).annotation = ''; + M.parameter(id_tmp).sboTerm = -1; + M.parameter(id_tmp).name = ar.model(m).p{jp}; + M.parameter(id_tmp).id = ar.model(m).p{jp}; + M.parameter(id_tmp).units = ''; + M.parameter(id_tmp).constant = constant; + M.parameter(id_tmp).isSetValue = isSetValue; + M.parameter(id_tmp).value = pvalue; + M.parameter(id_tmp).level = 2; + M.parameter(id_tmp).version = 4; - fu = arSym(ar.model(m).fp{jp}); - % replace time parameters with 'time' - fu = char(arSubs(arSym(fu), arSym(ar.model(m).t), arSym('time'))); - - ixrule = length(M.rule) + 1;% index of current rule - M.rule(ixrule).typecode = 'SBML_ASSIGNMENT_RULE'; - M.rule(ixrule).metaid = ''; - M.rule(ixrule).notes = ''; - M.rule(ixrule).annotation = ''; - M.rule(ixrule).sboTerm = -1; - M.rule(ixrule).formula = fu; - M.rule(ixrule).variable = ar.model(m).p{jp}; - M.rule(ixrule).species = ''; - M.rule(ixrule).compartment = ''; - M.rule(ixrule).name = ''; - M.rule(ixrule).units = ''; - M.rule(ixrule).level = 2; - M.rule(ixrule).version = 4; - - elseif ~isReplaced(jp) && isInit(jp) - - allPars(jp) = true; - constPars(jp) = false; - - end end - - % add parameters to model - for id = 1:length(allPars) - if allPars(id) == 1 - - % Is there a numeric value for the model parameter in ar.p? - qp = strcmp(ar.model(m).p(id), ar.pLabel); %R2013a compatible - if any(qp) - % if possible: get parameter value from ar.p - pvalue = ar.p(qp); - if(ar.qLog10(qp) == 1) - pvalue = 10^pvalue; + + + %% second: collect replacements of model parameters + + for jp = 1:length(ar.model(m).p) + if isReplaced(jp) && ~isCompSize(jp) + + if isInit(jp) + % initial values for species are already defined in "GetSpecies" + % and written to SBML in "GetInitialAssignments". + % However: An init parameter can also appears as a model parameter + % (e.g. in a rate equation or input). Then the CONDITION + % init_State "expression" + % must also be applied to the parameter, not just the species. + if ~ismember(ar.model(m).p{jp}, union(ar.model(m).pu, ar.model(m).pv)) + continue end - else - continue end + + % parameter CONDITIONS should be implemeted as initialAssignment + % reason: parameters in d2d are constant, CONDITIONS are applied before start of simulation - id_tmp = length(M.parameter) + 1; - M.parameter(id_tmp).typecode = 'SBML_PARAMETER'; - M.parameter(id_tmp).metaid = ''; - M.parameter(id_tmp).notes = ''; - M.parameter(id_tmp).annotation = ''; - M.parameter(id_tmp).sboTerm = -1; - M.parameter(id_tmp).name = ar.model(m).p{id}; - M.parameter(id_tmp).id = ar.model(m).p{id}; - M.parameter(id_tmp).units = ''; - M.parameter(id_tmp).constant = constPars(id); - M.parameter(id_tmp).isSetValue = 1; - M.parameter(id_tmp).value = pvalue; - M.parameter(id_tmp).level = 2; - M.parameter(id_tmp).version = 4; - end + ixInitAssign = length(M.initialAssignment) + 1;% index of current rule + M.initialAssignment(ixInitAssign).typecode = 'SBML_INITIAL_ASSIGNMENT'; + M.initialAssignment(ixInitAssign).metaid = ''; + M.initialAssignment(ixInitAssign).notes = ''; + M.initialAssignment(ixInitAssign).annotation = ''; + M.initialAssignment(ixInitAssign).sboTerm = -1; + M.initialAssignment(ixInitAssign).symbol = ar.model(m).p{jp}; + M.initialAssignment(ixInitAssign).math = ar.model(m).fp{jp}; + M.initialAssignment(ixInitAssign).level = 2; + M.initialAssignment(ixInitAssign).version = 4; + + % ixrule = length(M.rule) + 1;% index of current rule + % M.rule(ixrule).typecode = 'SBML_ASSIGNMENT_RULE'; + % M.rule(ixrule).metaid = ''; + % M.rule(ixrule).notes = ''; + % M.rule(ixrule).annotation = ''; + % M.rule(ixrule).sboTerm = -1; + % M.rule(ixrule).formula = ar.model(m).fp{jp}; + % M.rule(ixrule).variable = ar.model(m).p{jp}; + % M.rule(ixrule).species = ''; + % M.rule(ixrule).compartment = ''; + % M.rule(ixrule).name = ''; + % M.rule(ixrule).units = ''; + % M.rule(ixrule).level = 2; + % M.rule(ixrule).version = 4; + end end - %% second: add parameters from replacements + %% third: search numerical values for replacements - additionalPars = zeros(1, length(ar.pLabel)); - for condId = 1:length(ar.model(m).condition) - condPars = ismember(ar.pLabel,ar.model(m).condition(condId).p) ; - additionalPars = condPars|additionalPars; - end - % loop necessary? index unused - for inpId = 1:length(ar.model(m).u) - inpPars = ismember(ar.pLabel',ar.model(m).pu); - additionalPars = inpPars|additionalPars; - end - additionalPars = additionalPars & ~ismember(ar.pLabel,ar.model(m).p); + % optimization parameters that appear in replacements (and are not model parameters) + replParams = cellfun(@(x) any(contains(ar.model(m).fp(isReplaced), x)), ar.pLabel); + replParams = replParams & ~cellfun(@(x) ismember(x, ar.model(m).p), ar.pLabel); - - for id = 1:length(additionalPars) - if additionalPars(id) == 1 + for jp = 1:length(ar.pLabel) + + if replParams(jp) id_tmp = length(M.parameter) + 1; M.parameter(id_tmp).typecode = 'SBML_PARAMETER'; @@ -502,16 +496,16 @@ function arExportSBML_FullModel(m,name) M.parameter(id_tmp).notes = ''; M.parameter(id_tmp).annotation = ''; M.parameter(id_tmp).sboTerm = -1; - M.parameter(id_tmp).name = ar.pLabel{id}; - M.parameter(id_tmp).id = ar.pLabel{id}; + M.parameter(id_tmp).name = ar.pLabel{jp}; + M.parameter(id_tmp).id = ar.pLabel{jp}; M.parameter(id_tmp).units = ''; M.parameter(id_tmp).constant = 1; M.parameter(id_tmp).isSetValue = 1; M.parameter(id_tmp).level = 2; M.parameter(id_tmp).version = 4; - pvalue = ar.p(id); - if(ar.qLog10(id) == 1) + pvalue = ar.p(jp); + if(ar.qLog10(jp) == 1) pvalue = 10^pvalue; end M.parameter(id_tmp).value = pvalue; diff --git a/arFramework3/ImportExport/arParseSBML.m b/arFramework3/ImportExport/arParseSBML.m index 947337ee..9c5e2b65 100644 --- a/arFramework3/ImportExport/arParseSBML.m +++ b/arFramework3/ImportExport/arParseSBML.m @@ -685,22 +685,22 @@ elseif initValuesOpt(j) == 2 assignment_value = replacePowerFunction(m.initialAssignment(idx_assval(j)).math,0); assignment_value = arSym(assignment_value); - assignment_value = arSubs(assignment_value, arSym({m.species.name}), arSym(cellfun(@(x) ['init_',x], {m.species.name}, 'UniformOutput', false))); %replace state by init_state + assignment_value = arSubs(assignment_value, arSym({m.species.id}), arSym(cellfun(@(x) ['init_',x], {m.species.id}, 'UniformOutput', false))); %replace state by init_state - % assignment_value = arSubs(arSym(assignment_value), arSym(pars), arSym(par_value)); - %assignment_value = subs(assignment_value, specs, spec_value); - % assignment_value = arSubs(arSym(assignment_value), arSym(comps), arSym(comp_value)); - % assignment_value = eval(assignment_value); + % assignment_value = arSubs(arSym(assignment_value), arSym(pars), arSym(par_value)); + % assignment_value = subs(assignment_value, specs, spec_value); + % assignment_value = arSubs(arSym(assignment_value), arSym(comps), arSym(comp_value)); + % assignment_value = eval(assignment_value); - % ub = 1000; - % if(assignment_value>ub) - % ub = assignment_value * 10; - % end + % ub = 1000; + % if(assignment_value>ub) + % ub = assignment_value * 10; + % end fprintf(fid, '%s\t "%s"\n', ['init_' sym_check(m.time_symbol,m.species(j).id2)], ... assignment_value); end - + end fprintf(fid, '\nPARAMETERS\n'); diff --git a/arFramework3/ImportExport/arRenameModelCondPars.m b/arFramework3/ImportExport/arRenameModelCondPars.m new file mode 100644 index 00000000..b2bdc0b1 --- /dev/null +++ b/arFramework3/ImportExport/arRenameModelCondPars.m @@ -0,0 +1,261 @@ +function arRenameModelCondPars(m, targetDir, makeCopies) +% ARRENAMEMODELCONDITIONPARS removes dependencies between condition replacements +% by renaming the model parameters in the model.def (and if necessary data.def) file(s). +% +% Explanation: +% In d2d the replacements specified in CONDITIONS are applied exactly once and independently +% of each other. Let's illustrate this with the example of swapping two parameters in the model: +% myExpression = str2sym('x^y'); +% p = arSym({'x', 'y'}); +% fp = arSym({'y', 'x'}); +% arSubs(myExpression, p, fp) returns 'y^x' +% +% However, if we want to export the model as a SBML file, the conditions can only be encoded as +% "initialAssignment" rules. These rules are all combined to form a set of equations that are +% solved simultaneously. Therefore, the previous example cannot be encoded in SBML since there +% is a circular dependency: x depends on y and y depends on x. +% +% To resolve this, we want to rename the original parameters in the reactions, ODEs and inputs. +% In the example this would mean: +% myExpression = str2sym('x_model^y_model'); +% p = arSym({'x_model', 'y_model'}); +% fp = arSym({'y', 'x'}); +% arSubs(myExpression, p, fp) returns 'y^x' +% +% Finally, since there might also be conditions on the model parameters in the data.def files, +% we also need to rename the parameters in the data.def files. + +arguments + m (1,1) double + targetDir (1,1) string = '' + makeCopies (1,1) logical = true +end + + +global ar + + +% Check if model is SBML compatible +[qSBMLCompatible, pNeedsRenaming] = arCheckSBMLCompatibilty(m); +if qSBMLCompatible + fprintf('Model conditions already SBML compatible. No parameters need renaming.\n') + return +end + +%% Overwrite the project or create a new one +if isempty(targetDir) + % overwrite the project + targetDir = ar.info.path; +else + % copy the project to a new directory + targetDir = fullfile(targetDir); + if exist(targetDir, 'dir') + warning('Directory already exists: %s\n', targetDir) + end + copyfile(ar.info.path, targetDir) + makeCopies = false; +end + + +%% Read old and write new *.def files + +% find correct model.def file +callLoadModel = strcmp(ar.setup.commands, 'arLoadModel'); +if sum(callLoadModel) == 1 + modelDef = fullfile(ar.info.path, ar.setup.modelfiles{callLoadModel}); +elseif sum(callLoadModel) > 1 + % multiple models loades, search for the correct one + error('Multiple models loaded. Seraching the correct one is not implemented yet.') +else + error('No model.def file found in ar.setup.modelfiles.') +end + +% find all data.def files +callLoadData = strcmp(ar.setup.commands, 'arLoadData'); +dataFiles = ar.setup.datafiles(callLoadData); +dataDefs = cellfun(@(d) fullfile(ar.info.path, d{1}), dataFiles, 'UniformOutput', false); +dataDefs = unique(dataDefs); % remove duplicates (if multiple model.data share a data.def) +defFiles = [modelDef, dataDefs]; + +% flag for each parameter (did it appear in reactions, ODEs or inputs?) +qInModelReactions = false(1, length(pNeedsRenaming)); +qInDataReactions = false(length(dataDefs), length(pNeedsRenaming)); + +for fileIdx = 1:length(defFiles) + + didReplace = false; + defFile = defFiles{fileIdx}; + + % backup and target file + defFileTarget = strrep(defFile, ar.info.path, targetDir); + defFileCopy = sprintf('%s_original.def', defFile(1:end-4)); + defFileCopy = strrep(defFileCopy, ar.info.path, targetDir); + if makeCopies + copyfile(defFile, defFileCopy); + end + + % read file and split content into sections + fileContent = fileread(defFile); + defFileHeads = ["CONDITIONS", "PARAMETERS"]; + [sections, headings] = split(fileContent, defFileHeads); + + % modify file content + for ip = 1:length(pNeedsRenaming) + paramOld = pNeedsRenaming{ip}; + paramNew = [paramOld, '_model']; + + % modify the first section (before CONDITIONS) + % -> replace all occurrences of p by p_model + pattern1 = sprintf('\\<%s\\>', paramOld); + matches = regexp(sections{1}, pattern1, 'match'); + if ~isempty(matches) + didReplace = true; + if fileIdx == 1 + qInModelReactions(ip) = true; + else + qInDataReactions(fileIdx-1, ip) = true; + end + sections{1} = regexprep(sections{1}, pattern1, paramNew); + end + + % modify CONDITIONS section (line by line backwards) + % -> replace all occurrences of p by p_model in LHS of conditions + % -> also keep old conditions if p is a state init parameter + condSection = splitlines(sections{2}); + pattern2 = sprintf('^\\<%s\\>', paramOld); + for il = length(condSection):-1:1 + oldCond = condSection{il}; + matches = regexp(oldCond, pattern2, 'match'); + if ~isempty(matches) + didReplace = true; + if ~ismember(paramOld, ar.model(m).px0) + % normal parameter + % -> replace p by p_model in LHS of conditions + condSection{il} = regexprep(oldCond, pattern2, paramNew); + elseif qInModelReactions(ip) || qInDataReactions(fileIdx-1, ip) + % init param that also appears in reactions + % -> keep the old condition (for initial value of state) + % -> use p_model as new parameter in equations + condSection{il} = regexprep(oldCond, pattern2, paramNew); + condSection = [condSection(1:il); oldCond; condSection(il+1:end)]; + else + % init param that does not appear in reactions + % -> do nothing, no renaming necessary + end + end + end + sections(2) = join(condSection, newline); + end + + % console output + [~, fileName] = fileparts(defFile); + if didReplace + fprintf('%s.def: did replace model parameters.\n', fileName); + else + fprintf('%s.def: no replacements. No backup necessary.\n', fileName); + if makeCopies + delete(defFileCopy) + end + end + + % re-join the sections and write to the original file# + fileContent = join(sections, headings); + fid = fopen(defFileTarget, 'w'); + fprintf(fid, '%s', fileContent{:}); + fclose(fid); + +end + + +%% Read and modify data tables (if necessary) +dataTables = cellfun(@(d) fullfile(ar.info.path, d{2}), dataFiles, 'UniformOutput', false); +dataTables = unique(dataTables); % remove duplicates (if multiple model.data share a data table) + +for d = 1:length(dataTables) + + dataFile = dataTables{d}; + [~, dataName, ext] = fileparts(dataFile); + + % backup and target file + dataFileCopy = sprintf('%s_original%s', dataFile(1:end-4), ext); + dataFileCopy = strrep(dataFileCopy, ar.info.path, targetDir); + dataFileTarget = strrep(dataFile, ar.info.path, targetDir); + + % read, modify and write data file + didReplace = false; + switch ext + case {'.xls', '.xlsx'} + [~, ~, data] = xlsread(dataFile); + for ip = 1:length(pNeedsRenaming) + paramOld = pNeedsRenaming{ip}; + paramNew = [paramOld, '_model']; + qParam = strcmp(data(1,:), paramOld); + if any(qParam) + didReplace = true; + if ~ismember(paramOld, ar.model(m).px0) + % normal parameter + % -> replace p by p_model in column names + data{1, qParam} = paramNew; + elseif qInModelReactions(ip) || qInDataReactions(fileIdx-1, ip) + % init param that also appears in reactions + % -> keep the old column (for initial value of state) + % -> use p_model as new parameter in equations + idReplace = find(qParam); + dataInsert = data(:, idReplace); + dataInsert{1} = paramNew; + data = [data(:, 1:idReplace), dataInsert, data(:, idReplace+1:end)]; + else + % init param that does not appear in reactions + % -> do nothing, no renaming necessary + end + end + end + if didReplace + if makeCopies + copyfile(dataFile, dataFileCopy); + end + xlswrite(dataFileTarget, data); + end + case '.csv' + data = readtable(dataFile); + for ip = 1:length(pNeedsRenaming) + paramOld = pNeedsRenaming{ip}; + paramNew = [paramOld, '_model']; + varNames = data.Properties.VariableNames; + qParam = strcmp(varNames, paramOld); + if any(qParam) + didReplace = true; + if ~ismember(paramOld, ar.model(m).px0) + % normal parameter + % -> rename column + varNames{qParam} = paramNew; + data.Properties.VariableNames = varNames; + elseif qInModelReactions(ip) || qInDataReactions(fileIdx-1, ip) + % init param that also appears in reactions + % -> keep the old column (for initial value of state) + % -> use p_model as new parameter in equations + % -> add duplicate column with new name + idReplace = find(qParam); + dataInsert = data(:, idReplace); + dataInsert.Properties.VariableNames = {paramNew}; + data = [data(:, 1:idReplace), dataInsert, data(:, idReplace+1:end)]; + else + % init param that does not appear in reactions + % -> do nothing, no renaming necessary + end + end + end + if didReplace + if makeCopies + copyfile(dataFile, dataFileCopy); + end + writetable(data, dataFile); + end + end + if didReplace + fprintf('%s%s: did replace model parameters in column names.\n', dataName, ext); + else + fprintf('%s%s: no replacements. No backup necessary.\n', dataName, ext); + end + +end