CellMLToNektar.optimize.LinearityAnalyser Class Reference

Inheritance diagram for CellMLToNektar.optimize.LinearityAnalyser:

Public Member Functions
def	analyse_for_jacobian (self, doc, V=None)
def	find_linear_odes (self, state_vars, V, free_var)
def	rearrange_linear_odes (self, doc)
def	show (self, d)

Static Public Attributes
	LINEAR_KINDS = Enum('None', 'Linear', 'Nonlinear')

Private Member Functions
def	_get_rhs (self, expr)
def	_check_expr (self, expr, state_var, bad_vars)
def	_clone (self, expr)
def	_make_apply (self, operator, ghs, i, filter_none=True, preserve=False)
def	_transfer_lut (self, expr, gh, var)
def	_rearrange_expr (self, expr, var)

Private Attributes
	__expr

Detailed Description

Analyse linearity aspects of a model.

This class performs analyses to determine which ODEs have a linear
dependence on their state variable, discounting the presence of
the transmembrane potential.

This can be used to decouple the ODEs for more efficient solution,
especially in a multi-cellular context.

analyse_for_jacobian(doc) must be called before
rearrange_linear_odes(doc).

Definition at line 903 of file optimize.py.

Member Function Documentation

_check_expr()

def CellMLToNektar.optimize.LinearityAnalyser._check_expr (   self,   expr,   state_var,   bad_vars  )

private

The actual linearity checking function.

Recursively determine the type of dependence expr has on
state_var.  The presence of any members of bad_vars indicates
a non-linear dependency.

Return a member of the self.LINEAR_KINDS enum.

Definition at line 959 of file optimize.py.

References CellMLToNektar.optimize.LinearityAnalyser._check_expr(), CellMLToNektar.optimize.LinearityAnalyser._get_rhs(), CellMLToNektar.pycml.child_i(), CellMLToNektar.utilities.DEBUG(), and CellMLToNektar.optimize.LinearityAnalyser.LINEAR_KINDS.

Referenced by CellMLToNektar.optimize.LinearityAnalyser._check_expr(), and CellMLToNektar.optimize.LinearityAnalyser.find_linear_odes().

    def _check_expr(self, expr, state_var, bad_vars):
        """The actual linearity checking function.

        Recursively determine the type of dependence expr has on
        state_var.  The presence of any members of bad_vars indicates
        a non-linear dependency.

        Return a member of the self.LINEAR_KINDS enum.
        """
        kind = self.LINEAR_KINDS
        result = None
        if isinstance(expr, mathml_ci):
            var = expr.variable.get_source_variable(recurse=True)
            if var is state_var:
                result = kind.Linear
            elif var in bad_vars:
                result = kind.Nonlinear
            elif var.get_type(follow_maps=True) == VarTypes.Computed:
                # Recurse into defining expression
                src_var = var.get_source_variable(recurse=True)
                src_expr = self._get_rhs(src_var.get_dependencies()[0])
                DEBUG('find-linear-deps', "--recurse for", src_var.name,
                      "to", src_expr)
                result = self._check_expr(src_expr, state_var, bad_vars)
            else:
                result = kind.None
            # Record the kind of this variable, for later use when
            # rearranging linear ODEs
            var._cml_linear_kind = result
        elif isinstance(expr, mathml_cn):
            result = kind.None
        elif isinstance(expr, mathml_piecewise):
            # If any conditions have a dependence, then we're
            # nonlinear.  Otherwise, all the pieces must be the same
            # (and that's what we are) or we're nonlinear.
            pieces = getattr(expr, u'piece', [])
            conds = map(lambda p: child_i(p, 2), pieces)
            chld_exprs = map(lambda p: child_i(p, 1), pieces)
            if hasattr(expr, u'otherwise'):
                chld_exprs.append(child_i(expr.otherwise, 1))
            for cond in conds:
                if self._check_expr(cond, state_var, bad_vars) != kind.None:
                    result = kind.Nonlinear
                    break
            else:
                # Conditions all OK
                for e in chld_exprs:
                    res = self._check_expr(e, state_var, bad_vars)
                    if result is not None and res != result:
                        # We have a difference
                        result = kind.Nonlinear
                        break
                    result = res
        elif isinstance(expr, mathml_apply):
            # Behaviour depends on the operator
            operator = expr.operator().localName
            operands = expr.operands()
            if operator in ['plus', 'minus']:
                # Linear if any operand linear, and none non-linear
                op_kinds = map(lambda op: self._check_expr(op, state_var,
                                                           bad_vars),
                               operands)
                result = max(op_kinds)
            elif operator == 'divide':
                # Linear iff only numerator linear
                numer = operands.next()
                denom = operands.next()
                if self._check_expr(denom, state_var, bad_vars) != kind.None:
                    result = kind.Nonlinear
                else:
                    result = self._check_expr(numer, state_var, bad_vars)
            elif operator == 'times':
                # Linear iff only 1 linear operand
                op_kinds = map(lambda op: self._check_expr(op, state_var,
                                                           bad_vars),
                               operands)
                lin, nonlin = 0, 0
                for res in op_kinds:
                    if res == kind.Linear: lin += 1
                    elif res == kind.Nonlinear: nonlin += 1
                if nonlin > 0 or lin > 1:
                    result = kind.Nonlinear
                elif lin == 1:
                    result = kind.Linear
                else:
                    result = kind.None
            else:
                # Cannot be linear; may be no dependence at all
                result = max(map(lambda op: self._check_expr(op, state_var,
                                                             bad_vars),
                                 operands))
                if result == kind.Linear:
                    result = kind.Nonlinear
        else:
            # Either a straightforward container element
            try:
                child = child_i(expr, 1)
            except ValueError:
                # Assume it's just a constant
                result = kind.None
            else:
                result = self._check_expr(child, state_var, bad_vars)
        DEBUG('find-linear-deps', "Expression", expr, "gives result", result)
        return result

