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 class to_scipy
.
In order to run this code, you need, obviously, SciPy next to NumPy.

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 :mod:dsl for the syntax). This is a different approach then using the :mod: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 the :mod:cpp_exporter module.
Invocation is like
python m dda.scipy help
anywhere from the system. Run this to explore the usage of this command line interface.