Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Public Member Functions | Static Public Attributes | Private Member Functions | Private Attributes | List of all members
CellMLToNektar.optimize.LinearityAnalyser Class Reference
Inheritance diagram for CellMLToNektar.optimize.LinearityAnalyser:
Inheritance graph
[legend]
Collaboration diagram for CellMLToNektar.optimize.LinearityAnalyser:
Collaboration graph
[legend]

Public Member Functions

def analyse_for_jacobian
 
def find_linear_odes
 
def rearrange_linear_odes
 
def show
 

Static Public Attributes

tuple LINEAR_KINDS = Enum('None', 'Linear', 'Nonlinear')
 

Private Member Functions

def _get_rhs
 
def _check_expr
 
def _clone
 
def _make_apply
 
def _transfer_lut
 
def _rearrange_expr
 

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

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().

959 
960  def _check_expr(self, expr, state_var, bad_vars):
961  """The actual linearity checking function.
962 
963  Recursively determine the type of dependence expr has on
964  state_var. The presence of any members of bad_vars indicates
965  a non-linear dependency.
966 
967  Return a member of the self.LINEAR_KINDS enum.
968  """
969  kind = self.LINEAR_KINDS
970  result = None
971  if isinstance(expr, mathml_ci):
972  var = expr.variable.get_source_variable(recurse=True)
973  if var is state_var:
974  result = kind.Linear
975  elif var in bad_vars:
976  result = kind.Nonlinear
977  elif var.get_type(follow_maps=True) == VarTypes.Computed:
978  # Recurse into defining expression
979  src_var = var.get_source_variable(recurse=True)
980  src_expr = self._get_rhs(src_var.get_dependencies()[0])
981  DEBUG('find-linear-deps', "--recurse for", src_var.name,
982  "to", src_expr)
983  result = self._check_expr(src_expr, state_var, bad_vars)
984  else:
985  result = kind.None
986  # Record the kind of this variable, for later use when
987  # rearranging linear ODEs
988  var._cml_linear_kind = result
989  elif isinstance(expr, mathml_cn):
990  result = kind.None
991  elif isinstance(expr, mathml_piecewise):
992  # If any conditions have a dependence, then we're
993  # nonlinear. Otherwise, all the pieces must be the same
994  # (and that's what we are) or we're nonlinear.
995  pieces = getattr(expr, u'piece', [])
996  conds = map(lambda p: child_i(p, 2), pieces)
997  chld_exprs = map(lambda p: child_i(p, 1), pieces)
998  if hasattr(expr, u'otherwise'):
999  chld_exprs.append(child_i(expr.otherwise, 1))
1000  for cond in conds:
1001  if self._check_expr(cond, state_var, bad_vars) != kind.None:
1002  result = kind.Nonlinear
1003  break
1004  else:
1005  # Conditions all OK
1006  for e in chld_exprs:
1007  res = self._check_expr(e, state_var, bad_vars)
1008  if result is not None and res != result:
1009  # We have a difference
1010  result = kind.Nonlinear
1011  break
1012  result = res
1013  elif isinstance(expr, mathml_apply):
1014  # Behaviour depends on the operator
1015  operator = expr.operator().localName
1016  operands = expr.operands()
1017  if operator in ['plus', 'minus']:
1018  # Linear if any operand linear, and none non-linear
1019  op_kinds = map(lambda op: self._check_expr(op, state_var,
1020  bad_vars),
1021  operands)
1022  result = max(op_kinds)
1023  elif operator == 'divide':
1024  # Linear iff only numerator linear
1025  numer = operands.next()
1026  denom = operands.next()
1027  if self._check_expr(denom, state_var, bad_vars) != kind.None:
1028  result = kind.Nonlinear
1029  else:
1030  result = self._check_expr(numer, state_var, bad_vars)
1031  elif operator == 'times':
1032  # Linear iff only 1 linear operand
1033  op_kinds = map(lambda op: self._check_expr(op, state_var,
1034  bad_vars),
1035  operands)
1036  lin, nonlin = 0, 0
1037  for res in op_kinds:
1038  if res == kind.Linear: lin += 1
1039  elif res == kind.Nonlinear: nonlin += 1
1040  if nonlin > 0 or lin > 1:
1041  result = kind.Nonlinear
1042  elif lin == 1:
1043  result = kind.Linear
1044  else:
1045  result = kind.None
1046  else:
1047  # Cannot be linear; may be no dependence at all
1048  result = max(map(lambda op: self._check_expr(op, state_var,
1049  bad_vars),
1050  operands))
1051  if result == kind.Linear:
1052  result = kind.Nonlinear
1053  else:
1054  # Either a straightforward container element
1055  try:
1056  child = child_i(expr, 1)
1057  except ValueError:
1058  # Assume it's just a constant
1059  result = kind.None
1060  else:
1061  result = self._check_expr(child, state_var, bad_vars)
1062  DEBUG('find-linear-deps', "Expression", expr, "gives result", result)
1063  return result
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().