_clone()

def CellMLToNektar.optimize.LinearityAnalyser._clone (   self,   expr  )

private

Properly clone a MathML sub-expression.

Definition at line 1084 of file optimize.py.

References CellMLToNektar.optimize.LinearityAnalyser._make_apply().

Referenced by CellMLToNektar.optimize.LinearityAnalyser._make_apply(), and CellMLToNektar.optimize.LinearityAnalyser._rearrange_expr().

    def _clone(self, expr):
        """Properly clone a MathML sub-expression."""
        if isinstance(expr, mathml):
            clone = expr.clone_self(register=True)
        else:
            clone = mathml.clone(expr)
        return clone

_get_rhs()

def CellMLToNektar.optimize.LinearityAnalyser._get_rhs (   self,   expr  )

private

Return the RHS of an assignment expression.

Definition at line 953 of file optimize.py.

Referenced by CellMLToNektar.optimize.LinearityAnalyser._check_expr(), CellMLToNektar.optimize.LinearityAnalyser._rearrange_expr(), CellMLToNektar.optimize.LinearityAnalyser.find_linear_odes(), and CellMLToNektar.optimize.LinearityAnalyser.rearrange_linear_odes().

    def _get_rhs(self, expr):
        """Return the RHS of an assignment expression."""
        ops = expr.operands()
        ops.next()
        return ops.next()

_make_apply()

def CellMLToNektar.optimize.LinearityAnalyser._make_apply (   self,   operator,   ghs,   i,   filter_none = True,   preserve = False  )

private

Construct a new apply expression for g or h.

ghs is an iterable of (g,h) pairs for operands.

i indicates whether to construct g (0) or h (1).

filter_none indicates the behaviour of 0 under this operator.
If True, it's an additive zero, otherwise it's a
multiplicative zero.

Definition at line 1093 of file optimize.py.

References CellMLToNektar.optimize.LinearityAnalyser.__expr, and CellMLToNektar.optimize.LinearityAnalyser._clone().

