Jacobians¶
Myokit contains two classes designed to calculate Jacobian matrices. Given a
model f
with state y
such that dy/dt = f(y)
the Jacobian is the
matrix of partial derivatives df/dy
. The JacobianCalculator
takes
a point in the state space as input and returns the corresponding Jacobian.
The JacobianTracer
does the same, but takes a full
DataLog
as input.
- class myokit.JacobianTracer(model)¶
Given a model and a simulation log, this class can revisit every logged point and generate the corresponding Jacobian.
When created, a
JacobianTracer
will generate and compile a back-end that uses automatic differentiation to evaluate the partial derivatives of the model’s right-hand side function at any given point in time.For a model:
dx/dt = f(x, t, i)
with state x, time t and inputs i, the class will generate the partial derivatives of each component of f with respect to each component of x.
Methods to calculate eigenvalues and isolate dominant eigenvalues are provided.
N.B. The partial derivatives can not be calculated for the following functions:
floor
,ceil
,abs
, quotients and remainders. If these are encountered the resulting derivatives will be yielded asNaN
. However, in many cases, these functions will only occur as part of a condition in an if statement, so theNaN
’s won’t propagate to the final result.- dominant_eigenvalues(log=None, block=None)¶
Calculates the dominant eigenvalues of the jacobian matrix for each point in time. The returned value is 1d numpy array.
The “dominant eigenvalue” is defined as the eigenvalue with the largest magnitude (
sqrt(a + bi)
). Note that the returned values may be complex.If a
DataLog
suitable forjacobians()
is given aslog
, the jacobians are calculated on the fly. To re-use a set of jacobians generated earlier, pass in theDataBlock2d
generated byjacobians()
asblock
.
- jacobians(log)¶
Calculates the Jacobian matrix for each point in
log
and returns aDataBlock2d
containing all the information from the log as well as a 2d data series containing the jacobians (stored under the key “jacobians”).The given
DataLog
must contain logged results for all states in the model and any bound variables used by the model. Bound variables whose value does not appear in the log must be unbound before creating theJacobianTracer
. Only results from one-dimensional simulations are supported.
- largest_eigenvalues(log=None, block=None)¶
Calculates the largest eigenvalues of the jacobian matrix at each point in time. The returned value is 1d numpy array.
The “largest eigenvalue” is defined as the eigenvalue with the most positive real part. Note that the returned values may be complex.
If a
DataLog
suitable forjacobians()
is given aslog
, the jacobians are calculated on the fly. To re-use a set of jacobians generated earlier, pass in theDataBlock2d
generated byjacobians()
asblock
.
- class myokit.JacobianCalculator(model)¶
Given a cell model, this class can calculate Jacobian matrices for any point in the state space.
The given model will be cloned before use. No inputs are provided by the jacobian calculator, so all default values will be used.
N.B. The partial derivatives can not be calculated for the following functions:
floor
,ceil
,abs
, quotients and remainders. If these are encountered the resulting derivatives will be yielded asNaN
. However, in many cases, these functions will only occur as part of a condition in an if statement, so theNaN
’s won’t propagate to the final result.- calculate(state)¶
Calculates both the derivatives
f
and the JacobianJ
at the given statex
and returns a tuple(f(x), J(x))
.The order of state variables must be that specified by the model (i.e. the one obtained from calling
Model.states()
).
- newton_root(x=None, accuracy=0, max_iter=50, damping=1)¶
Uses Newton’s method to search for a stable point.
An initial guess can be given as
x
, if no guess is provided the model’s initial state is used.Search is halted if the next iteration doesn’t change the result, when
max|f(x)| < accuracy
or aftermax_iter
iterations. The accuracy and maximum iterations criteria can be disabled by setting the parameters to 0.A damping factor can be applied to every step by setting a damping factor
damping
to some value between0
and1
. Withdamping=1
the full step suggested by Newton’s method is taken. With any smaller value only a fraction of the suggested step is made.Returns a tuple
(x*, f*, j*, e*)
, wherex*
is a root,f*
is the derivative vector atx*
,j*
is the Jacobian matrix at this point ande* = max|f(x*)|
.