Skip to content

Commit

Permalink
Fix log() issues #412 (#416)
Browse files Browse the repository at this point in the history
* Don't use log() with two arguments and confusing arguments order (Fixes #412)

* Disallow ambiguous 'log(x)' in sbml import, should use log(base, x) or ln(x), log10(x), ... explicitely (#412)

* Replace log(base, x) by ln(x)/ln(base) in observables (#412)

* Revert "Disallow ambiguous 'log(x)' in sbml import, should use log(base, x) or ln(x), log10(x), ... explicitely (#412)"

This reverts commit 635e48f. This was meaningless since l3ParserSettings parseLog option does not affect formulaToL3StringWithSettings but only parseL3FormulaWithSettings.

* Refactor(python) naming
  • Loading branch information
dweindl authored Aug 28, 2018
1 parent af50a9f commit c039ffd
Showing 1 changed file with 59 additions and 10 deletions.
69 changes: 59 additions & 10 deletions python/amici/sbml_import.py
Original file line number Diff line number Diff line change
Expand Up @@ -851,9 +851,15 @@ def computeModelEquationsObjectiveFunction(self, observables=None,

# add user-provided observables or make all species observable
if(observables):
# Replace logX(.) by log(X,.) since symengine cannot parse the former
# Replace logX(.) by log(., X) since symengine cannot parse the former
# also replace symengine-incompatible sbml log(basis, x)
for observable in observables:
observables[observable]['formula'] = re.sub(r'(^|\W)log(\d+)\(', r'\1log(\2, ', observables[observable]['formula'])
observables[observable]['formula'] = re.sub(r'(^|\W)log(\d+)\(', r'\g<1>1/ln(\2)*ln(', observables[observable]['formula'])
repl = replaceLogAB(observables[observable]['formula'])
if repl != observables[observable]['formula']:
print('Replaced "%s" by "%s", assuming first argument to log() was the basis.' % (observables[observable]['formula'], repl))
observables[observable]['formula'] = repl

self.observables = sp.DenseMatrix([observables[observable]['formula'] for observable in observables])
observableNames = [observables[observable]['name'] if 'name' in observables[observable].keys()
else 'y' + str(index)
Expand Down Expand Up @@ -1480,12 +1486,12 @@ class TemplateAmici(Template):
delimiter = 'TPL_'


def assignmentRules2observables(sbml, filter = lambda *_: True):
def assignmentRules2observables(sbml_model, filter_function=lambda *_: True):
"""Turn assignment rules into observables.
Arguments:
sbml: an sbml Model instance
filter: callback function taking assignment variable as input and returning True/False to indicate if the respective rule should be turned into an observable
sbml_model: an sbml Model instance
filter_function: callback function taking assignment variable as input and returning True/False to indicate if the respective rule should be turned into an observable
Returns:
A dictionary(observableId:{'name':observableNamem,'formula':formulaString})
Expand All @@ -1494,14 +1500,14 @@ def assignmentRules2observables(sbml, filter = lambda *_: True):
"""
observables = {}
for p in sbml.getListOfParameters():
for p in sbml_model.getListOfParameters():
parameterId = p.getId()
if filter(parameterId):
observables[parameterId] = {'name': p.getName(), 'formula': sbml.getAssignmentRuleByVariable(parameterId).getFormula()}
if filter_function(parameterId):
observables[parameterId] = {'name': p.getName(), 'formula': sbml_model.getAssignmentRuleByVariable(parameterId).getFormula()}

for parameterId in observables:
sbml.removeRuleByVariable(parameterId)
sbml.removeParameter(parameterId)
sbml_model.removeRuleByVariable(parameterId)
sbml_model.removeParameter(parameterId)

return observables

Expand Down Expand Up @@ -1545,4 +1551,47 @@ def constantSpeciesToParameters(sbml_model):
reaction.removeProduct(speciesId)

return transformable

def replaceLogAB(x):
"""
Replace log(a, b) in the given string by ln(b)/ln(a)
Works for nested parentheses and nested 'log's. This can be used to circumvent
the incompatible argument order between symengine (log(x, basis)) and libsbml
(log(basis, x)).
Arguments:
x: string to replace
Returns:
string with replaced 'log's
Raises:
"""

match = re.search(r'(^|\W)log\(', x)
if not match:
return x

# index of 'l' of 'log'
logStart = match.start() if match.end() - match.start() == 4 else match.start() + 1
level = 0 # parenthesis level
posComma = -1 # position of comma in log(a,b)
for i in range(logStart + 4, len(x)):
if x[i] == '(':
level += 1
elif x[i] == ')':
level -= 1
if level == -1: break
elif x[i] == ',' and level == 0:
posComma = i

if posComma < 0:
# was log(a), not log(a,b), so nothing to replace
return x

replacement = '{prefix}ln({a})/ln({basis}){suffix}'.format(prefix = x[:logStart],
suffix = x[i+1:],
basis = x[logStart+4 : posComma],
a = x[posComma+1 : i])
return replaceLogAB(replacement)

0 comments on commit c039ffd

Please sign in to comment.