Referenced by CellMLToNektar.optimize.LinearityAnalyser._clone().

                    preserve=False):
        """Construct a new apply expression for g or h.

        ghs is an iterable of (g,h) pairs for operands.
        
        i indicates whether to construct g (0) or h (1).
        
        filter_none indicates the behaviour of 0 under this operator.
        If True, it's an additive zero, otherwise it's a
        multiplicative zero.
        """
        # Find g or h operands
        ghs_i = map(lambda gh: gh[i], ghs)
        if not filter_none and None in ghs_i:
            # Whole expr is None
            new_expr = None
        else:
            # Filter out None subexprs
            if operator == u'minus':
                # Do we need to retain a unary minus?
                if len(ghs_i) == 1 or ghs_i[0] is None:
                    # Original was -a or 0-a
                    retain_unary_minus = True
                else:
                    # Original was a-0 or a-b
                    retain_unary_minus = False
            else:
                # Only retain if we're told to preserve as-is
                retain_unary_minus = preserve
            ghs_i = filter(None, ghs_i)
            if ghs_i:
                if len(ghs_i) > 1 or retain_unary_minus:
                    new_expr = mathml_apply.create_new(
                        self.__expr, operator, ghs_i)
                else:
                    new_expr = self._clone(ghs_i[0]) # Clone may be unneeded
            else:
                new_expr = None
        return new_expr

_rearrange_expr()

def CellMLToNektar.optimize.LinearityAnalyser._rearrange_expr (   self,   expr,   var  )

private

Rearrange an expression into the form g + h*var.

Performs a post-order traversal of this expression's tree,
and returns a pair (g, h)

Definition at line 1156 of file optimize.py.

References CellMLToNektar.optimize.LinearityAnalyser._clone(), CellMLToNektar.optimize.LinearityAnalyser._get_rhs(), CellMLToNektar.optimize.LinearityAnalyser._rearrange_expr(), CellMLToNektar.optimize.LinearityAnalyser._transfer_lut(), CellMLToNektar.pycml.child_i(), and CellMLToNektar.optimize.LinearityAnalyser.LINEAR_KINDS.

Referenced by CellMLToNektar.optimize.LinearityAnalyser._rearrange_expr(), and CellMLToNektar.optimize.LinearityAnalyser.rearrange_linear_odes().

    def _rearrange_expr(self, expr, var):
        """Rearrange an expression into the form g + h*var.

        Performs a post-order traversal of this expression's tree,
        and returns a pair (g, h)
        """