1084 
1085  def _clone(self, expr):
1086  """Properly clone a MathML sub-expression."""
1087  if isinstance(expr, mathml):
1088  clone = expr.clone_self(register=True)
1089  else:
1090  clone = mathml.clone(expr)
1091  return clone
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().

953 
954  def _get_rhs(self, expr):
955  """Return the RHS of an assignment expression."""
956  ops = expr.operands()
957  ops.next()
958  return ops.next()
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().

1093 
1094  preserve=False):
1095  """Construct a new apply expression for g or h.
1096 
1097  ghs is an iterable of (g,h) pairs for operands.
1098 
1099  i indicates whether to construct g (0) or h (1).
1100 
1101  filter_none indicates the behaviour of 0 under this operator.
1102  If True, it's an additive zero, otherwise it's a
1103  multiplicative zero.
1104  """
1105  # Find g or h operands
1106  ghs_i = map(lambda gh: gh[i], ghs)
1107  if not filter_none and None in ghs_i:
1108  # Whole expr is None
1109  new_expr = None
1110  else:
1111  # Filter out None subexprs
1112  if operator == u'minus':
1113  # Do we need to retain a unary minus?
1114  if len(ghs_i) == 1 or ghs_i[0] is None:
1115  # Original was -a or 0-a
1116  retain_unary_minus = True
1117  else:
1118  # Original was a-0 or a-b
1119  retain_unary_minus = False
1120  else:
1121  # Only retain if we're told to preserve as-is
1122  retain_unary_minus = preserve
1123  ghs_i = filter(None, ghs_i)
1124  if ghs_i:
1125  if len(ghs_i) > 1 or retain_unary_minus:
1126  new_expr = mathml_apply.create_new(
1127  self.__expr, operator, ghs_i)
1128  else:
1129  new_expr = self._clone(ghs_i[0]) # Clone may be unneeded
1130  else:
1131  new_expr = None
1132  return new_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(), and CellMLToNektar.pycml.child_i().

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

