Merge 969178b into 4f9f4a1

pysb/__init__.py

@@ -3,7 +3,7 @@
 
 __all__ = ['Observable', 'Initial', 'MatchOnce', 'Model', 'Monomer',
            'Parameter', 'Compartment', 'Rule', 'Expression', 'ANY', 'WILD',
-           'Annotation', 'MultiState']
+           'Annotation', 'MultiState', 'Tag']
 
 try:
     import reinteract         # fails if reinteract not installed

pysb/bng.py

@@ -15,7 +15,6 @@
 import collections
 import pysb.pathfinder as pf
 from pysb.logging import get_logger, EXTENDED_DEBUG
-import logging
 
 try:
     from cStringIO import StringIO
@@ -735,6 +734,13 @@ def generate_equations(model, cleanup=True, verbose=False, **kwargs):
 def _parse_netfile(model, lines):
     """ Parse species, rxns and groups from a BNGL net file """
     try:
+        while 'begin parameters' not in next(lines):
+            pass
+        while True:
+            line = next(lines)
+            if 'end parameters' in line: break
+            _parse_parameter(model, line)
+
         while 'begin species' not in next(lines):
             pass
         model.species = []
@@ -770,6 +776,68 @@ def _parse_netfile(model, lines):
         pass
 
 
+def _parse_parameter(model, line):
+    index, pname, pval, hash, ptype = line.strip().split()
+    par_names = model.components.keys()
+    if pname not in par_names:
+        if ptype == 'Constant' and pname not in model._derived_parameters.keys():
+            try:
+                p = pysb.core.Parameter(pname, pval, _export=False)
+            except ValueError:
+                p = pysb.core.Parameter(pname,
+                                        eval(pval.replace('^', '**')),
+                                        _export=False)
+            model._derived_parameters.add(p)
+        elif ptype == 'ConstantExpression' and \
+                pname not in model._derived_expressions.keys():
+            pval = _fix_booleans(pval)
+            p = pysb.core.Expression(pname, sympy.sympify(pval),
+                                     _export=False)
+            model._derived_expressions.add(p)
+        else:
+            raise ValueError('Unknown type {} for parameter {}'.format(
+                ptype, pname))
+
+
+_RE_NUMBER = '\d+(?:\.\d+)?'
+_RE_EQUAL_INEQ = '{0}\s*([<>]=?|!=|==)\s*{0}'.format(_RE_NUMBER)
+_RE_AND = '({0})\s*&&\s*({0})'.format(_RE_NUMBER)
+_RE_OR = '({0})\s*\|\|\s*({0})'.format(_RE_NUMBER)
+_RE_NUMBER_PARENS = '\(\s*({0})\s*\)'.format(_RE_NUMBER)
+
+
+def _fix_booleans(pval):
+    """ Evaluate boolean expressions to integers
+
+    Inequalities (e.g. "(0>1)") don't parse to integers in sympy, so expressions
+    from BNG like "(1>0)*10" cause a sympy parse error. To avoid this, we
+    replace the boolean expressions using a regex before calling sympy.
+    """
+    while True:
+        orig = pval
+
+        for match in re.finditer(_RE_OR, pval):
+            pval = pval.replace(
+                match.group(0),
+                '1' if float(match.group(1)) or float(match.group(2)) else '0'
+            )
+        for match in re.finditer(_RE_AND, pval):
+            pval = pval.replace(
+                match.group(0),
+                '1' if float(match.group(1)) and float(match.group(2)) else '0'
+            )
+        for match in re.finditer(_RE_EQUAL_INEQ, pval):
+            pval = pval.replace(
+                match.group(0),
+                '1' if eval(match.group(0)) else '0'
+            )
+        for match in re.finditer(_RE_NUMBER_PARENS, pval):
+            pval = pval.replace(match.group(0), match.group(1))
+
+        if orig == pval:
+            return pval
+
+
 def _parse_species(model, line):
     """Parse a 'species' line from a BNGL net file."""
     index, species, value = line.strip().split()
@@ -807,7 +875,7 @@ def _parse_species(model, line):
                     site_name, condition = ss, None
                 site_conditions[site_name].append(condition)
 