#        import inspect
#        depth = len(inspect.stack())
#        print ' '*depth, "_rearrange_expr", prid(expr, True), var.name, expr
        gh = None
        if isinstance(expr, mathml_ci):
            # Variable
            ci_var = expr.variable.get_source_variable(recurse=True)
            if var is ci_var:
                gh = (None, mathml_cn.create_new(expr,
                                                 u'1', u'dimensionless'))
            else:
                if ci_var._cml_linear_kind == self.LINEAR_KINDS.None:
                    # Just put the <ci> in g, but with full name
                    gh = (mathml_ci.create_new(expr, ci_var.fullname()), None)
                    gh[0]._set_variable_obj(ci_var)
                else:
                    # ci_var is a linear function of var, so rearrange
                    # its definition
                    if not hasattr(ci_var, '_cml_linear_split'):
                        ci_defn = ci_var.get_dependencies()[0]
                        ci_var._cml_linear_split = self._rearrange_expr(
                            self._get_rhs(ci_defn), var)
                    gh = ci_var._cml_linear_split
        elif isinstance(expr, mathml_piecewise):
            # The tests have to move into both components of gh:
            # "if C1 then (a1,b1) elif C2 then (a2,b2) else (a0,b0)"
            # maps to "(if C1 then a1 elif C2 then a2 else a0,
            #           if C1 then b1 elif C2 then b2 else b0)"
            # Note that no test is a function of var.
            # First rearrange child expressions
            pieces = getattr(expr, u'piece', [])
            cases = map(lambda p: child_i(p, 1), pieces)
            cases_ghs = map(lambda c: self._rearrange_expr(c, var), cases)
            if hasattr(expr, u'otherwise'):
                ow_gh = self._rearrange_expr(child_i(expr.otherwise, 1), var)
            else:
                ow_gh = (None, None)
            # Now construct the new expression
            conds = map(lambda p: self._clone(child_i(p, 2)), pieces)
            def piecewise_branch(i):
                pieces_i = zip(map(lambda gh: gh[i], cases_ghs),
                               conds)
                pieces_i = filter(lambda p: p[0] is not None,
                                  pieces_i) # Remove cases that are None
                ow = ow_gh[i]
                if pieces_i:
                    new_expr = mathml_piecewise.create_new(
                        expr, pieces_i, ow)
                elif ow:
                    new_expr = ow
                else:
                    new_expr = None
                return new_expr
            gh = (piecewise_branch(0), piecewise_branch(1))
            self._transfer_lut(expr, gh, var)
        elif isinstance(expr, mathml_apply):
            # Behaviour depends on the operator
            operator = expr.operator().localName
            operands = expr.operands()
            self.__expr = expr # For self._make_apply
            if operator in ['plus', 'minus']:
                # Just split the operation into each component
                operand_ghs = map(lambda op: self._rearrange_expr(op, var),
                                  operands)
                g = self._make_apply(operator, operand_ghs, 0)
                h = self._make_apply(operator, operand_ghs, 1)
                gh = (g, h)
            elif operator == 'divide':
                # (a, b) / (c, 0) = (a/c, b/c)
                numer = self._rearrange_expr(operands.next(), var)
                denom = self._rearrange_expr(operands.next(), var)
                assert denom[1] is None
                denom_g = denom[0]
                g = h = None
                if numer[0]:
                    g = mathml_apply.create_new(expr, operator,
                                                [numer[0], denom_g])
                if numer[1]:
                    if g:
                        denom_g = self._clone(denom_g)
                    h = mathml_apply.create_new(expr, operator,
                                                [numer[1], denom_g])
                gh = (g, h)
            elif operator == 'times':
                # (a1,b1)*(a2,b2) = (a1*a2, b1*a2 or a1*b2 or None)
                # Similarly for the nary case - at most one b_i is not None
                operand_ghs = map(lambda op: self._rearrange_expr(op, var),
                                  operands)
                g = self._make_apply(operator, operand_ghs, 0,
                                     filter_none=False)
                # Find non-None b_i, if any
                for i, ab in enumerate(operand_ghs):
                    if ab[1] is not None:
                        operand_ghs[i] = (ab[1], None)
                        # Clone the a_i to avoid objects having 2 parents
                        for j, ab in enumerate(operand_ghs):
                            if j != i:
                                operand_ghs[j] = (self._clone(operand_ghs[j][0]), None)
                        h = self._make_apply(operator, operand_ghs, 0, filter_none=False)
                        break
                else:
                    h = None
                gh = (g, h)
            else:
                # (a, None) op (b, None) = (a op b, None)
                operand_ghs = map(lambda op: self._rearrange_expr(op, var),
                                  operands)
                g = self._make_apply(operator, operand_ghs, 0, preserve=True)
                gh = (g, None)
            self._transfer_lut(expr, gh, var)
        else:
            # Since this expression is linear, there can't be any
            # occurrence of var in it; all possible such cases are covered
            # above.  So just clone it into g.
            gh = (self._clone(expr), None)
#        print ' '*depth, "Re-arranged", prid(expr, True), "to", prid(gh[0], True), ",", prid(gh[1], True)
        return gh
    

_transfer_lut()

def CellMLToNektar.optimize.LinearityAnalyser._transfer_lut (   self,   expr,   gh,   var  )

private

Transfer lookup table annotations from expr to gh.

gh is a pair (g, h) s.t. expr = g + h*var.

If expr can be represented by a lookup table, then the lookup
variable cannot be var, since if it were, then expr would have
a non-linear dependence on var.  Hence h must be 0, since
otherwise expr would contain a (state) variable other than the
lookup variable, and hence not be a candidate for a table.
Thus expr=g, so we transfer the annotations to g.