1156 
1157  def _rearrange_expr(self, expr, var):
1158  """Rearrange an expression into the form g + h*var.
1159 
1160  Performs a post-order traversal of this expression's tree,
1161  and returns a pair (g, h)
1162  """
1163 # import inspect
1164 # depth = len(inspect.stack())
1165 # print ' '*depth, "_rearrange_expr", prid(expr, True), var.name, expr
1166  gh = None
1167  if isinstance(expr, mathml_ci):
1168  # Variable
1169  ci_var = expr.variable.get_source_variable(recurse=True)
1170  if var is ci_var:
1171  gh = (None, mathml_cn.create_new(expr,
1172  u'1', u'dimensionless'))
1173  else:
1174  if ci_var._cml_linear_kind == self.LINEAR_KINDS.None:
1175  # Just put the <ci> in g, but with full name
1176  gh = (mathml_ci.create_new(expr, ci_var.fullname()), None)
1177  gh[0]._set_variable_obj(ci_var)
1178  else:
1179  # ci_var is a linear function of var, so rearrange
1180  # its definition
1181  if not hasattr(ci_var, '_cml_linear_split'):
1182  ci_defn = ci_var.get_dependencies()[0]
1183  ci_var._cml_linear_split = self._rearrange_expr(
1184  self._get_rhs(ci_defn), var)
1185  gh = ci_var._cml_linear_split
1186  elif isinstance(expr, mathml_piecewise):
1187  # The tests have to move into both components of gh:
1188  # "if C1 then (a1,b1) elif C2 then (a2,b2) else (a0,b0)"
1189  # maps to "(if C1 then a1 elif C2 then a2 else a0,
1190  # if C1 then b1 elif C2 then b2 else b0)"
1191  # Note that no test is a function of var.
1192  # First rearrange child expressions
1193  pieces = getattr(expr, u'piece', [])
1194  cases = map(lambda p: child_i(p, 1), pieces)
1195  cases_ghs = map(lambda c: self._rearrange_expr(c, var), cases)
1196  if hasattr(expr, u'otherwise'):
1197  ow_gh = self._rearrange_expr(child_i(expr.otherwise, 1), var)
1198  else:
1199  ow_gh = (None, None)
1200  # Now construct the new expression
1201  conds = map(lambda p: self._clone(child_i(p, 2)), pieces)
1202  def piecewise_branch(i):
1203  pieces_i = zip(map(lambda gh: gh[i], cases_ghs),
1204  conds)
1205  pieces_i = filter(lambda p: p[0] is not None,
1206  pieces_i) # Remove cases that are None
1207  ow = ow_gh[i]
1208  if pieces_i:
1209  new_expr = mathml_piecewise.create_new(
1210  expr, pieces_i, ow)
1211  elif ow:
1212  new_expr = ow
1213  else:
1214  new_expr = None
1215  return new_expr
1216  gh = (piecewise_branch(0), piecewise_branch(1))
1217  self._transfer_lut(expr, gh, var)
1218  elif isinstance(expr, mathml_apply):
1219  # Behaviour depends on the operator
1220  operator = expr.operator().localName
1221  operands = expr.operands()
1222  self.__expr = expr # For self._make_apply
1223  if operator in ['plus', 'minus']:
1224  # Just split the operation into each component
1225  operand_ghs = map(lambda op: self._rearrange_expr(op, var),
1226  operands)
1227  g = self._make_apply(operator, operand_ghs, 0)
1228  h = self._make_apply(operator, operand_ghs, 1)
1229  gh = (g, h)
1230  elif operator == 'divide':
1231  # (a, b) / (c, 0) = (a/c, b/c)
1232  numer = self._rearrange_expr(operands.next(), var)
1233  denom = self._rearrange_expr(operands.next(), var)
1234  assert denom[1] is None
1235  denom_g = denom[0]
1236  g = h = None
1237  if numer[0]:
1238  g = mathml_apply.create_new(expr, operator,
1239  [numer[0], denom_g])
1240  if numer[1]:
1241  if g:
1242  denom_g = self._clone(denom_g)
1243  h = mathml_apply.create_new(expr, operator,
1244  [numer[1], denom_g])
1245  gh = (g, h)
1246  elif operator == 'times':
1247  # (a1,b1)*(a2,b2) = (a1*a2, b1*a2 or a1*b2 or None)
1248  # Similarly for the nary case - at most one b_i is not None
1249  operand_ghs = map(lambda op: self._rearrange_expr(op, var),
1250  operands)
1251  g = self._make_apply(operator, operand_ghs, 0,
1252  filter_none=False)
1253  # Find non-None b_i, if any
1254  for i, ab in enumerate(operand_ghs):
1255  if ab[1] is not None:
1256  operand_ghs[i] = (ab[1], None)
1257  # Clone the a_i to avoid objects having 2 parents
1258  for j, ab in enumerate(operand_ghs):
1259  if j != i:
1260  operand_ghs[j] = (self._clone(operand_ghs[j][0]), None)
1261  h = self._make_apply(operator, operand_ghs, 0, filter_none=False)
1262  break
1263  else:
1264  h = None
1265  gh = (g, h)
1266  else:
1267  # (a, None) op (b, None) = (a op b, None)
1268  operand_ghs = map(lambda op: self._rearrange_expr(op, var),
1269  operands)
1270  g = self._make_apply(operator, operand_ghs, 0, preserve=True)
1271  gh = (g, None)
1272  self._transfer_lut(expr, gh, var)
1273  else:
1274  # Since this expression is linear, there can't be any
1275  # occurrence of var in it; all possible such cases are covered
1276  # above. So just clone it into g.
1277  gh = (self._clone(expr), None)
1278 # print ' '*depth, "Re-arranged", prid(expr, True), "to", prid(gh[0], True), ",", prid(gh[1], True)
1279  return gh
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().

1133 
1134  def _transfer_lut(self, expr, gh, var):
1135  """Transfer lookup table annotations from expr to gh.
1136 
1137  gh is a pair (g, h) s.t. expr = g + h*var.
1138 
1139  If expr can be represented by a lookup table, then the lookup
1140  variable cannot be var, since if it were, then expr would have
1141  a non-linear dependence on var. Hence h must be 0, since
1142  otherwise expr would contain a (state) variable other than the
1143  lookup variable, and hence not be a candidate for a table.
1144  Thus expr=g, so we transfer the annotations to g.
1145  """
1146  if expr.getAttributeNS(NSS['lut'], u'possible', '') != u'yes':
1147  return
1148  # Paranoia check that our reasoning is correct
1149  g, h = gh
1150  assert h is None
1151  # Transfer the annotations into g
1152  LookupTableAnalyser.copy_lut_annotations(expr, g)
1153  # Make sure g has a reference to its component, for use by code generation.
1154  g._cml_component = expr.component
1155  return
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().

