compute propagator for each block of block-diagonal system matrix

doc/index.rst

@@ -474,6 +474,12 @@ If the imaginary unit :math:`i` is found in any of the entries in :math:`\mathbf
 
 In some cases, elements of :math:`\mathbf{P}` may contain fractions that have a factor of the form :python:`param1 - param2` in their denominator. If at a later stage, the numerical value of :python:`param1` is chosen equal to that of :python:`param2`, a numerical singularity (division by zero) occurs. To avoid this issue, it is necessary to eliminate either :python:`param1` or :python:`param2` in the input, before the propagator matrix is generated. ODE-toolbox will detect conditions (in this example, :python:`param1 = param2`) under which these singularities can occur. If any conditions were found, log warning messages will be emitted during the computation of the propagator matrix. A condition is only reported if the system matrix :math:`A` is defined under that condition, ensuring that only those conditions are returned that are purely an artifact of the propagator computation.
 
+To speed up processing, the final system matrix :math:`\mathbf{A}` is rewritten as a block-diagonal matrix :math:`\mathbf{A} = \text{diag}(\mathbf{A}_1, \mathbf{A}_2, \dots, \mathbf{A}_k)`, where each of :math:`\mathbf{A}_1, \mathbf{A}_2, \dots, \mathbf{A}_k` is square. Then, the propagator matrix is computed for each individual block separately, making use of the following identity:
+
+.. math::
+
+ \exp(\mathbf{A}\cdot h) = \text{diag}(\exp(\mathbf{A}_1 \cdot h), \exp(\mathbf{A}_2 \cdot h), \dots, \exp(\mathbf{A}_k \cdot h))
+
 
 Computing the update expressions
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -572,26 +578,6 @@ The file `test_mixed_integrator_numeric.py <https://github.com/nest/ode-toolbox/
  <img src="https://raw.githubusercontent.com/nest/ode-toolbox/master/doc/fig/test_mixed_integrator_numeric.png" alt="g_in, g_in__d, g_ex, g_ex__d, V_m timeseries plots" width="620" height="451">
 
 
-Caching of results
-------------------
-
-.. admonition:: TODO
-
- Not implemented yet.
-
-Some operations on SymPy expressions can be quite slow (see the section :ref:`Working with large expressions`\ ).
-
-Even dynamical systems of moderate size can require a few minutes of processing time, in large part due to SymPy calls, and solver selection.
-
-To speed up processing, a caching mechanism analyses the final system matrix :math:`\mathbf{A}` and rewrites it as a block-diagonal matrix :math:`\mathbf{A} = \text{diag}(\mathbf{B}_1, \mathbf{B}_2, \dots, \mathbf{B}_k)`, where each of :math:`\mathbf{B}_1, \mathbf{B}_2, \dots, \mathbf{B}_k` is square.
-
-For propagators, we note that
-
-.. math::
-
- e^{\mathbf{A}t} = \text{diag}(e^{\mathbf{B}_1t}, e^{\mathbf{B}_2t}, \dots, e^{\mathbf{B}_kt})
-
-
 API documentation
 -----------------
 

odetoolbox/system_of_shapes.py

@@ -24,6 +24,7 @@
 import logging
 import numpy as np
 import scipy
+import scipy.linalg
 import scipy.sparse
 import sympy
 import sympy.matrices
@@ -41,6 +42,16 @@ def off_diagonal_is_zero(row: int, A) -> bool:
  return True
 
 
+def extract_blocks(array):
+ prev = -1
+ for i in range(len(array) - 1):
+ if array[i + 1][i] == 0 and array[i][i + 1] == 0:
+ yield array[prev + 1: i + 1, prev + 1: i + 1]
+ prev = i
+
+ yield array[prev + 1: len(array), prev + 1: len(array)]
+
+
 class PropagatorGenerationException(Exception):
  """
  Thrown in case an error occurs while generating propagators.
@@ -181,26 +192,40 @@ def get_sub_system(self, symbols):
  return SystemOfShapes(x_sub, A_sub, b_sub, c_sub, shapes_sub)
 
 
- def generate_propagator_solver(self):
- r"""
- Generate the propagator matrix and symbolic expressions for propagator-based updates; return as JSON.
+ def _generate_propagator_matrix(self, A):
+ r"""Generate the propagator matrix by matrix exponentiation.
+
+ XXX: the default custom simplification expression does not work well with sympy 1.4 here. Consider replacing sympy.simplify() with _custom_simplify_expr() if sympy 1.4 support is dropped.
  """
- #
- # generate the propagator matrix
- #
 
- P = sympy.simplify(sympy.exp(self.A_ * sympy.Symbol(Config().output_timestep_symbol))) # XXX: the default custom simplification expression does not work well with sympy 1.4 here. Consider replacing sympy.simplify() with _custom_simplify_expr() if sympy 1.4 support is dropped.
+ # naive: calculate propagators in one step
+ # P = sympy.simplify(sympy.exp(A * sympy.Symbol(Config().output_timestep_symbol)))
 
+ # optimized: be explicit about block diagonal elements; much faster!
+ blocks = list(extract_blocks(np.array(A)))
+ propagators = [sympy.simplify(sympy.exp(sympy.Matrix(block) * sympy.Symbol(Config().output_timestep_symbol))) for block in blocks]
+ P = sympy.Matrix(scipy.linalg.block_diag(*propagators))
+
+ # check the result
  if sympy.I in sympy.preorder_traversal(P):
  raise PropagatorGenerationException("The imaginary unit was found in the propagator matrix. This can happen if the dynamical system that was passed to ode-toolbox is unstable, i.e. one or more state variables will diverge to minus or positive infinity.")
 
- condition = SingularityDetection.find_singularities(P, self.A_)
+ condition = SingularityDetection.find_singularities(P, A)
  if condition:
  logging.warning("Under certain conditions, the propagator matrix is singular (contains infinities).")
  logging.warning("List of all conditions that result in a singular propagator:")
  for cond in condition:
  logging.warning("\t" + r" ∧ ".join([str(k) + " = " + str(v) for k, v in cond.items()]))
 
+ return P
+
+
+ def generate_propagator_solver(self):
+ r"""
+ Generate the propagator matrix and symbolic expressions for propagator-based updates; return as JSON.
+ """
+
+ P = self._generate_propagator_matrix(self.A_)
 
  #
  # generate symbols for each nonzero entry of the propagator matrix