Definition at line 1133 of file optimize.py.

Referenced by CellMLToNektar.optimize.LinearityAnalyser._rearrange_expr().

    def _transfer_lut(self, expr, gh, var):
        """Transfer lookup table annotations from expr to gh.

        gh is a pair (g, h) s.t. expr = g + h*var.

        If expr can be represented by a lookup table, then the lookup
        variable cannot be var, since if it were, then expr would have
        a non-linear dependence on var.  Hence h must be 0, since
        otherwise expr would contain a (state) variable other than the
        lookup variable, and hence not be a candidate for a table.
        Thus expr=g, so we transfer the annotations to g.
        """
        if expr.getAttributeNS(NSS['lut'], u'possible', '') != u'yes':
            return
        # Paranoia check that our reasoning is correct
        g, h = gh
        assert h is None
        # Transfer the annotations into g
        LookupTableAnalyser.copy_lut_annotations(expr, g)
        # Make sure g has a reference to its component, for use by code generation.
        g._cml_component = expr.component
        return

analyse_for_jacobian()

def CellMLToNektar.optimize.LinearityAnalyser.analyse_for_jacobian	(		self,
			doc,
			V = None
	)

Analyse the model for computing a symbolic Jacobian.

Determines automatically which variables will need to be solved
for using Newton's method, and stores their names in
doc.model._cml_nonlinear_system_variables, as a list of variable
objects.

Also stores doc.model._cml_linear_vars, doc.model._cml_free_var,
doc.model._cml_transmembrane_potential.

Definition at line 918 of file optimize.py.

References CellMLToNektar.utilities.DEBUG(), and CellMLToNektar.optimize.LinearityAnalyser.find_linear_odes().

    def analyse_for_jacobian(self, doc, V=None):
        """Analyse the model for computing a symbolic Jacobian.

        Determines automatically which variables will need to be solved
        for using Newton's method, and stores their names in
        doc.model._cml_nonlinear_system_variables, as a list of variable
        objects.

        Also stores doc.model._cml_linear_vars, doc.model._cml_free_var,
        doc.model._cml_transmembrane_potential.
        """
        # TODO: Add error checking and tidy
        stvs = doc.model.find_state_vars()
        if V is None:
            Vcname, Vvname = 'membrane', 'V'
            V = doc.model.get_variable_by_name(Vcname, Vvname)
        V = V.get_source_variable(recurse=True)
        doc.model._cml_transmembrane_potential = V
        free_var = doc.model.find_free_vars()[0]
        lvs = self.find_linear_odes(stvs, V, free_var)
        # Next 3 lines for benefit of rearrange_linear_odes(doc)
        lvs.sort(key=lambda v: v.fullname())
        doc.model._cml_linear_vars = lvs
        doc.model._cml_free_var = free_var
        # Store nonlinear vars in canonical order
        nonlinear_vars = list(set(stvs) - set([V]) - set(lvs))
        nonlinear_vars.sort(key=lambda v: v.fullname())
        doc.model._cml_nonlinear_system_variables = nonlinear_vars
        # Debugging
        f = lambda var: var.fullname()
        DEBUG('linearity-analyser', 'V=', V.fullname(), '; free var=',
              free_var.fullname(), '; linear vars=', map(f, lvs),
              '; nonlinear vars=', map(f, nonlinear_vars))
        return

find_linear_odes()

def CellMLToNektar.optimize.LinearityAnalyser.find_linear_odes	(		self,
			state_vars,
			V,
			free_var
	)

Identify linear ODEs.

For each ODE (except that for V), determine whether it has a
linear dependence on the dependent variable, and thus can be
updated directly, without using Newton's method.

We also require it to not depend on any other state variable,
except for V.

Definition at line 1064 of file optimize.py.

References CellMLToNektar.optimize.LinearityAnalyser._check_expr(), CellMLToNektar.optimize.LinearityAnalyser._get_rhs(), and CellMLToNektar.optimize.LinearityAnalyser.LINEAR_KINDS.

