DDA SciPy interface (to generic ODE solvers)

The dda.scipy module allows in-python 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.

Note

In order to run this code, you need, obviously, SciPy next to NumPy.

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)

Evaluate a symbol within the context of an already evaluated values dictionary given.

class dda.scipy.to_scipy(state)

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) system dy/dt = f(y). Here, y are the evolution quantities, i.e. a vector which is composed automatically from the linearized DDA system (see ast.State.name_computing_elements() and ast.State.variable_ordering()). Furthermore, this class prepares the initial values y0 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.94334984e-08])
>>> 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 per-integral 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 high-level 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 signature dda.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 scale k_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 high-level integrators do adaptive timestepping anyway. The fine-tuning of timestep sizes is something which is only paid respect to in the cpp_exporter module.

evaluate_const(var)

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 to cpp_exporter.lookup_const(var).

evaluate_state(evolution_vector, copy=False)

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() or rhs() 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 first rhs evaluation). If you set copy=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 set copy=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 run-time compilation (JIT/VM) in pure python. Needless to say, this won’t give a great performance!

There are plenty of low-hanging 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)

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)

The ODE right hand side in dy / dt = f(y). y is a numpy vector, and f(y) returns a similarly sized (numpy) vector which we call rhs 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().

rhst(t, evolution_vector)

Syntactic sugar for scipy integrators who want a signature rhst(t,y). Will just call rhs(y) instead.

solve(tfinal, **kwargs)

Basically passes all arguments to scipy.integrate.solve_ivp. See documentation for to_scipy for usage example.

Currently, it is hardcodedly tspan=[0,tfinal]. All other (keyword) arguments are passed to solve_ivp.

dda.scipy.cli_scipy()

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 the dda.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 dda.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).