-        site_conditions = {k: v[0] if len(v) == 1 else tuple(v)
+        site_conditions = {k: v[0] if len(v) == 1 else pysb.core.MultiState(*v)
                            for k, v in site_conditions.items()}
         monomer = model.monomers[monomer_name]
         monomer_compartment = model.compartments.get(monomer_compartment_name)
@@ -839,7 +907,8 @@ def _parse_reaction(model, line, reaction_cache):
     is_reverse = tuple(bool(i) for i in is_reverse)
     r_names = ['__s%d' % r for r in reactants]
     rate_param = [model.parameters.get(r) or model.expressions.get(r) or
-                  float(r) for r in rate]
+                  model._derived_parameters.get(r) or
+                  model._derived_expressions.get(r) or float(r) for r in rate]
     combined_rate = sympy.Mul(*[sympy.S(t) for t in r_names + rate_param])
     reaction = {
         'reactants': reactants,

pysb/builder.py

@@ -178,6 +178,12 @@ def initial(self, *args, **kwargs):
         self.model.add_initial(i)
         return i
 
+    def tag(self, *args, **kwargs):
+        """Adds a tag to the Builder's model instance."""
+        t = Tag(*args, _export=False, **kwargs)
+        self.model.add_component(t)
+        return t
+
     def __getitem__(self, index):
         """Returns the component with the given string index
         from the instance of the model contained by the Builder."""

pysb/core.py

@@ -517,6 +517,7 @@ def __init__(self, monomer, site_conditions, compartment):
         self.site_conditions = site_conditions
         self.compartment = compartment
         self._graph = None
+        self._tag = None
 
     def is_concrete(self):
         """
@@ -595,7 +596,9 @@ def __call__(self, conditions=None, **kwargs):
         # updated according to our args (as in Monomer.__call__).
         site_conditions = self.site_conditions.copy()
         site_conditions.update(extract_site_conditions(conditions, **kwargs))
-        return MonomerPattern(self.monomer, site_conditions, self.compartment)
+        mp = MonomerPattern(self.monomer, site_conditions, self.compartment)
+        mp._tag = self._tag
+        return mp
 
     def __add__(self, other):
         if isinstance(other, MonomerPattern):
@@ -649,6 +652,18 @@ def __pow__(self, other):
         else:
             return NotImplemented
 
+    def __matmul__(self, other):
+        if not isinstance(other, Tag):
+            return NotImplemented
+
+        if self._tag is not None:
+            raise TagAlreadySpecifiedError()
+
+        # Need to upgrade to a ComplexPattern
+        cp_new = as_complex_pattern(self)
+        cp_new._tag = other
+        return cp_new
+
     def __repr__(self):
         value = '%s(' % self.monomer.name
         sites_unique = list(collections.OrderedDict.fromkeys(
@@ -661,6 +676,8 @@ def __repr__(self):
         value += ')'
         if self.compartment is not None:
             value += ' ** ' + self.compartment.name
+        if self._tag is not None:
+            value = '{} @ {}'.format(self._tag.name, value)
         return value
 
 
@@ -698,6 +715,7 @@ def __init__(self, monomer_patterns, compartment, match_once=False):
         self.compartment = compartment
         self.match_once = match_once
         self._graph = None
+        self._tag = None
 
     def is_concrete(self):
         """
@@ -988,6 +1006,8 @@ def __radd__(self, other):
             return NotImplemented
 
     def __mod__(self, other):
+        if self._tag is not None:
+            raise ValueError('Tag should be specified at the end of the complex')
         if isinstance(other, MonomerPattern):
             return ComplexPattern(self.monomer_patterns + [other], self.compartment, self.match_once)
         elif isinstance(other, ComplexPattern):
@@ -1031,12 +1051,25 @@ def __pow__(self, other):
         else:
             return NotImplemented
 
+    def __matmul__(self, other):
+        if not isinstance(other, Tag):
+            return NotImplemented
+
+        if self._tag is not None:
+            raise TagAlreadySpecifiedError()
+
+        cp_new = self.copy()
+        cp_new._tag = other
+        return cp_new
+
     def __repr__(self):
         ret = ' % '.join([repr(p) for p in self.monomer_patterns])
         if self.compartment is not None:
             ret = '(%s) ** %s' % (ret, self.compartment.name)
         if self.match_once:
             ret = 'MatchOnce(%s)' % ret
+        if self._tag:
+            ret = '{} @ {}'.format(ret, self._tag.name)
         return ret
 
 
@@ -1358,6 +1391,44 @@ def __init__(self, name, rule_expression, rate_forward, rate_reverse=None,
 
         Component.__init__(self, name, _export)
 
+        # Get tags from rule expression
+        tags = set()
+        for rxn_pat in (rule_expression.reactant_pattern,
+                        rule_expression.product_pattern):
+            if rxn_pat.complex_patterns:
+                for cp in rxn_pat.complex_patterns:
+                    if cp is not None:
+                        if cp._tag:
+                            tags.add(cp._tag)
+                        tags.update(mp._tag for mp in cp.monomer_patterns
+                                    if mp._tag is not None)
+
+        # Check that tags defined in rates are used in the expression
+        tags_rates = (self._check_rate_tags('forward', tags) +
+                      self._check_rate_tags('reverse', tags))
+
+        missing = tags.difference(set(tags_rates))
+        if missing:
+            names = [t.name for t in missing]
+            warnings.warn(
+                'Rule "{}": Tags {} defined in rule expression but not used in '
+                'rates'.format(self.name, ', '.join(names)), UserWarning)
+
+    def _check_rate_tags(self, direction, tags):
+        rate = self.rate_forward if direction == 'forward' else \
+            self.rate_reverse
+        if not isinstance(rate, Expression):
+            return []
+        tags_rate = rate.tags()
+        missing = set(tags_rate).difference(tags)
+        if missing:
+            names = [t.name for t in missing]
+            raise ValueError(
+                'Rule "{}": Tag(s) {} defined in {} rate but not in '
+                'expression'.format(self.name, ', '.join(names), direction))
+
+        return tags_rate
+
     def is_synth(self):
         """Return a bool indicating whether this is a synthesis rule."""
         return len(self.reactant_pattern.complex_patterns) == 0 or \
@@ -1485,6 +1556,13 @@ def __repr__(self):
     def __str__(self):
         return repr(self)
 
+    def __call__(self, tag):
+        if not isinstance(tag, Tag):
+            raise ValueError('Observables are only callable with a Tag '
+                             'instance, for use within local Expressions')
+
+        return sympy.Function(self.name)(tag)
+
 
 class Expression(Component, sympy.Symbol):
 
@@ -1543,6 +1621,13 @@ def get_value(self):
     def func(self):
         return sympy.Symbol
 
+    @property
+    def is_local(self):
+        return len(self.expr.atoms(Tag)) > 0
+
+    def tags(self):
+        return sorted(self.expr.atoms(Tag), key=lambda tag: tag.name)
+
     def __repr__(self):
         ret = '%s(%s, %s)' % (self.__class__.__name__, repr(self.name),
                               repr(self.expr))
@@ -1551,6 +1636,39 @@ def __repr__(self):
     def __str__(self):
         return repr(self)
 
+    def __call__(self, tag):
+        if not isinstance(tag, Tag):
+            raise ValueError('Expressions are only callable with a Tag '
+                             'instance, for use within local Expressions')
+
+        return sympy.Function(self.name)(tag)
+
+
+class Tag(Component, sympy.Symbol):
+    """ Tag for labelling MonomerPatterns and ComplexPatterns """
+    def __new__(cls, name, _export=True):
+        return super(sympy.Symbol, cls).__new__(cls, name)
+
+    def __getnewargs__(self):
+        return self.name, False
+
+    def __init__(self, name, _export=True):
+        Component.__init__(self, name, _export)
+
+    def __matmul__(self, other):
+        if not isinstance(other, MonomerPattern):
+            return NotImplemented
+
+        if other._tag is not None:
+            raise TagAlreadySpecifiedError()
+
+        new_mp = other()
+        new_mp._tag = self
+        return new_mp
+
+    def __repr__(self):
+        return "{}({})".format(self.__class__.__name__, repr(self.name))
+
 
 class Initial(object):
     """
@@ -1659,7 +1777,7 @@ class Model(object):
     """
 
     _component_types = (Monomer, Compartment, Parameter, Rule, Observable,
-                        Expression)
+                        Expression, Tag)
 
     def __init__(self, name=None, base=None, _export=True):
         self.name = name
@@ -1671,6 +1789,7 @@ def __init__(self, name=None, base=None, _export=True):
         self.rules = ComponentSet()
         self.observables = ComponentSet()
         self.expressions = ComponentSet()
+        self.tags = ComponentSet()
         self.initials = []
         self.annotations = []
         self._odes = OdeView(self)
@@ -1820,9 +1939,12 @@ def expressions_constant(self):
                             if e.is_constant_expression())
         return cset
 
-    def expressions_dynamic(self):
+    def expressions_dynamic(self, include_local=True):
         """Return a ComponentSet of non-constant expressions."""
-        return self.expressions - self.expressions_constant()
+        cset = self.expressions - self.expressions_constant()
+        if not include_local:
+            cset = ComponentSet(e for e in cset if not e.is_local)
+        return cset
 
     @property
     def odes(self):
@@ -1990,6 +2112,8 @@ def reset_equations(self):
         self.reactions = []
         self.reactions_bidirectional = []
         self._stoichiometry_matrix = None
+        self._derived_parameters = ComponentSet()
+        self._derived_expressions = ComponentSet()
         for obs in self.observables:
             obs.species = []
             obs.coefficients = []
@@ -2053,6 +2177,11 @@ class UnknownSiteError(ValueError):
 class CompartmentAlreadySpecifiedError(ValueError):
     pass
 
+
+class TagAlreadySpecifiedError(ValueError):
+    pass
+
+
 class ModelNotDefinedError(RuntimeError):
     """SelfExporter method was called before a model was defined."""
     def __init__(self):