Referenced by CellMLToNektar.optimize.LinearityAnalyser.analyse_for_jacobian().

    def find_linear_odes(self, state_vars, V, free_var):
        """Identify linear ODEs.

        For each ODE (except that for V), determine whether it has a
        linear dependence on the dependent variable, and thus can be
        updated directly, without using Newton's method.

        We also require it to not depend on any other state variable,
        except for V.
        """
        kind = self.LINEAR_KINDS
        candidates = set(state_vars) - set([V])
        linear_vars = []
        for var in candidates:
            ode_expr = var.get_ode_dependency(free_var)
            if self._check_expr(self._get_rhs(ode_expr), var,
                                candidates - set([var])) == kind.Linear:
                linear_vars.append(var)
        return linear_vars

rearrange_linear_odes()

def CellMLToNektar.optimize.LinearityAnalyser.rearrange_linear_odes	(		self,
			doc
	)

Rearrange the linear ODEs so they can be updated directly
on solving.

Each ODE du/dt = f(u, t) can be written in the form
    du/dt = g(t) + h(t)u.
A backward Euler update step is then as simple as
    u_n = (u_{n-1} + g(t)dt) / (1 - h(t)dt)
(assuming that the transmembrane potential has already been
solved for at t_n.

Stores the results in doc.model._cml_linear_update_exprs, a
mapping from variable object u to pair (g, h).

Definition at line 1280 of file optimize.py.

References CellMLToNektar.optimize.LinearityAnalyser._get_rhs(), and CellMLToNektar.optimize.LinearityAnalyser._rearrange_expr().

    def rearrange_linear_odes(self, doc):
        """Rearrange the linear ODEs so they can be updated directly
        on solving.

        Each ODE du/dt = f(u, t) can be written in the form
            du/dt = g(t) + h(t)u.
        A backward Euler update step is then as simple as
            u_n = (u_{n-1} + g(t)dt) / (1 - h(t)dt)
        (assuming that the transmembrane potential has already been
        solved for at t_n.

        Stores the results in doc.model._cml_linear_update_exprs, a
        mapping from variable object u to pair (g, h).
        """

        odes = map(lambda v: v.get_ode_dependency(doc.model._cml_free_var),
                   doc.model._cml_linear_vars)
        result = {}
        for var, ode in itertools.izip(doc.model._cml_linear_vars, odes):
            # Do the re-arrangement for this variable.
            # We do this by a post-order traversal of its defining ODE,
            # constructing a pair (g, h) for each subexpression recursively.
            # First, get the RHS of the ODE
            rhs = self._get_rhs(ode)
            # And traverse
            result[var] = self._rearrange_expr(rhs, var)
        # Store result in model
        doc.model._cml_linear_update_exprs = result
        return result

show()

def CellMLToNektar.optimize.LinearityAnalyser.show	(		self,
			d
	)

Print out a more readable report on a rearrangement,
as given by self.rearrange_linear_odes.

Definition at line 1310 of file optimize.py.

    def show(self, d):
        """Print out a more readable report on a rearrangement,
        as given by self.rearrange_linear_odes."""
        for var, expr in d.iteritems():
            print var.fullname()
            print "G:", expr[0].xml()
            print "H:", expr[1].xml()
            print "ODE:", var.get_ode_dependency(
                var.model._cml_free_var).xml()
            print

Member Data Documentation

__expr

CellMLToNektar.optimize.LinearityAnalyser.__expr

private

Definition at line 1221 of file optimize.py.

Referenced by CellMLToNektar.optimize.LinearityAnalyser._make_apply().

LINEAR_KINDS

CellMLToNektar.optimize.LinearityAnalyser.LINEAR_KINDS = Enum('None', 'Linear', 'Nonlinear')

static

Definition at line 916 of file optimize.py.

Referenced by CellMLToNektar.optimize.LinearityAnalyser._check_expr(), CellMLToNektar.optimize.LinearityAnalyser._rearrange_expr(), and CellMLToNektar.optimize.LinearityAnalyser.find_linear_odes().