918 
919  def analyse_for_jacobian(self, doc, V=None):
920  """Analyse the model for computing a symbolic Jacobian.
921 
922  Determines automatically which variables will need to be solved
923  for using Newton's method, and stores their names in
924  doc.model._cml_nonlinear_system_variables, as a list of variable
925  objects.
926 
927  Also stores doc.model._cml_linear_vars, doc.model._cml_free_var,
928  doc.model._cml_transmembrane_potential.
929  """
930  # TODO: Add error checking and tidy
931  stvs = doc.model.find_state_vars()
932  if V is None:
933  Vcname, Vvname = 'membrane', 'V'
934  V = doc.model.get_variable_by_name(Vcname, Vvname)
935  V = V.get_source_variable(recurse=True)
936  doc.model._cml_transmembrane_potential = V
937  free_var = doc.model.find_free_vars()[0]
938  lvs = self.find_linear_odes(stvs, V, free_var)
939  # Next 3 lines for benefit of rearrange_linear_odes(doc)
940  lvs.sort(key=lambda v: v.fullname())
941  doc.model._cml_linear_vars = lvs
942  doc.model._cml_free_var = free_var
943  # Store nonlinear vars in canonical order
944  nonlinear_vars = list(set(stvs) - set([V]) - set(lvs))
945  nonlinear_vars.sort(key=lambda v: v.fullname())
946  doc.model._cml_nonlinear_system_variables = nonlinear_vars
947  # Debugging
948  f = lambda var: var.fullname()
949  DEBUG('linearity-analyser', 'V=', V.fullname(), '; free var=',
950  free_var.fullname(), '; linear vars=', map(f, lvs),
951  '; nonlinear vars=', map(f, nonlinear_vars))
952  return
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().

1064 
1065  def find_linear_odes(self, state_vars, V, free_var):
1066  """Identify linear ODEs.
1067 
1068  For each ODE (except that for V), determine whether it has a
1069  linear dependence on the dependent variable, and thus can be
1070  updated directly, without using Newton's method.
1071 
1072  We also require it to not depend on any other state variable,
1073  except for V.
1074  """
1075  kind = self.LINEAR_KINDS
1076  candidates = set(state_vars) - set([V])
1077  linear_vars = []
1078  for var in candidates:
1079  ode_expr = var.get_ode_dependency(free_var)
1080  if self._check_expr(self._get_rhs(ode_expr), var,
1081  candidates - set([var])) == kind.Linear:
1082  linear_vars.append(var)
1083  return linear_vars
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().

1280 
1281  def rearrange_linear_odes(self, doc):
1282  """Rearrange the linear ODEs so they can be updated directly
1283  on solving.
1284 
1285  Each ODE du/dt = f(u, t) can be written in the form
1286  du/dt = g(t) + h(t)u.
1287  A backward Euler update step is then as simple as
1288  u_n = (u_{n-1} + g(t)dt) / (1 - h(t)dt)
1289  (assuming that the transmembrane potential has already been
1290  solved for at t_n.
1291 
1292  Stores the results in doc.model._cml_linear_update_exprs, a
1293  mapping from variable object u to pair (g, h).
1294  """
1295 
1296  odes = map(lambda v: v.get_ode_dependency(doc.model._cml_free_var),
1297  doc.model._cml_linear_vars)
1298  result = {}
1299  for var, ode in itertools.izip(doc.model._cml_linear_vars, odes):
1300  # Do the re-arrangement for this variable.
1301  # We do this by a post-order traversal of its defining ODE,
1302  # constructing a pair (g, h) for each subexpression recursively.
1303  # First, get the RHS of the ODE
1304  rhs = self._get_rhs(ode)
1305  # And traverse
1306  result[var] = self._rearrange_expr(rhs, var)
1307  # Store result in model
1308  doc.model._cml_linear_update_exprs = result
1309  return result
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.

1310 
1311  def show(self, d):
1312  """Print out a more readable report on a rearrangement,
1313  as given by self.rearrange_linear_odes."""
1314  for var, expr in d.iteritems():
1315  print var.fullname()
1316  print "G:", expr[0].xml()
1317  print "H:", expr[1].xml()
1318  print "ODE:", var.get_ode_dependency(
1319  var.model._cml_free_var).xml()
1320  print
1321 

Member Data Documentation

CellMLToNektar.optimize.LinearityAnalyser.__expr
private

Definition at line 1221 of file optimize.py.

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

tuple 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(), and CellMLToNektar.optimize.LinearityAnalyser.find_linear_odes().

Generated on Wed Mar 8 2017 12:29:22 for Nektar++ by   doxygen 1.8.8