DDA SciPy interface (to generic ODE solvers)¶
The dda.scipy
module allows inpython evaluation of DDA systems
as well as their solution with ODE Integrators in Python,
such as scipy.integrate.
For the usage and examples, see the main class to_scipy
.
Warning
This module exposes a solver of a DDA system which is quite different to the
cpp_exporter
. In particular,
The solver is required to be told the solution time span or final time in time units, not iteration indices.
The solver only spills out the evolved (integration) quantities and not any derived quantities. You can recompute them at any timestep, but there is currently no code helping you in order to archieve this result. This can result in confusion when you cannot query the fields you asked for (in particular in the CLI frontend).
The SciPy time integrator tries to find an optimal (and minimal) time step, yielding in a quite “rough” solution. You can turn on dense output in order to tell the SciPy solver to integrate between these time steps, yielding a more smooth output with more datapoints.

dda.scipy.
evaluate_values
(smbl, values)[source]¶ Evaluate a symbol within the context of an already evaluated values dictionary given.

class
dda.scipy.
to_scipy
(state)[source]¶ The SciPy exporter. When initializing this class with your DDA system, it will setup a function
f(y)
which can be evaluated as any right hand side in the ordinary differential equation (ODE) systemdy/dt = f(y)
. Here,y
are the evolution quantities, i.e. a vector which is composed automatically from the linearized DDA system (seeast.State.name_computing_elements()
andast.State.variable_ordering()
). Furthermore, this class prepares the initial valuesy0
for the integration.You can evaluate these quantities in any python context, i.e. with any scientific python ODE solver library. For the time being, this class provides a convenience method
solve()
which calls scipy.integrate.solve_ivp. There is no other scipy dependence in this code.Usage example:
>>> from dda.computing_elements import neg,int,mult >>> from dda import State, symbols >>> x, y, z = symbols("x, y, z") >>> dda_state = State({x: int(y, 1, 1), y: mult(2,x), z: neg(int(z, 1, 1)) }) >>> py_state = to_scipy(dda_state) >>> py_state.state # state has been linearized before processing State({'int_1': Int(z), 'x': Int(y), 'y': mult(2, x), 'z': neg(int_1)}) >>> py_state.vars.evolved # evolved variables are therefore not [x,z] but [int_1,x] ['int_1', 'x'] >>> py_state.y0 # initial values array([1, 1]) >>> py_state.dt # same timestep for all integrals (but see note below) 1 >>> py_state.rhs(py_state.y0) # evaluation of f(y) on y0 array([1, 2]) >>> y1 = py_state.rhs(py_state.y0) * py_state.dt # a single Euler integration timestep, for instance >>> y1 array([1, 2]) >>> sol = py_state.solve(10) # ODE integration with SciPy >>> sol.t # integration went from 0>10 with 17 timesteps array([ 0. , ..., 10. ]) >>> sol.y[:,1] # the first solution is ~exp(+t)>inf, the second exp(t)>0 array([2.20269685e+04, 1.94334984e08]) >>> from pylab import plot, legend, show # plotting example >>> for i,fieldname in enumerate(py_state.vars.evolved): plot(sol.t, sol.y[i], label=fieldname) >>> legend(); show()
Warning
Due to the way how widespread ODE integrators work, the perintegral step size
dt
is required to be the same for every integration which is part of the DDA system. That is, the following generates an error:>>> from dda import dda, symbols >>> a,b=symbols("a,b") >>> state = State({ a: dda.int(a, 0.2, 0), b: dda.int(b, 0.5, 0) }) >>> to_scipy(state) Traceback (most recent call last): ... ValueError: Scipy requires all timesteps to be the same, however dt_(['a', 'b']) = [0.2 0.5]
Most highlevel integrators available in scientific Python toolkits (such as scipy) assume the overall system to have a single timestep size
Δt
(which is also quite natural from a mathematical perspective). The signaturedda.int(f,dt,ic)
is thus quirky from a mathematical or numerical viewpoint. It is written in such a way because analog computing integrators have a tunable time scalek_0 ~= 1/dt
which however can also be consumed in the integrand itself:dda.int(f,1/k_0,ic) == dda.int(f/k_0,1,ic)
.Furthermore, most highlevel integrators do adaptive timestepping anyway. The finetuning of timestep sizes is something which is only paid respect to in the
cpp_exporter
module.
evaluate_const
(var)[source]¶ Translate
const(foo)
by stripping foo or some bar by looking up if it is an explicit constant. Dynamical variables are not allowed here. This is somewhat similar but different tocpp_exporter.lookup_const(var)
.

evaluate_state
(evolution_vector, copy=False)[source]¶ Recomputes the full state from the evolution state vector. Returns a dictionary with same keys as
self.state
and scalars (floats) as values.This will especially compute the aux variables, while for the evaluation variables the RHS of
dy / dt = f(y)
is computed.Note
As a user, you most likely want to call
reconstruct_state()
orrhs()
instead of this function.For optimization purpose, numerical state evaluation is always carried out on the
evaluation_default_values
member (which also hold the initial values for the firstrhs
evaluation). If you setcopy=True
, a shallow copy (which is equal to a deep copy for a dict holding floats) is returned. In external calls, you should probably always setcopy=True
.Note
The implementation of this function currently evaluates the (prepared) DDA sytem by recursive calls with the help of a variable assignment directory. This is basically a runtime compilation (JIT/VM) in pure python. Needless to say, this won’t give a great performance!
There are plenty of lowhanging fruits to provide optimized versions of this code: One could call the efficient (but still scalar) C++ implementation which
cpp_exporter
provides by methods provided by Cython. One could also map the DDA abstract syntax tree (AST) to the python one and use some unparser code to evaluate the DDA code as pure python (see for instance Python: Modify AST and write back python code).For the time being, this code remains as a pure demonstration code. Thanks to using the linearized state, there should be no troubles with call stack overflows, however cyclic dependencies may not be properly resolved and can result in infinite recursions (stack overflow).

reconstruct_state
(evolution_vector, copy=True)[source]¶ Given the evolution vector sizes, this computes the full state. That is, this function differs from
evaluate()
at all evaluation quantities where the values of the evolution vector itself are put in place.

rhs
(evolution_vector)[source]¶ The ODE right hand side in
dy / dt = f(y)
.y
is a numpy vector, andf(y)
returns a similarly sized (numpy) vector which we callrhs
here:>>> ode = to_scipy(State({ Symbol("x"): Symbol("int")(Symbol("x"),0.1,1) })) >>> y1 = ode.y0 + ode.rhs(ode.y0) * ode.dt # perform some euler integration step >>> y1 array([0.9])
Usually, you want to pass this function to some scipy integrator. See also
ft()
.


dda.scipy.
cli_scipy
()[source]¶ A Command Line Interface (CLI) for
dda.scipy
.This CLI API basically solves a DDA file (see
dda.dsl
for the syntax). This is a different approach then using thedda.cpp_exporter
Instead of code generation (and the need for a C++ compiler), this evaluates the DDA file within python. The disadvantage is that this is damned slow, the advantage is that the time integrator is much better then the selfmade one in thedda.cpp_exporter
module. And there is no need for a C++ compiler at all, all is (more or less) pure python.Invocation is like
python m dda.scipy help
anywhere from the system. Run this to explore the usage of this command line interface.The output will be CSV (file or stdout), in terms of one line per integration step (called dense solution in scipy ODESolver language).