Merge pull request #35 from bjodah/native

Native

.drone.yml

@@ -1,14 +1,11 @@
 build:
-  image: bjodah/bjodahimgdev:latest
+  image: bjodah/bjodahimg16:v1.2
   commands:
-    - conda config --add channels bjodah
-    - conda config --add channels conda-forge
-    - conda config --add channels symengine/label/dev
-    - conda config --set show_channel_urls True
+    - bash -c '[[ $(python setup.py --version) =~ ^[0-9]+.* ]]'
     - ./scripts/ci.sh pyodesys
-    - ./scripts/render_notebooks.sh
-    - ./scripts/generate_docs.sh
+#    - ./scripts/ci_conda.sh pyodesys
     - ./scripts/prepare_deploy.sh
+    - if grep "DO-NOT-MERGE!" -R . --exclude ".drone.yml"; then exit 1; fi
 
 deploy:
 

.gitignore

@@ -10,3 +10,4 @@ examples/*.wav
 examples/*.html
 examples/thumbs/
 *.egg-info/
+gks.svg

CHANGES.rst

@@ -1,3 +1,8 @@
+v0.7.0
+======
+- Generate (multi-threaded) native code (against pyodeint, pycvodes, pygslodeiv2)
+- OdeSys.internal_* -> OdeSys._internal
+
 v0.6.0
 ======
 - depend on package ``sym`` for symbolic backend

MANIFEST.in

@@ -1,3 +1,4 @@
+graft pyodesys/native/sources
 include AUTHORS
 include CHANGES.rst
 include LICENSE

README.rst

@@ -41,7 +41,7 @@ Documentation
 -------------
 Autogenerated API documentation for latest stable release is found here:
 `<https://pythonhosted.org/pyodesys>`_
-(and development docs for the current master branch are found here:
+(and the development version for the current master branch are found here:
 `<http://hera.physchem.kth.se/~pyodesys/branches/master/html>`_).
 
 
@@ -60,10 +60,7 @@ alternatively you may also use `pip`:
 
    $ python -m pip install --user pyodesys[all]
 
-see `setup.py`_ for optional requirements.
-
-Source distribution is available here:
-`<https://pypi.python.org/pypi/pyodesys>`_
+see `setup.py <setup.py>`_ for optional requirements.
 
 Example
 -------
@@ -76,7 +73,7 @@ The classic van der Pol oscillator (see `examples/van_der_pol.py <examples/van_d
    ...     return [y[1], -y[0] + p[0]*y[1]*(1 - y[0]**2)]
    ...
    >>> odesys = SymbolicSys.from_callback(f, 2, 1)
-   >>> xout, yout, info = odesys.integrate(10, [1, 0], [1], integrator='odeint')
+   >>> xout, yout, info = odesys.integrate(10, [1, 0], [1], integrator='odeint', nsteps=1000)
    >>> _ = odesys.plot_result()
    >>> import matplotlib.pyplot as plt; plt.show()  # doctest: +SKIP
 

conda-recipe/meta.yaml

@@ -1,9 +1,12 @@
+{% set name = "pyodesys" %}
+{% set version = "0.7.0.git" %}
+
 package:
-  name: pyodesys
-  version: 0.6.0.git
+  name: {{ name|lower }}
+  version: {{ version }}
 
 source:
-  path: ../
+  git_url: ../
 
 build:
   number: 0
@@ -14,20 +17,35 @@ requirements:
   run:
     - python
     - numpy
+    - sym
     - sympy
     - scipy
+    - appdirs
 
 test:
   imports:
-    - pyodesys
+    - {{ name }}
   requires:
     - pytest
     - scipy
     - pyodeint
     - pygslodeiv2
     - pycvodes
+    - pycompilation
+    - pycodeexport
+    - python-symengine
+    - symcxx
+    - pysym
+  commands:
+    - python -m pytest --pyargs {{ name }}
 
 about:
-  home: https://github.com/bjodah/pyodesys
-  license: BSD
-  summary: Straightforward numerical integration of ODE systems from SymPy
+  home: https://github.com/bjodah/{{ name }}
+  license: BSD 2-Clause
+  license_file: LICENSE
+  summary: 'Straightforward numerical integration of systems of ordinary differential equations'
+  doc_url: https://bjodah.github.io/{{ name }}
+
+extra:
+  recipe-maintainers:
+    - bjodah

examples/interactive-plotting.ipynb → examples/.interactive-plotting.ipynb

@@ -1,5 +1,13 @@
 {
  "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Interactive plotting\n",
+    "This notebook uses widgets allowing the user to vary an input parameter (you need to run the notebook locally for the widget to show). For performance reasons we are using [symengine](https://github.com/symengine/symengine) as the backend."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -10,23 +18,25 @@
    "source": [
     "import math\n",
     "from ipywidgets import interact\n",
-    "from pyodesys.symbolic import SymbolicSys"
+    "from pyodesys.symbolic import SymbolicSys\n",
+    "use_gr = False"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "collapsed": true
+    "collapsed": false
    },
    "outputs": [],
    "source": [
-    "#from matplotlib.pyplot import plot\n",
-    "#%matplotlib inline\n",
-    "\n",
-    "import gr\n",
-    "from gr.pygr import plot\n",
-    "gr.inline()"
+    "if use_gr:\n",
+    "    import gr\n",
+    "    from gr.pygr import plot\n",
+    "    gr.inline()    \n",
+    "else:  # use matplotlib\n",
+    "    from matplotlib.pyplot import plot\n",
+    "    %matplotlib inline"
    ]
   },
   {
@@ -70,7 +80,8 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "collapsed": false
+    "collapsed": false,
+    "scrolled": false
    },
    "outputs": [],
    "source": [
@@ -94,7 +105,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.4.3"
+   "version": "3.5.2"
   }
  },
  "nbformat": 4,

examples/transformations.ipynb

@@ -1,5 +1,13 @@
 {
  "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Solving a transformed system of first order ordinary differential equations\n",
+    "In this notebook we explore how we can use different convenient subclasses of ``SymbolicSys`` to reformulate a problem (for other stability/accuracy/performance characteristics)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -136,7 +144,7 @@
     "        ratios = odesys.stiffness()\n",
     "        plt.subplot(1, 3, 3)\n",
     "        plt.yscale('linear')\n",
-    "        plt.plot(odesys.internal_xout, ratios)\n",
+    "        plt.plot(odesys._internal[0], ratios)\n",
     "    info.pop('internal_xout')\n",
     "    info.pop('internal_yout')\n",
     "    return len(xout), info, rmsd(yout[-1, :])"
@@ -262,7 +270,7 @@
    },
    "outputs": [],
    "source": [
-    "integrate_and_plot(tsys, 'cvode', [1e-12, tend], y0aug, k, atol=1e-7, rtol=1e-7, first_step=1e-4)"
+    "integrate_and_plot(tsys, 'cvode', [1e-12, tend], y0aug, k, atol=1e-7, rtol=1e-7, first_step=1e-4, nsteps=5000)"
    ]
   },
   {
@@ -274,14 +282,14 @@
    },
    "outputs": [],
    "source": [
-    "integrate_and_plot(tsys, 'odeint', [1e-12, tend], y0aug, k, atol=1e-8, rtol=1e-8, first_step=1e-4)"
+    "integrate_and_plot(tsys, 'odeint', [1e-12, tend], y0aug, k, atol=1e-8, rtol=1e-8, first_step=1e-4, nsteps=50000)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Substantially more work was needed to solve this transformed system, also the accuracy suffered. So this transformation was not of any help for this particular transformation applied to the current problem and choice of initial conditions and length of integration. There may be situations though where it is useful."
+    "Substantially more work was needed to solve this transformed system, also the accuracy suffered. So this transformation was not of any help for this particular problem, choice of initial conditions and length of integration. There may be situations where it is useful though."
    ]
   },
   {
@@ -316,7 +324,8 @@
    },
    "outputs": [],
    "source": [
-    "integrate_and_plot(tsys2, 'cvode', [1e-12, tend], y0aug[1:], [y0aug[0], 7, 3, 2, 1], atol=1e-7, rtol=1e-7, first_step=1e-4, stiffness=True)"
+    "integrate_and_plot(tsys2, 'cvode', [1e-12, tend], y0aug[1:], [y0aug[0], 7, 3, 2, 1],\n",
+    "                   atol=1e-7, rtol=1e-7, first_step=1e-4, stiffness=True, nsteps=5000)"
    ]
   },
   {
@@ -335,7 +344,8 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Using PartiallySolvedSys"
+    "# Using PartiallySolvedSys\n",
+    "Sometimes we can solve a system paritally by integrating some of the dependent variables. The ODE system then needs to be reformulated in terms of the new anayltic function. ``PartiallySolvedSystem`` helps us with this task:"
    ]
   },
   {
@@ -366,7 +376,8 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Using fixed RK4 stepper"
+    "## Using fixed RK4 stepper\n",
+    "pyodesys provides an example class of an integrator (``RK4_example_integrator``) in ``pyodesys.integrators``. It can be used as a model integrator when designing custom steppers."
    ]
   },
   {
@@ -400,7 +411,7 @@
    },
    "outputs": [],
    "source": [
-    "integrate_and_plot(tsys, RK4_example_integartor, [1e-12, tend], y0aug, k, atol=1e-8, rtol=1e-8, first_step=1e-1)"
+    "integrate_and_plot(psys, RK4_example_integartor, [1e-12, tend], y0aug, [7, 3, 2, 1], atol=1e-8, rtol=1e-8, first_step=1e-1)"
    ]
   }
  ],
@@ -420,7 +431,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.4.3"
+   "version": "3.5.2"
   }
  },
  "nbformat": 4,

examples/van_der_pol_formulations.ipynb

@@ -25,7 +25,7 @@
     "import numpy as np\n",
     "import sympy as sp\n",
     "import matplotlib.pyplot as plt\n",
-    "from pyodesys.symbolic import SymbolicSys\n",
+    "from pyodesys.native.cvode import NativeCvodeSys as SymbolicSys\n",
     "sp.init_printing()\n",
     "%matplotlib inline\n",
     "print(sp.__version__)"
@@ -35,6 +35,8 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
+    "Note that we imported ``NativeCvodeSys`` as ``SymbolicSys``, this speed up the time of integration by more than an order of magnitude due to using compiled C++ code for our mathematical expressions.\n",
+    "\n",
     "One way to reduce the order of our second order differential equation is to formulate a system of first order ODEs, using:\n",
     "\n",
     "$$ y_1 = \\dot y_0 $$\n",
@@ -61,7 +63,7 @@
    "source": [
     "vdp1 = lambda x, y, p: [y[1], -y[0] + p[0]*y[1]*(1 - y[0]**2)]\n",
     "mu_val = 2.5\n",
-    "y0_1 = [0, 1]\n",
+    "y0_1 = [0.0, 1.0]\n",
     "y0_1, (y0_1[0], vdp1(0, y0_1, [mu_val])[0])"
    ]
   },
@@ -194,7 +196,7 @@
    "outputs": [],
    "source": [
     "tau = lambda mu: 1.6137056388801094*mu + 4.676*mu**(-1./3)\n",
-    "for mu in [20, 40, 80]:\n",
+    "for mu in [20, 40, 60]:\n",
     "    solve_and_plot(odesys1, y0_1, np.linspace(0, 5*tau(mu), 500), mu)"
    ]
   },
@@ -234,7 +236,7 @@
    },
    "outputs": [],
    "source": [
-    "solve_and_plot(odesys2, y0_2, tend, mu_val)"
+    "solve_and_plot(odesys2, y0_2, tend, mu_val, nsteps=2000)"
    ]
   },
   {
@@ -283,7 +285,7 @@
    },
    "outputs": [],
    "source": [
-    "solve_and_plot(odesys2, calc_y0_2(y0_1, mu_val), tend, mu_val)\n",
+    "solve_and_plot(odesys2, calc_y0_2(y0_1, mu_val), tend, mu_val, nsteps=2000)\n",
     "J = odesys2.get_jac()\n",
     "J"
    ]
@@ -329,7 +331,7 @@
    "outputs": [],
    "source": [
     "eigvals = np.array([(eig_cbs[0](*(tuple(yvals)+(mu_val + 0j,))),\n",
-    "                     eig_cbs[1](*(tuple(yvals)+(mu_val + 0j,)))) for yvals in odesys2.internal_yout])"
+    "                     eig_cbs[1](*(tuple(yvals)+(mu_val + 0j,)))) for yvals in odesys2._internal[1]])"
    ]
   },
   {
@@ -340,8 +342,8 @@
    },
    "outputs": [],
    "source": [
-    "plt.plot(odesys2.internal_xout, odesys2.stiffness(), label='from SVD')\n",
-    "plt.plot(odesys2.internal_xout, np.abs(eigvals[:,0])/np.abs(eigvals[:,1]), label='analytic')\n",
+    "plt.plot(odesys2._internal[0], odesys2.stiffness(), label='from SVD')\n",
+    "plt.plot(odesys2._internal[0], np.abs(eigvals[:,0])/np.abs(eigvals[:,1]), label='analytic')\n",
     "plt.legend()"
    ]
   },
@@ -350,7 +352,7 @@
    "metadata": {},
    "source": [
     "## Audio\n",
-    "Plotting is instructive from a mathematical view, but these equations were often investigated by to *listening* to audio amplified by electrical circuits modelled by the equation:"
+    "Plotting is instructive from a mathematical standpoint, but these equations were often investigated by to *listening* to audio amplified by electrical circuits modelled by the equation. So let us generate some audio."
    ]
   },
   {
@@ -502,7 +504,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.4.3"
+   "version": "3.5.2"
   }
  },
  "nbformat": 4,