1pyodesys
2========
3
4.. image:: http://hera.physchem.kth.se:8090/api/badges/bjodah/pyodesys/status.svg
5   :target: http://hera.physchem.kth.se:8090/bjodah/pyodesys
6   :alt: Build status on Drone
7.. image:: https://circleci.com/gh/bjodah/pyodesys.svg?style=svg
8   :target: https://circleci.com/gh/bjodah/pyodesys
9   :alt: Build status on CircleCI
10.. image:: https://img.shields.io/pypi/v/pyodesys.svg
11   :target: https://pypi.python.org/pypi/pyodesys
12   :alt: PyPI version
13.. image:: https://img.shields.io/badge/python-3.7,3.8-blue.svg
14   :target: https://www.python.org/
15   :alt: Python version
16.. image:: https://img.shields.io/pypi/l/pyodesys.svg
17   :target: https://github.com/bjodah/pyodesys/blob/master/LICENSE
18   :alt: License
19.. image:: http://img.shields.io/badge/benchmarked%20by-asv-green.svg?style=flat
20   :target: http://hera.physchem.kth.se/~pyodesys/benchmarks
21   :alt: airspeedvelocity
22.. image:: http://hera.physchem.kth.se/~pyodesys/branches/master/htmlcov/coverage.svg
23   :target: http://hera.physchem.kth.se/~pyodesys/branches/master/htmlcov
24   :alt: coverage
25.. image:: http://joss.theoj.org/papers/10.21105/joss.00490/status.svg
26   :target: https://doi.org/10.21105/joss.00490
27   :alt: Journal of Open Source Software DOI
28
29``pyodesys`` provides a straightforward way
30of numerically integrating systems of ordinary differential equations (initial value problems).
31It unifies the interface of several libraries for performing the numerical integration as well as
32several libraries for symbolic representation. It also provides a convenience class for
33representing and integrating ODE systems defined by symbolic expressions, e.g. `SymPy <http://www.sympy.org>`_
34expressions. This allows the user to write concise code and rely on ``pyodesys`` to handle the subtle differences
35between libraries.
36
37The numerical integration is performed using either:
38
39- `scipy.integrate.ode <http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html>`_
40- pygslodeiv2_
41- pyodeint_
42- pycvodes_
43
44.. _pygslodeiv2: https://github.com/bjodah/pygslodeiv2
45.. _pyodeint: https://github.com/bjodah/pyodeint
46.. _pycvodes: https://github.com/bjodah/pycvodes
47
48
49Note that implicit steppers require a user supplied callback for calculating the Jacobian.
50``pyodesys.SymbolicSys`` derives the Jacobian automatically.
51
52The symbolic representation is usually in the form of SymPy expressions, but the user may
53choose another symbolic back-end (see `sym <https://github.com/bjodah/sym>`_).
54
55When performance is of utmost importance, e.g. in model fitting where results are needed
56for a large set of initial conditions and parameters, the user may transparently
57rely on compiled native code (classes in ``pyodesys.native.native_sys`` can generate optimal C++ code).
58The major benefit is that there is no need to manually rewrite the corresponding expressions in another
59programming language.
60
61Documentation
62-------------
63Auto-generated API documentation for latest stable release is found here:
64`<https://bjodah.github.io/pyodesys/latest>`_
65(and the development version for the current master branch is found here:
66`<http://hera.physchem.kth.se/~pyodesys/branches/master/html>`_).
67
68
69Installation
70------------
71Simplest way to install pyodesys and its (optional) dependencies is to use the
72`conda package manager <http://conda.pydata.org/docs/>`_:
73
74::
75
76   $ conda install -c bjodah pyodesys pytest
77   $ python -m pytest --pyargs pyodesys
78
79Optional dependencies
80~~~~~~~~~~~~~~~~~~~~~
81If you used ``conda`` to install pyodesys_ you can skip this section.
82But if you use ``pip`` you may want to know that the default installation
83of ``pyodesys`` only requires SciPy::
84
85   $ pip install pyodesys
86   $ pytest --pyargs pyodesys -rs
87
88The above command should finish without errors but with some skipped tests.
89The reason for why some tests are skipped should be because missing optional solvers.
90To install the optional solvers you will first need to install third party libraries for
91the solvers and then their python bindings. The 3rd party requirements are as follows:
92
93- pygslodeiv2_ (requires GSL_ >=1.16)
94- pyodeint_ (requires boost_ >=1.72.0)
95- pycvodes_ (requires SUNDIALS_ >=5.3.0)
96
97.. _GSL: https://www.gnu.org/software/gsl/
98.. _boost: http://www.boost.org/
99.. _SUNDIALS: https://computation.llnl.gov/projects/sundials
100
101if you want to see what packages need to be installed on a Debian based system you may look at this
102`Dockerfile <scripts/environment/Dockerfile>`_.
103
104If you manage to install all three external libraries you may install pyodesys with the option "all"::
105
106  $ pip install pyodesys[all]
107  $ pytest --pyargs pyodesys -rs
108
109now there should be no skipped tests. If you try to install pyodesys on a machine where you do not have
110root permissions you may find the flag ``--user`` helpful when using pip. Also if there are multiple
111versions of python installed you may want to invoke python for an explicit version of python, e.g.::
112
113  $ python3.6 -m pip install --user pyodesys[all]
114
115see `setup.py <setup.py>`_ for the exact list of requirements.
116
117Using Docker
118~~~~~~~~~~~~
119If you have `Docker <https://www.docker.com>`_ installed, you may use it to host a jupyter
120notebook server::
121
122  $ ./scripts/host-jupyter-using-docker.sh . 8888
123
124the first time you run the command some dependencies will be downloaded. When the installation
125is complete there will be a link visible which you can open in your browser. You can also run
126the test suite using the same docker-image::
127
128  $ ./scripts/host-jupyter-using-docker.sh . 0
129
130there will be one skipped test (due to symengine missing in this pip installed environment) and
131quite a few instances of RuntimeWarning.
132
133Examples
134--------
135The classic van der Pol oscillator (see `examples/van_der_pol.py <examples/van_der_pol.py>`_)
136
137.. code:: python
138
139   >>> from pyodesys.symbolic import SymbolicSys
140   >>> def f(t, y, p):
141   ...     return [y[1], -y[0] + p[0]*y[1]*(1 - y[0]**2)]
142   ...
143   >>> odesys = SymbolicSys.from_callback(f, 2, 1)
144   >>> xout, yout, info = odesys.integrate(10, [1, 0], [1], integrator='odeint', nsteps=1000)
145   >>> _ = odesys.plot_result()
146   >>> import matplotlib.pyplot as plt; plt.show()  # doctest: +SKIP
147
148.. image:: https://raw.githubusercontent.com/bjodah/pyodesys/master/examples/van_der_pol.png
149
150If the expression contains transcendental functions you will need to provide a ``backend`` keyword argument:
151
152.. code:: python
153
154   >>> import math
155   >>> def f(x, y, p, backend=math):
156   ...     return [backend.exp(-p[0]*y[0])]  # analytic: y(x) := ln(kx + kc)/k
157   ...
158   >>> odesys = SymbolicSys.from_callback(f, 1, 1)
159   >>> y0, k = -1, 3
160   >>> xout, yout, info = odesys.integrate(5, [y0], [k], integrator='cvode', method='bdf')
161   >>> _ = odesys.plot_result()
162   >>> import matplotlib.pyplot as plt
163   >>> import numpy as np
164   >>> c = 1./k*math.exp(k*y0)  # integration constant
165   >>> _ = plt.plot(xout, np.log(k*(xout+c))/k, '--', linewidth=2, alpha=.5, label='analytic')
166   >>> _ = plt.legend(loc='best'); plt.show()  # doctest: +SKIP
167
168.. image:: https://raw.githubusercontent.com/bjodah/pyodesys/master/examples/lnx.png
169
170If you already have symbolic expressions created using e.g. SymPy you can create your system from those:
171
172.. code:: python
173
174   >>> import sympy as sp
175   >>> t, u, v, k  = sp.symbols('t u v k')
176   >>> dudt = v
177   >>> dvdt = -k*u  # differential equations for a harmonic oscillator
178   >>> odesys = SymbolicSys([(u, dudt), (v, dvdt)], t, [k])
179   >>> result = odesys.integrate(7, {u: 2, v: 0}, {k: 3}, integrator='gsl', method='rk8pd', atol=1e-11, rtol=1e-12)
180   >>> _ = plt.subplot(1, 2, 1)
181   >>> _ = result.plot()
182   >>> _ = plt.subplot(1, 2, 2)
183   >>> _ = plt.plot(result.xout, 2*np.cos(result.xout*3**0.5) - result.yout[:, 0])
184   >>> plt.show()  # doctest: +SKIP
185
186.. image:: https://raw.githubusercontent.com/bjodah/pyodesys/master/examples/harmonic.png
187
188You can also refer to the dependent variables by name instead of index:
189
190.. code:: python
191
192   >>> odesys = SymbolicSys.from_callback(
193   ...     lambda t, y, p: {
194   ...         'x': -p['a']*y['x'],
195   ...         'y': -p['b']*y['y'] + p['a']*y['x'],
196   ...         'z': p['b']*y['y']
197   ...     }, names='xyz', param_names='ab', dep_by_name=True, par_by_name=True)
198   ...
199   >>> t, ic, pars = [42, 43, 44], {'x': 7, 'y': 5, 'z': 3}, {'a': [11, 17, 19], 'b': 13}
200   >>> for r, a in zip(odesys.integrate(t, ic, pars, integrator='cvode'), pars['a']):
201   ...     assert np.allclose(r.named_dep('x'), 7*np.exp(-a*(r.xout - r.xout[0])))
202   ...     print('%.2f ms ' % (r.info['time_cpu']*1e3))  # doctest: +SKIP
203   ...
204   10.54 ms
205   11.55 ms
206   11.06 ms
207
208Note how we generated a list of results for each value of the parameter ``a``. When using a class
209from ``pyodesys.native.native_sys`` those integrations are run in separate threads (bag of tasks
210parallelism):
211
212.. code:: python
213
214   >>> from pyodesys.native import native_sys
215   >>> native = native_sys['cvode'].from_other(odesys)
216   >>> for r, a in zip(native.integrate(t, ic, pars), pars['a']):
217   ...     assert np.allclose(r.named_dep('x'), 7*np.exp(-a*(r.xout - r.xout[0])))
218   ...     print('%.2f ms ' % (r.info['time_cpu']*1e3))  # doctest: +SKIP
219   ...
220   0.42 ms
221   0.43 ms
222   0.42 ms
223
224For this small example we see a 20x (serial) speedup by using native code. Bigger systems often see 100x speedup.
225Since the latter is run in parallel the (wall clock) time spent waiting for the results is in practice
226further reduced by a factor equal to the number of cores of your CPU (number of threads used is set by
227the environment variable ``ANYODE_NUM_THREADS``).
228
229For further examples, see `examples/ <https://github.com/bjodah/pyodesys/tree/master/examples>`_, and rendered
230jupyter notebooks here: `<http://hera.physchem.kth.se/~pyodesys/branches/master/examples>`_
231
232Run notebooks using binder
233~~~~~~~~~~~~~~~~~~~~~~~~~~
234Using only a web-browser (and an internet connection) it is possible to explore the
235notebooks here: (by the courtesy of the people behind mybinder)
236
237.. image:: http://mybinder.org/badge.svg
238   :target: https://mybinder.org/v2/gh/bjodah/pyodesys/v0.11.6?filepath=index.ipynb
239   :alt: Binder
240
241
242Citing
243------
244If you make use of pyodesys in e.g. academic work you may cite the following peer-reviewed publication:
245
246.. image:: http://joss.theoj.org/papers/10.21105/joss.00490/status.svg
247   :target: https://doi.org/10.21105/joss.00490
248   :alt: Journal of Open Source Software DOI
249
250Depending on what underlying solver you are using you should also cite the appropriate paper
251(you can look at the list of references in the JOSS article). If you need to reference,
252in addition to the paper, a specific point version of pyodesys (for e.g. reproducibility)
253you can get per-version DOIs from the zenodo archive:
254
255.. image:: https://zenodo.org/badge/43131469.svg
256   :target: https://zenodo.org/badge/latestdoi/43131469
257   :alt: Zenodo DOI
258
259
260Licenseing
261----------
262The source code is Open Source and is released under the simplified 2-clause BSD license. See `LICENSE <LICENSE>`_ for further details.
263
264Contributing
265------------
266Contributors are welcome to suggest improvements at https://github.com/bjodah/pyodesys (see further details `here <CONTRIBUTORS.rst>`_).
267
268Author
269------
270Original author: Björn I. Dahlgren (gmail address: bjodah).
271See file `AUTHORS <AUTHORS>`_ for a list of all authors.
272