diff --git a/.ci/run_ci.sh b/.ci/run_ci.sh index a9c330e..2314893 100755 --- a/.ci/run_ci.sh +++ b/.ci/run_ci.sh @@ -14,41 +14,46 @@ if [ -e /etc/profile.d/boost.sh ]; then export CPATH=$BOOST_ROOT/include fi source $(compgen -G "/opt-3/cpython-v3.*-apt-deb/bin/activate") -python -m pip install --cache-dir $CI_WORKSPACE/cache-ci/pip_cache --upgrade-strategy=eager --upgrade cython "git+https://github.com/bjodah/pycompilation@master#egg=pycompilation" +python3 -m pip install --cache-dir $CI_WORKSPACE/cache-ci/pip_cache --upgrade-strategy=eager --upgrade cython "git+https://github.com/bjodah/pycompilation@master#egg=pycompilation" # REPO_TEMP_DIR="$(mktemp -d)" # trap 'rm -rf -- "$REPO_TEMP_DIR"' EXIT # cp -ra . "$REPO_TEMP_DIR/." # cd "$REPO_TEMP_DIR" mkdir -p $HOME/.config/pip/ -echo -e "[global]\nno-cache-dir = false\ndownload-cache = $CI_WORKSPACE/cache-ci/pip_cache" >$HOME/.config/pip/pip.conf -python -m pip install mako cython -python -m pip install --no-build-isolation "git+https://github.com/bjodah/symcxx#egg=symcxx" "git+https://github.com/bjodah/pysym#egg=pysym" # unofficial backends, symengine is tested in the conda build + +echo -e "[global]\nno-cache-dir = false\ndownload-cache = $CI_WORKSPACE/cache-ci/pip_cache" >$HOME/.config/pip/pip.conf # CI_WORKSPACE or pwd? +python3 -m pip install mako cython +python3 -m pip install --no-build-isolation "git+https://github.com/bjodah/symcxx#egg=symcxx" "git+https://github.com/bjodah/pysym#egg=pysym" # unofficial backends, symengine is tested in the conda build +python3 -m pip install "git+https://github.com/bjodah/sym@jun21#egg=sym" + # (cd ./tmp/pycvodes; SUND_CFLAGS="-isystem $SUNDBASE/include $CFLAGS" SUND_LDFLAGS="-Wl,--disable-new-dtags -Wl,-rpath,$SUNDBASE/lib -L$SUNDBASE/lib $LDFLAGS" -CFLAGS="$SUND_CFLAGS $CXXFLAGS" CXXFLAGS="$SUND_CFLAGS $CXXFLAGS" LDFLAGS=$SUND_LDFLAGS python -m pip install --no-build-isolation pycvodes -CPLUS_INCLUDE_PATH="$BOOST_ROOT/include" python -m pip install --no-build-isolation "git+https://github.com/bjodah/pyodeint#egg=pyodeint" -python -m pip install --no-build-isolation "git+https://github.com/bjodah/pygslodeiv2#egg=pygslodeiv2" +git clean -xfd # -e tmp/ + +CFLAGS="$SUND_CFLAGS $CXXFLAGS" CXXFLAGS="$SUND_CFLAGS $CXXFLAGS" LDFLAGS=$SUND_LDFLAGS python3 -m pip install --no-build-isolation "git+https://github.com/bjodah/pycvodes@may21#egg=pycvodes" +python3 -m pip install --no-build-isolation "git+https://github.com/bjodah/pyodeint@sep21#egg=pyodeint" +python3 -m pip install --no-build-isolation "git+https://github.com/bjodah/pygslodeiv2#egg=pygslodeiv2" -python setup.py sdist -PKG_VERSION=$(python setup.py --version) +python3 setup.py sdist +PKG_VERSION=$(python3 setup.py --version) export PYODESYS_CVODE_FLAGS=$SUND_CFLAGS export PYODESYS_CVODE_LDFLAGS=$SUND_LDFLAGS -(cd dist/; python -m pip install "$PKG_NAME-$PKG_VERSION.tar.gz[all]"; python -m pytest -v -x --pyargs $PKG_NAME) -python -m pip uninstall --yes $PKG_NAME -python -m pip install -e .[all] -python -m pytest -xv -k test_integrate_chained_robertson pyodesys/tests/test_robertson.py +(cd dist/; python3 -m pip install "$PKG_NAME-$PKG_VERSION.tar.gz[all]"; python3 -m pytest -v -x --pyargs $PKG_NAME) +python3 -m pip uninstall --yes $PKG_NAME +python3 -m pip install -e .[all] +python3 -m pytest -xv -k test_integrate_chained_robertson pyodesys/tests/test_robertson.py export PYTHONHASHSEED=$(python3 -c "import random; print(random.randint(1,2**32-1))") -PYTHON="python -R" ./scripts/run_tests.sh --cov $PKG_NAME --cov-report html +PYTHON="python3 -R" ./scripts/run_tests.sh --cov $PKG_NAME --cov-report html ( export PYTHONPATH=$(pwd); ./scripts/render_notebooks.sh; cd $PKG_NAME/tests; jupyter nbconvert --log-level=INFO --to=html --ExecutePreprocessor.enabled=True --ExecutePreprocessor.timeout=600 *.ipynb ) ./scripts/generate_docs.sh # Test package without any 3rd party libraries that are in extras_require: -python -m pip install virtualenv -python -m virtualenv venv +python3 -m pip install virtualenv +python3 -m virtualenv venv git archive -o dist/$PKG_NAME-head.zip HEAD # test pip installable zip (symlinks break) set +u -(source ./venv/bin/activate; cd dist/; python -m pip install pytest $PKG_NAME-head.zip; python -m pytest --pyargs $PKG_NAME) +(source ./venv/bin/activate; cd dist/; python3 -m pip install pytest $PKG_NAME-head.zip; python3 -m pytest --pyargs $PKG_NAME) diff --git a/.woodpecker.yaml b/.woodpecker.yaml index 494165b..58ad5fc 100644 --- a/.woodpecker.yaml +++ b/.woodpecker.yaml @@ -26,7 +26,7 @@ steps: environment: CC: gcc-14 CXX: g++-14 - CXXFLAGS: -std=c++20 + CXXFLAGS: -std=c++20 -DBOOST_UBLAS_TYPE_CHECK=0 -DNDEBUG OMP_NUM_THREADS: 1 ANYODE_NUM_THREADS: 2 CPLUS_INCLUDE_PATH: /usr/include/suitesparse:/opt-4/boost-1.91.0-release/include @@ -48,7 +48,7 @@ steps: environment: CC: gcc-14 CXX: g++-14 - CXXFLAGS: -std=c++20 + CXXFLAGS: -std=c++20 -DBOOST_UBLAS_TYPE_CHECK=0 -DNDEBUG OMP_NUM_THREADS: 1 ANYODE_NUM_THREADS: 2 CPLUS_INCLUDE_PATH: /usr/include/suitesparse:/opt-4/boost-1.91.0-release/include diff --git a/CHANGES.rst b/CHANGES.rst index 16579dd..a99b1ba 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,3 +1,7 @@ +v0.14.6 +======= +- Bump dependencies + v0.14.5 ======= - Update setup.yp diff --git a/examples/_const_step_size.ipynb b/examples/_const_step_size.ipynb index fcc68ba..3aee5e5 100644 --- a/examples/_const_step_size.ipynb +++ b/examples/_const_step_size.ipynb @@ -3,9 +3,17 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, "outputs": [], "source": [ "from itertools import product, repeat\n", @@ -15,7 +23,7 @@ "from pyodesys import ODESys\n", "from pyodesys.integrators import (RK4_example_integrator, EulerForward_example_integrator,\n", " EulerBackward_example_integrator, Trapezoidal_example_integrator,\n", - " Midpoint_example_integrator)\n", + " Midpoint_example_integrator, BDF2FVC_example_integrator)\n", "from pyodesys.tests.bateman import bateman_full\n", "%matplotlib inline" ] @@ -23,13 +31,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "choice = RadioButtons(\n", - " options=['decay', 'sine'],\n", + " options=['sine', 'decay'],\n", " description='Case:',\n", " disabled=False\n", ")\n", @@ -39,9 +45,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "if choice.value == 'decay':\n", @@ -69,57 +73,80 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, + "outputs": [], + "source": [ + "def demo(integrator, slc=slice(None)):\n", + " xout1 = np.linspace(x0, xend, n0+1)[slc]\n", + " res1 = odesys.integrate(xout1, y0, p, first_step=xend/n0, integrator=integrator)\n", + "\n", + " fig, axes = plt.subplots(1, 2, figsize=(14, 4))\n", + " fig.suptitle(integrator.__name__)\n", + " res1.plot(ls=('-',), ax=axes[0])\n", + " axes[0].set_title(\"values\")\n", + "\n", + " xplt = np.linspace(xout1[0], xout1[-1])\n", + " yref1 = ref(res1.xout)\n", + " res1.plot(x=xplt, y=ref(xplt), ls=(':',), ax=axes[0])\n", + " res1.plot(y=res1.yout - yref1, ax=axes[1])\n", + " axes[1].set_title(\"errors\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, "outputs": [], "source": [ - "xout1 = np.linspace(x0, xend, n0 + 1)\n", - "res1 = odesys.integrate(xout1, y0, p, first_step=xend/n0, integrator=EulerForward_example_integrator)\n", - "yref1 = ref(res1.xout)\n", - "fig, axes = plt.subplots(1, 2, figsize=(14, 4))\n", - "res1.plot(ls=('-',), ax=axes[0])\n", - "xplt = np.linspace(x0, xend)\n", - "res1.plot(x=xplt, y=ref(xplt), ls=(':',), ax=axes[0])\n", - "res1.plot(y=res1.yout - yref1, ax=axes[1])" + "demo(EulerForward_example_integrator, slc=slice(0,13))" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, + "outputs": [], + "source": [ + "demo(BDF2FVC_example_integrator, slc=slice(0,13))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, "outputs": [], "source": [ "def plot_order_of_convergence(Integrator, nruns=6):\n", - " res, nruns = [], 7\n", + " res = []\n", " for i in range(nruns):\n", " n = n0 * 2**i\n", " x = np.linspace(0, 1, n)\n", " res.append(odesys.integrate(x, y0, p, integrator=Integrator))\n", "\n", - " vals = {}\n", + " vals, refs = {}, {}\n", " for ix, iy in product(range(n0), range(len(y0))):\n", " vals[ix, iy] = [res[i].yout[ix*(2**i), iy] for i in range(nruns)]\n", + " refs[ix, iy] = [ref(res[i].xout)[ix*(2**i), iy] for i in range(nruns)]\n", "\n", - " logn = np.log(n0 * 2**np.arange(nruns))\n", + " lgn = np.log10([len(r.xout) for r in res])\n", " fig, axes = plt.subplots(1, len(y0), figsize=(14, 4))\n", + " fig.suptitle(Integrator.__name__)\n", " for iy, ax in enumerate(axes):\n", " for ix in range(1, n0):\n", - " lnabserr = np.log(np.abs(np.array(vals[ix, iy]) - yref1[ix, iy]))\n", - " ax.plot(logn, lnabserr, '.-', c=(ix/n0, 0, 0))\n", - " popt = np.polyfit(logn, lnabserr, 1)\n", - " ax.plot(logn[[0,-1]], np.polyval(popt, logn[[0,-1]]), label='y=%.2f x %+.2f' % tuple(popt))\n", - " ax.legend()" + " lgabserr = np.log10(np.abs(np.array(vals[ix, iy]) - refs[ix, iy]))\n", + " ax.plot(lgn, lgabserr, '.-', c=(ix/n0, 0, 0))\n", + " popt = np.polyfit(lgn, lgabserr, 1)\n", + " ax.plot(lgn[[0,-1]], np.polyval(popt, lgn[[0,-1]]), label='err_y%d(n)=%.2f x %+.2f' % ((iy,)+tuple(popt)))\n", + " ax.legend()\n", + " ax.set_xlabel(r'$\\lg{n_\\mathrm{t}}$')\n", + " ax.set_ylabel(r'$\\lg{|e|}$')\n", + " ax.set_ylim([-13, 2])" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "plot_order_of_convergence(EulerForward_example_integrator)" @@ -128,10 +155,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "plot_order_of_convergence(EulerBackward_example_integrator)" @@ -140,9 +164,16 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, + "outputs": [], + "source": [ + "plot_order_of_convergence(BDF2FVC_example_integrator)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, "outputs": [], "source": [ "plot_order_of_convergence(Midpoint_example_integrator)" @@ -151,9 +182,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "plot_order_of_convergence(Trapezoidal_example_integrator)" @@ -162,16 +191,32 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "plot_order_of_convergence(RK4_example_integrator)" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { + "@webio": { + "lastCommId": null, + "lastKernelId": null + }, "kernelspec": { "display_name": "Python 3", "language": "python", @@ -187,16 +232,17 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.5.2" + "version": "3.8.1+" }, + "name": "_const_step_size.ipynb", "widgets": { "state": { "7d4a5a6f8d584513bb39fffa7cdc282a": { - "views": [ - { - "cell_index": 1 - } - ] + "views": { + "cell_index": [ + 1 + ] + } } }, "version": "1.2.0" diff --git a/examples/_extend_by_integration.ipynb b/examples/_extend_by_integration.ipynb index 6ed7264..3fb9a65 100644 --- a/examples/_extend_by_integration.ipynb +++ b/examples/_extend_by_integration.ipynb @@ -3,9 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -20,9 +18,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "linf, linj = get_ode_exprs()\n", @@ -34,9 +30,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "psysA = PartiallySolvedSystem.from_linear_invariants(linsys, preferred=[0], description='A ')\n", @@ -47,9 +41,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "tend, iv, pars = 1e18, [1, 0, 0], [0.04, 1e4, 3e7]\n", @@ -75,9 +67,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "linsys.autonomous_exprs, psysA.autonomous_exprs, psysC.autonomous_exprs" @@ -86,9 +76,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "linsys.autonomous_interface, psysA.autonomous_interface, psysC.autonomous_interface" @@ -97,9 +85,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "psysAr = PartiallySolvedSystem.from_linear_invariants(linsys, preferred=[0], roots=[1000*linsys.dep[0] - linsys.dep[2]])\n", @@ -109,9 +95,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "resAroot = psysAr.integrate(tend, iv, pars, return_on_root=True, **integrate_kw)\n", diff --git a/examples/_lower_bounds.ipynb b/examples/_lower_bounds.ipynb index 2a2e755..f18e5a2 100644 --- a/examples/_lower_bounds.ipynb +++ b/examples/_lower_bounds.ipynb @@ -3,9 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -20,9 +18,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "linf, linj = get_ode_exprs()\n", @@ -34,9 +30,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "psysA = PartiallySolvedSystem.from_linear_invariants(linsys, preferred=[0], description='A ')\n", @@ -47,9 +41,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "tend, iv, pars = 1e18, [1, 0, 0], [0.04, 1e4, 3e7]\n", @@ -76,9 +68,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "results = integrate_and_plot_systems(odes)" @@ -87,9 +77,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "_ = integrate_and_plot_systems([SymbolicSys.from_other(odesys, lower_bounds=[0]*odesys.ny) for odesys in odes])" @@ -98,9 +86,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [] } diff --git a/examples/_native_override_chemical_kinetics.ipynb b/examples/_native_override_chemical_kinetics.ipynb index 75886a4..9d77a9e 100644 --- a/examples/_native_override_chemical_kinetics.ipynb +++ b/examples/_native_override_chemical_kinetics.ipynb @@ -114,13 +114,11 @@ " const indextype ny = get_ny();\n", " std::vector f(ny);\n", " realtype tot=0.0;\n", - " rhs(x, y, &f[0]);\n", + " auto flag_rhs = rhs(x, y, &f[0]); if (flag_rhs != AnyODE::Status::success) { return AnyODE::Status::unrecoverable_error; }\n", " for (indextype i=0; inrev++;\n", - " return AnyODE::Status::success;\n", "\"\"\"\n", "}\n", "native_extend={\n", diff --git a/examples/_nonlinear_3body_system.ipynb b/examples/_nonlinear_3body_system.ipynb index 23b6a03..3f34823 100644 --- a/examples/_nonlinear_3body_system.ipynb +++ b/examples/_nonlinear_3body_system.ipynb @@ -3,9 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "from IPython.display import display\n", @@ -25,9 +23,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "def dydt(t, y, params=(), be=None):\n", @@ -58,9 +54,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "odesys = SymbolicSys.from_callback(dydt, 18, 3)" @@ -69,9 +63,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "y0 = [0, 0, 0, 0, 0, 1, 0, 1, 0] + [0]*9\n", @@ -82,9 +74,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "_ = res.plot(title_info=1)" @@ -93,9 +83,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "fig1 = plt.figure()\n", @@ -108,9 +96,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -122,9 +108,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "nlvls = 300\n", diff --git a/examples/_params_by_name.ipynb b/examples/_params_by_name.ipynb index 5eab775..e1f8d9e 100644 --- a/examples/_params_by_name.ipynb +++ b/examples/_params_by_name.ipynb @@ -3,9 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", @@ -16,9 +14,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "def decay(t, y, p):\n", @@ -28,9 +24,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "odesys = SymbolicSys.from_callback(decay, 1, 1, par_by_name=True, param_names=['lmb'])" @@ -39,9 +33,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "res = odesys.integrate(3, [1], {'lmb': 7, 'irrelevant_key': 'irrelevant_value'}, integrator='odeint')\n", @@ -51,9 +43,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "results = odesys.integrate(3, [1], {'lmb': [1, 5, 25]}, integrator='odeint')\n", @@ -65,9 +55,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "def decay2(t, y, p):\n", @@ -77,9 +65,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "odesys2 = SymbolicSys.from_callback(decay2, 2, 2, par_by_name=True, param_names=['foo', 'bar'])\n", @@ -90,9 +76,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "results2 = odesys2.integrate(3, [2, 1], {'foo': 2, 'bar': [1, 2, 3]}, integrator='cvode')\n", diff --git a/examples/_predefined_subexpression.ipynb b/examples/_predefined_subexpression.ipynb new file mode 100644 index 0000000..7476f52 --- /dev/null +++ b/examples/_predefined_subexpression.ipynb @@ -0,0 +1,435 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "57f784db", + "metadata": {}, + "source": [ + "## Analytic solution\n", + "Consider this system of chemical reactions:\n", + "\n", + "
\n",
+    "A -> B; k=a*f(t)\n",
+    "B -> C; k=b*f(t)\n",
+    "\n",
+    "f(t) = 1/(1s + t)\n",
+    "0 <= t\n",
+    "
\n", + "\n", + "the system of ODEs is then defined as:\n", + "$$\n", + "\\frac{\\rm{d}A}{\\rm{d}t} = -a f(t) A \\\\\n", + "\\frac{\\rm{d}B}{\\rm{d}t} = -b f(t) B + a f(t) A \\\\\n", + "\\frac{\\rm{d}C}{\\rm{d}t} = +b f(t) B\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "84232383", + "metadata": {}, + "outputs": [], + "source": [ + "from sympy import symbols, Function, Tuple, Eq, dsolve, cse" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "699f042d", + "metadata": {}, + "outputs": [], + "source": [ + "t, a, b, A0, B0, C0 = symbols('t a b A0 B0 C0', real=True, nonnegative=True)\n", + "funcs_anon = f, A, B, C = symbols('f A B C', cls=Function, nonnegative=True)\n", + "funcs_t = ft, At, Bt, Ct = [_(t) for _ in funcs_anon]\n", + "funcs_dt = dAdt, dBdt, dCdt = [_.diff(t) for _ in funcs_t[1:]]\n", + "e_ft = 1/(1+t)\n", + "r1 = a*e_ft*At\n", + "r2 = b*e_ft*Bt\n", + "rhss = e_dAdt, e_dBdt, e_dCdt = -r1, -r2 + r1, r2\n", + "odes3 = Tuple(*[Eq(_1, _2) for _1, _2 in zip(funcs_dt, rhss)])\n", + "odes3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4d4618a7", + "metadata": {}, + "outputs": [], + "source": [ + "sA = dsolve(odes3[0], func=At, ics={A(0): A0})\n", + "assert sA.subs(t, 0).rhs - A0 == 0\n", + "assert sA.rhs.diff(t) - odes3[0].rhs.subs({At: sA.rhs}) == 0\n", + "sA" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e178347c", + "metadata": {}, + "outputs": [], + "source": [ + "odes2 = odes3[1:].subs({sA.lhs: sA.rhs})\n", + "odes2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d98b1e2f", + "metadata": {}, + "outputs": [], + "source": [ + "# SymPy's dsolve (effectively?) hangs for below input:\n", + "#sB = dsolve(odes2[0], func=Bt, ics={B(0): B0})\n", + "# so we prescribe the analytic solution directly:\n", + "_c1 = B0 - a*A0/(b-a)\n", + "sB = Eq(Bt, (a*A0*(t+1)**-a / (b-a) + _c1*(t+1)**-b))#.simplify())\n", + "assert sB.subs(t, 0).rhs - B0 == 0\n", + "assert sB.rhs.diff(t).expand().factor() - odes2[0].rhs.subs({Bt: sB.rhs}).expand().factor() == 0\n", + "sB" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fdf775c9", + "metadata": {}, + "outputs": [], + "source": [ + "sB.rhs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b6edd3b7", + "metadata": {}, + "outputs": [], + "source": [ + "odes1 = odes2[1:].subs({sB.lhs: sB.rhs})\n", + "odes1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8acfbf96", + "metadata": {}, + "outputs": [], + "source": [ + "_sC = (t+1)**(-a-b)*(b*(B0*(t+1)**a + A0*(t+1)**b) - a*(t+1)**a*(A0+B0))/(a-b)\n", + "sC = Eq(Ct, (_sC + C0 - _sC.subs(t, 0)))\n", + "assert sC.subs(t, 0).rhs - C0 == 0\n", + "assert sC.rhs.diff(t).expand().factor() - odes1[0].rhs.subs({Ct: sC.rhs}).expand().factor() == 0\n", + "sC" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0693d6ca", + "metadata": {}, + "outputs": [], + "source": [ + "analytic_rhss = [sA.rhs, sB.rhs, sC.rhs]\n", + "cses, red = cse(analytic_rhss)\n", + "cses, red" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ab3c476f", + "metadata": {}, + "outputs": [], + "source": [ + "[r.subs(cses[::-1]) - ref for r, ref in zip(red, analytic_rhss)]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8a4652f8", + "metadata": {}, + "outputs": [], + "source": [ + "from pyodesys.symbolic import SymbolicSys\n", + "#help(SymbolicSys)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "466c59ed", + "metadata": {}, + "outputs": [], + "source": [ + "def mk_odesys(autonomous):\n", + " ss = tau, cA, cB, cC = symbols('tau c_A c_B c_C')\n", + " _subs = {At: cA, Bt: cB, Ct: cC, t: tau}\n", + " exprs = [(cs, e.subs(_subs)) for cs, e in zip(ss[1:], rhss)]\n", + " names = \"[A] [B] [C]\".split()\n", + " if autonomous:\n", + " exprs += [(tau,1)]\n", + " names += [\"tau\"]\n", + " indep = None\n", + " else:\n", + " indep = tau\n", + " odesys = SymbolicSys(exprs, indep, params=[a,b], names=names)\n", + " if autonomous:\n", + " assert odesys.autonomous_exprs\n", + " return odesys" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5089cd51", + "metadata": {}, + "outputs": [], + "source": [ + "odesystems = [mk_odesys(auto) for auto in [False, True]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6ec0229f", + "metadata": {}, + "outputs": [], + "source": [ + "ics = [1,0,0]\n", + "pvals = [1.3e5, 2.7]\n", + "res = []\n", + "int_kw = dict(integrator='cvode', atol=1e-12, rtol=1e-14, nsteps=5000)\n", + "for auto, odesys in enumerate(odesystems):\n", + " res.append(odesys.integrate((0,1), ics + ([0] if auto else []), pvals, **int_kw))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "86db8205", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2677b792", + "metadata": {}, + "outputs": [], + "source": [ + "from pyodesys.plotting import plot_result\n", + "#help(plot_result)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a041c101", + "metadata": {}, + "outputs": [], + "source": [ + "from sympy import lambdify\n", + "ref = lambdify([t, A0, B0, C0, a, b], analytic_rhss)\n", + "ref_vals = [ref(r.xout, *ics[:3], *pvals) for r in res]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "129dff5a", + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib import pyplot as plt\n", + "fig, all_axes = plt.subplots(2, 2, figsize=(16, 8), dpi=200)\n", + "for i, r in enumerate(res):\n", + " r.plot(ax=all_axes[i, 0]\n", + " #, info_vlines_kw=True\n", + " )\n", + " for j, rv in enumerate(ref_vals[i]):\n", + " all_axes[i, 1].plot(r.xout, rv - r.yout[:, j], label=\"error in [%s]\" % \"ABC\"[j])\n", + " all_axes[i, 1].legend()\n", + " #axes[1].set_xlim(res.xout[[0,-1]])\n", + " #axes[1].vlines(res.xout, 0, 1e-13, transform=axes[1].get_yaxis_transform(), linewidth=0.5)\n", + " for ax in all_axes[i, :]:\n", + " ax.set_xscale('log')\n", + " print({k: v for k,v in r.info.items() if not k.startswith('internal')})" + ] + }, + { + "cell_type": "markdown", + "id": "16662d12", + "metadata": {}, + "source": [ + "### Deferred below\n", + "Better CSEs but implementation not yet complete." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0901a71d", + "metadata": {}, + "outputs": [], + "source": [ + "from collections import defaultdict\n", + "\n", + "def my_cse(exprs):\n", + " cses, red = cse(exprs)\n", + " new_keys, new_values = [], []\n", + " backsubs = {}\n", + " denoms = defaultdict(list)\n", + " for x, se in cses:\n", + " _a, _b = se.as_coeff_Mul()\n", + " if _a == -1 and _b.is_Symbol:\n", + " backsubs[x] = se\n", + " continue\n", + " \n", + " #if -se in new_values:\n", + " # backsubs[x] = -new_keys[new_values.index(-se)]\n", + " # continue\n", + " \n", + " new_keys.append(x)\n", + " se = se.subs(backsubs)\n", + " new_values.append(se)\n", + " \n", + " numer, denom = se.as_numer_denom()\n", + " if not denom.is_number:\n", + " denoms[denom].append((x, numer))\n", + " taken = set()\n", + " for denom, pairs1 in denoms.items():\n", + " if denom in taken:\n", + " continue\n", + " pairs2 = denoms.get(-denom, None)\n", + " if pairs2 is None:\n", + " continue\n", + " else:\n", + " taken.add(-denom)\n", + " \n", + " def _cost(numer):\n", + " return int(numer.is_number), abs(numer) if numer.is_number else numer.count_ops()\n", + " pairs12 = (pairs1, pairs2)\n", + " scores = [sorted([(x, _cost(nur)) for x, nur in pairs], key=lambda _: _[1]) for pairs in pairs12]\n", + " x1, sc1 = scores[0][0]\n", + " x2, sc2 = scores[1][0]\n", + " use1 = sc1 < sc2\n", + " x = (x1, x2)[use1]\n", + " pairs12d = [dict(pairs) for pairs in pairs12]\n", + " print(pairs12d)\n", + " nur = pairs12d[use1][x]\n", + " kk, vv = [], []\n", + " for k, v in zip(new_keys, new_values):\n", + " kk.append(k)\n", + " if k != x and k in pairs12d[0]:\n", + " side = 0\n", + " elif k != x and k in pairs12d[1]:\n", + " side = 1\n", + " else:\n", + " vv.append(v)\n", + " continue\n", + " factor = -1 if (use1 ^ side) else 1\n", + " print(x, k, v, factor)\n", + " expr = factor*x*pairs12d[side][k]/pairs12d[use1][x]\n", + " print(expr)\n", + " vv.append(expr)\n", + " \n", + " new_keys = kk\n", + " new_values = vv\n", + " \n", + " return list(zip(new_keys, new_values)), [_.subs(backsubs) for _ in red]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5cac0790", + "metadata": {}, + "outputs": [], + "source": [ + "cses2, red2 = my_cse(analytic_rhss)\n", + "assert all(r.subs(cses2[::-1]) - ref == 0 for r, ref in zip(red2, analytic_rhss))\n", + "cses2, red2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b5dda788", + "metadata": {}, + "outputs": [], + "source": [ + "cses2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5206ab28", + "metadata": {}, + "outputs": [], + "source": [ + "numer, denom = cses2[3][1].as_numer_denom()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f1a160f0", + "metadata": {}, + "outputs": [], + "source": [ + "denom.as_numer_denom()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3892e8d4", + "metadata": {}, + "outputs": [], + "source": [ + "denom.is_number" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18b1a73b", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/_predefined_subexpression2.ipynb b/examples/_predefined_subexpression2.ipynb new file mode 100644 index 0000000..83e553f --- /dev/null +++ b/examples/_predefined_subexpression2.ipynb @@ -0,0 +1,173 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "57f784db", + "metadata": {}, + "source": [ + "## Analytic solution\n", + "Consider this system of chemical reactions:\n", + "\n", + "
\n",
+    "2 A -> B; r1=a*f(t)*A*A\n",
+    "2 B -> C; r2=b*f(t)*B*B\n",
+    "\n",
+    "f(t) = 1/(k + t)\n",
+    "0 <= t\n",
+    "
\n", + "\n", + "the system of ODEs is then defined as:\n", + "$$\n", + "\\frac{\\rm{d}A}{\\rm{d}t} = -2 a f(t) A^2 \\\\\n", + "\\frac{\\rm{d}B}{\\rm{d}t} = -2 b f(t) B^2 + a f(t) A^2 \\\\\n", + "\\frac{\\rm{d}C}{\\rm{d}t} = +b f(t) B^2\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "84232383", + "metadata": {}, + "outputs": [], + "source": [ + "from sympy import symbols, Function, Tuple, Eq, dsolve, cse" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "699f042d", + "metadata": {}, + "outputs": [], + "source": [ + "t, a, b, A0, B0, C0, k = symbols('t a b A0 B0 C0 k', real=True, nonnegative=True)\n", + "funcs_anon = f, A, B, C = symbols('f A B C', cls=Function, nonnegative=True)\n", + "funcs_t = ft, At, Bt, Ct = [_(t) for _ in funcs_anon]\n", + "funcs_dt = dAdt, dBdt, dCdt = [_.diff(t) for _ in funcs_t[1:]]\n", + "e_ft = 1/(k+t)\n", + "r1 = a*e_ft*At**2\n", + "r2 = b*e_ft*Bt**2\n", + "rhss = e_dAdt, e_dBdt, e_dCdt = -2*r1, -2*r2 + r1, r2\n", + "odes3 = Tuple(*[Eq(_1, _2) for _1, _2 in zip(funcs_dt, rhss)])\n", + "odes3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4d4618a7", + "metadata": {}, + "outputs": [], + "source": [ + "sA = dsolve(odes3[0], func=At, ics={A(0): A0})\n", + "assert sA.subs(t, 0).rhs.expand() - A0 == 0\n", + "assert sA.rhs.diff(t) - odes3[0].rhs.subs({At: sA.rhs}) == 0\n", + "sA" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7e4f8556", + "metadata": {}, + "outputs": [], + "source": [ + "sA.subs(t, 0).expand()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e178347c", + "metadata": {}, + "outputs": [], + "source": [ + "odes2 = odes3[1:].subs({sA.lhs: sA.rhs})\n", + "odes2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c405d444", + "metadata": {}, + "outputs": [], + "source": [ + "from sympy import sqrt, S, log" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d98b1e2f", + "metadata": {}, + "outputs": [], + "source": [ + "# we prescribe the analytic solution directly:\n", + "c1 = symbols('c1')\n", + "_1a = 1/(sqrt(2)*(k+t))\n", + "_1b = sqrt(a)*A0*sqrt(b)*c1\n", + "_1c = -(sqrt(2*a/b)+sqrt(2*(a-2*b)/b))\n", + "_1d = (-2*a*A0*log(k+t)+2*a*A0*log(k)+1)\n", + "_1e = -(sqrt(b)*(-sqrt(2*a/b)-sqrt(2*(a-2*b)/b)))/(2*sqrt(2)*sqrt(a))-1\n", + "_1 = _1a*_1b*_1c*_1d**_1e\n", + "_2 = A0/(k+t)*(sqrt(a*(a-2*b)) - a)\n", + "_3a = (-2*a*A0*log(k+t)+2*a*A0*log(k)+1)**(-(sqrt(b)*(-sqrt(2*a/b)-sqrt(2*(a-2*b)/b)))/sqrt(8*a))\n", + "_3b1 = -2*a*A0*log(k+t)+2*a*A0*log(k)+1\n", + "_3b2 = S(1)/2 - sqrt(a-2*b)/2/sqrt(a)\n", + "_3b = _3b1**_3b2\n", + "_3 = 2*b*(c1*(_3a + _3b))\n", + "_sB = ((k+t)*(_1+_2))/_3\n", + "_sB" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6255ddd7", + "metadata": {}, + "outputs": [], + "source": [ + "_sB.subs(t,0).expand().factor().simplify()" + ] + }, + { + "cell_type": "markdown", + "id": "ef75848a", + "metadata": {}, + "source": [ + "There's probably a typo in the analytic solution... not worth the time to investigate now." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0f686b95", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/_radiolysis.ipynb b/examples/_radiolysis.ipynb index f3b8ea7..3a36091 100644 --- a/examples/_radiolysis.ipynb +++ b/examples/_radiolysis.ipynb @@ -4,7 +4,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false, "deletable": true, "editable": true }, @@ -23,7 +22,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false, "deletable": true, "editable": true }, @@ -36,7 +34,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false, "deletable": true, "editable": true }, @@ -62,7 +59,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false, "deletable": true, "editable": true }, @@ -75,7 +71,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false, "deletable": true, "editable": true }, @@ -88,7 +83,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false, "deletable": true, "editable": true }, @@ -102,7 +96,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false, "deletable": true, "editable": true }, @@ -115,7 +108,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false, "deletable": true, "editable": true }, @@ -128,7 +120,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false, "deletable": true, "editable": true }, @@ -147,7 +138,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false, "deletable": true, "editable": true }, @@ -166,7 +156,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false, "deletable": true, "editable": true }, @@ -182,7 +171,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": true, "deletable": true, "editable": true }, diff --git a/examples/_robertson.ipynb b/examples/_robertson.ipynb index 68828b7..b0b3a42 100644 --- a/examples/_robertson.ipynb +++ b/examples/_robertson.ipynb @@ -3,9 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -20,9 +18,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "def integrate_and_plot(ax, title=None, info_vlines_kw=None, **kwargs):\n", @@ -46,9 +42,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "kw = dict(nsteps=18000, atol={'A': 1e-8, 'B': 1e-10, 'C': 1e-4}, rtol=1e-6,\n", @@ -61,9 +55,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "_ = integrate_and_plot(plt.subplot(1, 1, 1), wrapping_class=native_sys['cvode'], **kw)" @@ -72,9 +64,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "def integrate_and_plot_reduced(**kwargs):\n", @@ -95,9 +85,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "rlin = integrate_and_plot_reduced()" @@ -106,9 +94,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "kw_log = dict(atol=1e-7, logc=True, logt=True, dep_scaling=1)\n", @@ -118,9 +104,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "rloglin = integrate_and_plot_reduced(atol=1e-7, logc=True, logt=False, dep_scaling=1)" @@ -129,9 +113,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "rlinlog = integrate_and_plot_reduced(atol=1e-7, logc=False, logt=True, dep_scaling=1)" @@ -140,10 +122,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "rlog2 = integrate_and_plot_reduced(powsimp=True, **kw_log)" @@ -152,9 +131,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "import sympy as sp\n", @@ -164,9 +141,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "rlog1[0][-1].exprs" @@ -175,9 +150,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "rlog1[0][-1].get_jac()" @@ -186,10 +159,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "rlog2[0][-1].exprs" @@ -198,9 +168,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "rlog2[0][-1].get_jac()" @@ -209,9 +177,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "rlog2[0][-1]._native._written_files" @@ -220,9 +186,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "def rhs_code_in_odesys(odesys):\n", @@ -242,9 +206,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "print(''.join(rhs_code_in_odesys(rlog1[0][-1])))" @@ -253,10 +215,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "print(''.join(rhs_code_in_odesys(rlog2[0][-1])))" @@ -265,27 +224,21 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [] } diff --git a/examples/_scalar_ode.ipynb b/examples/_scalar_ode.ipynb new file mode 100644 index 0000000..1fa819e --- /dev/null +++ b/examples/_scalar_ode.ipynb @@ -0,0 +1,181 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "0e2171ed", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from pyodesys.symbolic import SymbolicSys\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d0ab26a9", + "metadata": {}, + "outputs": [], + "source": [ + "sys1 = SymbolicSys.from_callback(lambda t, y, p: [-y[0]/t], 1)\n", + "sys1.exprs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bd0b5fe4", + "metadata": {}, + "outputs": [], + "source": [ + "def vary(sys, kws, *, t_end=1e20):\n", + " res = [sys.integrate((1, t_end), [1.0], atol=1e-50, nsteps=64000#, record_order=True\n", + " , **kw) for kw in kws]\n", + " return res" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4796d001", + "metadata": {}, + "outputs": [], + "source": [ + "kws = [\n", + " dict(integrator='cvode', method='adams', rtol=1e-10),\n", + " dict(integrator='cvode', method='bdf', rtol=1e-10, first_step=1e-10),\n", + " dict(integrator='gsl', method='msbdf', rtol=1e-10),\n", + " dict(integrator='gsl', method='bsimp', rtol=1e-10),\n", + " dict(integrator='odeint', method='bs', rtol=1e-10)\n", + "]\n", + "\n", + "res1 = vary(sys1, kws)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d1323811", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_variation(res, kws, *, start_idx=0, fig_kw=None, plot_order=False):\n", + " if fig_kw is None:\n", + " fig_kw = dict(figsize=(16,6), dpi=150, sharey=True)\n", + " \n", + " if plot_order:\n", + " from mpl_toolkits.axes_grid1 import host_subplot\n", + " import mpl_toolkits.axisartist as aa\n", + " \n", + " fig = plt.figure(**fig_kw)\n", + " axes = [host_subplot(1, len(res), i+1, figure=fig, axes_class=aa.Axes) for i in range(len(res))]\n", + " else:\n", + " fig, axes = plt.subplots(1, len(res), **fig_kw)\n", + " \n", + " for ax, r, kw in zip(axes, res, kws):\n", + " #r.plot(ax=ax)\n", + " rx = 1/r.xout\n", + " relerr = np.abs(r.yout.squeeze() - rx)/rx\n", + " ax.plot(r.xout[start_idx:], np.log10(relerr[start_idx:])\n", + " , label=f\"{kw['integrator']}, {kw['method']}\"\n", + " , linewidth=0.5#, marker='x'\n", + " )\n", + " ax.set_xscale('log')\n", + " ax.set_yscale('linear')\n", + " ax.set_ylabel('log10 |relative error|')\n", + "\n", + " random_walk = np.sqrt(1.0+np.array(range(r.xout.size))) # non-systematic errors (no bias)\n", + " #ax.plot(r.xout, random_walk*kw['rtol'], label='prognosis')\n", + "\n", + " ax.legend()\n", + " ax.set_ylim([-16, -2])\n", + "\n", + " ax2 = ax.twinx()\n", + " ax2.plot(r.xout[:-1], np.diff(r.xout), color='k', linewidth=0.5\n", + " #, ls='None', marker='.', markersize=0.5, alpha=0.3\n", + " )\n", + " ax2.set_yscale('log')\n", + " ax2.set_ylim([1e-12, r.xout[-1]])\n", + " ax2.set_ylabel('step size')\n", + "\n", + " if plot_order:\n", + " ax3 = ax.twinx() \n", + " ax3.axis[\"right\"] = ax3.new_fixed_axis(loc=\"right\", offset=(60, 0))\n", + " p3, = ax3.plot(r.xout, r.info['orders'], linewidth=0.5, color='tab:orange')\n", + " ax3.axis[\"right\"].label.set_color(p3.get_color())\n", + "\n", + " ax2.axis[\"right\"].toggle(all=True)\n", + " ax3.axis[\"right\"].toggle(all=True)\n", + " \n", + " ax.set_title(f\"{r.info['n_steps']=}\\n{r.info['time_cpu']=:.3f} s\\n \")\n", + "\n", + " fig.tight_layout(pad=1.08, rect=(-.05, -.05, 1.05, 1.05))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "33d96a20", + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "plot_variation(res1, kws)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5c876e70", + "metadata": {}, + "outputs": [], + "source": [ + "sys2 = sys1.as_autonomous()\n", + "sys2.exprs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3e6a2083", + "metadata": {}, + "outputs": [], + "source": [ + "plot_variation(vary(sys2, kws), kws)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9dc8c5b1", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/_spring_error_convergence.ipynb b/examples/_spring_error_convergence.ipynb index 339833c..62cb5f7 100644 --- a/examples/_spring_error_convergence.ipynb +++ b/examples/_spring_error_convergence.ipynb @@ -3,9 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "\n", @@ -20,9 +18,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "def dydt(t, y, params, backend):\n", @@ -32,9 +28,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "A0, spring_const = 3, 5" @@ -43,9 +37,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "analytic_xmin = spring_const**-.5 * np.pi\n", @@ -57,9 +49,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "integrate_kw = dict(integrator='cvode', method='adams', nsteps=2000)" @@ -68,9 +58,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "odesys = SymbolicSys.from_callback(dydt, 2, 1)\n", @@ -83,9 +71,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "%load_ext autoreload\n", @@ -95,9 +81,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "from pyodesys.convergence import integrate_tolerance_series, fit_factory" @@ -106,9 +90,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "def plot_tolerance_series(tols, relative=False, **kwargs):\n", @@ -136,10 +118,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "plot_tolerance_series(np.logspace(-6, -4, 7))" @@ -148,9 +127,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "plot_tolerance_series(np.logspace(-6, -4, 27))" @@ -159,9 +136,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "plot_tolerance_series(np.logspace(-10, -7, 50))" @@ -170,9 +145,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "plot_tolerance_series(np.logspace(-10, -8, 6))" @@ -181,10 +154,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "plot_tolerance_series(np.logspace(-10, -8, 12))" @@ -193,9 +163,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "plot_tolerance_series(np.logspace(-10, -8, 24))" @@ -204,10 +172,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "plot_tolerance_series(np.logspace(-10, -8, 6), fit=fit_factory())" @@ -216,9 +181,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "plot_tolerance_series(np.logspace(-10, -8, 10), relative=True)" @@ -227,9 +190,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "from pyodesys.plotting import plot_result" @@ -238,10 +199,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "plot_tolerance_series(np.logspace(-2.5, -1.5, 12), fit=fit_factory(3))" @@ -250,9 +208,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [] } diff --git a/external/anyode b/external/anyode index 1e8d24f..156635a 160000 --- a/external/anyode +++ b/external/anyode @@ -1 +1 @@ -Subproject commit 1e8d24f00662cb1732cdfabd886aea1667117411 +Subproject commit 156635a11ba3a1f33027284e0eb2d0a96ca5e50a diff --git a/pyodesys/_release.py b/pyodesys/_release.py index 0321c5d..96ee97e 100644 --- a/pyodesys/_release.py +++ b/pyodesys/_release.py @@ -1 +1 @@ -__version__ = '0.14.0.dev0+git' +__version__ = '0.14.99.dev0+git' diff --git a/pyodesys/integrators.py b/pyodesys/integrators.py index 03cb2e0..1dac0a9 100644 --- a/pyodesys/integrators.py +++ b/pyodesys/integrators.py @@ -170,3 +170,67 @@ def integrate_predefined(rhs, jac, y0, xout, **kwargs): yout.append((ynew + y + euler_fw_dy)/2) x_old = x return np.array(yout), {'nfev': (len(xout)-1)} + + +class BDF2FVC_example_integrator(EulerBackward_example_integrator): + + @staticmethod + def integrate_predefined(rhs, jac, y0, tout, tol_iter=1e-12, iter_max=20, + predictor=None, + **kwargs): + if kwargs: + warnings.warn("Ignoring keyword-argumtents: %s" % ', '.join(kwargs.keys())) + yout = [y0[:], Trapezoidal_example_integrator.integrate_predefined( + rhs, jac, y0, tout[:2], **kwargs + )[0][1,:]] + ny = len(y0) + assert len(yout) == 2 and yout[1].shape == (ny,) + t_old = tout[1] + h_old = tout[1] - tout[0] + f = np.empty(ny) + J = np.empty((ny, ny)) + I = np.eye(ny) + for i, t in enumerate(tout[2:], 2): + jac(t_old, yout[-1], J) + h = t - t_old + rho = h/h_old + # https://computing.llnl.gov/projects/parallel-time-integration-multigrid/2017_BDF_Paper_v1.pdf + # Page 9, Table 2, FVC (see "Experiments on Temporal Variable Step BDF 2 Algorithms, + # Anja Katrin Denner (MSc thesis)" for derivation): + beta0 = (rho+1)/(2*rho+1) + alpha1 = -(rho+1)**2/(2*rho+1) + alpha2 = rho**2/(2*rho+1) + gamma = beta0*h + # α₀y₀ + α₁y₁ + α₂y₂ = hβ₀f(t,y₀) + # α₀ = 1, γ=hβ₀ => + # g = y₀ + α₁y₁ + α₂y₂ - γf(t, y₀) = 0 + # Find find a root of g(y₀) + # 1. y₀ ← y₁ + hf (predictor) + # 2. loop: y₀ ← y₀ - J⁻¹g (corrector) + + lu_piv = lu_factor(I - gamma*J) + #print(i, t, yout[-1]) # ,f) + if predictor is None: + ynew = yout[-1] + else: + pred = predictor.integrate_predefined(rhs, jac, yout[-1], tout[i-1:i+1], **kwargs)[0] + assert len(pred) == 2 and pred[1].shape == (ny,) + ynew = pred[1] # predictor + #print(f" rho={rho}, a1={alpha1}, a2={alpha2}, gamma={gamma}, ynew={ynew}") + norm_delta_ynew = float('inf') + for iiter in range(iter_max): + rhs(t, ynew, f) + g = ynew + alpha1*yout[-1] + alpha2*yout[-2] - gamma*f + delta_ynew = lu_solve(lu_piv, g) + ynew -= delta_ynew + norm_delta_ynew = np.sqrt(np.sum(np.square(delta_ynew))/ny) + #print(f" {iiter:4d} g={g} d={delta_ynew} |d|={norm_delta_ynew}") + if norm_delta_ynew < tol_iter: + break + else: + print("FAILURE") + return np.array(yout), dict(success=False) + yout.append(ynew) + t_old = t + h_old = h + return np.array(yout), {'nfev': (len(tout)-1), 'success': True} diff --git a/pyodesys/native/_base.py b/pyodesys/native/_base.py index f6f1d29..dee04fb 100644 --- a/pyodesys/native/_base.py +++ b/pyodesys/native/_base.py @@ -1,22 +1,29 @@ # -*- coding: utf-8 -*- from __future__ import (absolute_import, division, print_function) +from collections import defaultdict from datetime import datetime as dt from functools import reduce import logging from operator import add from pathlib import Path +import copy import os import shutil import sys import sysconfig import tempfile +import sympy +from sympy.codegen.ast import CodeBlock, Assignment, float64 import numpy as np + from ..symbolic import SymbolicSys from .. import __version__ +from .symcse.groupwise import GroupwiseCSE + try: import appdirs except ImportError: @@ -48,12 +55,57 @@ 'cplus': True, } + +def get_compile_kwargs(kwargs): + kw = copy.deepcopy(_compile_kwargs) + for k in "define include_dirs libraries flags ldflags".split(): + if k in kwargs: + if k in kw: + kw[k].extend(kwargs.pop(k)) + else: + kw[k] = kwargs.pop(k) + else: + if k not in kw: + kw[k] = [] + + if options := os.environ.get("PYODESYS_OPTIONS"): + kw['options'] = options.split(',') + return kw + + _ext_suffix = sysconfig.get_config_var('EXT_SUFFIX') _obj_suffix = '.o' # os.path.splitext(_ext_suffix)[0] + '.o' # '.obj' + + +def _r(s): + if isinstance(s, sympy.Symbol): + return s + else: + return sympy.Symbol(s, real=True) + + +class _AssignerGW: + def __init__(self, k, gw, wrap_out=lambda x: x): + self.k = k + self.gw = gw + self.n = len(gw.exprs(k)) + self.wrap_out = wrap_out + + def all(self, **kwargs): + return "\n".join(self(i, **kwargs) for i in range(self.n)) + + def expr_is_zero(self, i): + return self.gw.exprs(self.k)[i] == 0 + + def __call__(self, i, assign_to=lambda i: _r("out[%s]" % i)): + out = self.wrap_out(self.gw.exprs(self.k)[i]) + return self.gw.render(Assignment(_r(assign_to(i)), out)) + + class _NativeCodeBase(Cpp_Code): - """ Base class for generated code. + """Base class for generated code. Note kwargs ``namespace_override`` which allows the user to customize the variables used when rendering the template. @@ -68,9 +120,13 @@ class _NativeCodeBase(Cpp_Code): obj_files = ('odesys_anyode.o',) _save_temp = False - namespace_default = {'p_anon': None} + namespace_default = { + 'p_anon': None, + 'p_error_on_inf': True, + 'p_error_on_nan': True, + } namespace = { - 'p_includes': ['"odesys_anyode.hpp"'], + 'p_includes': {'"odesys_anyode.hpp"'}, 'p_support_recoverable_error': False, 'p_jacobian_set_to_zero_by_solver': False, 'p_realtype': 'double', @@ -82,16 +138,24 @@ class _NativeCodeBase(Cpp_Code): # `namespace_override` is set in init # `namespace_extend` is set in init - def __init__(self, odesys, *args, **kwargs): + def __init__(self, odesys, *args, groupwise_kw=None, assigner_kws=None, + types=None, **kwargs): + self.groupwise_kw = groupwise_kw + self.types = types + self.assigner_kws = assigner_kws or defaultdict(dict) if Cpp_Code is object: raise ModuleNotFoundError("failed to import Cpp_Code from pycodeexport") if compile_sources is None: raise ModuleNotFoundError("failed to import compile_sources from pycompilation") if odesys.nroots > 0 and not self._support_roots: raise ValueError("%s does not support nroots > 0" % self.__class__.__name__) + self.namespace_override = kwargs.pop('namespace_override', {}) self.namespace_extend = kwargs.pop('namespace_extend', {}) - self.tempdir_basename = '_pycodeexport_pyodesys_%s' % self.__class__.__name__ + self.tempdir_basename = '_pycodeexport_pyodesys_%s_%s' % ( + self.__class__.__name__, + ''.join(filter(lambda x: str.isalnum(x) or x in '_-', str(odesys.description).replace(' ', '_'))) + ) self.obj_files = self.obj_files + ('%s%s' % (self.wrapper_name, _obj_suffix),) self.so_file = '%s%s' % (self.wrapper_name, _ext_suffix) _wrapper_src0 = _native_sources_dir / ('%s.pyx' % self.wrapper_name) @@ -106,7 +170,6 @@ def __init__(self, odesys, *args, **kwargs): prebuild = {_wrapper_src: _wrapper_obj} self.build_files = self.build_files + tuple(prebuild.values()) - self.odesys = odesys for _src, _dest in prebuild.items(): if not os.path.exists(_dest): @@ -121,109 +184,79 @@ def __init__(self, odesys, *args, **kwargs): shutil.rmtree(tmpdir) if not os.path.exists(_dest): raise OSError("Failed to place prebuilt file at: %s" % _dest) - super(_NativeCodeBase, self).__init__(*args, logger=logger, **kwargs) + self.compensated_summation = kwargs.pop("compensated_summation", os.environ.get( + "PYODESYS_COMPENSATED_SUMMATION", "0") == "1") + super().__init__(*args, logger=logger, **kwargs) + + # def _ccode(self, expr, subsd): + # expr_x = expr.xreplace(subsd) + # return self.odesys.be.ccode(expr_x) def variables(self): - ny = self.odesys.ny if self.odesys.band is not None: raise NotImplementedError("Banded jacobian not yet implemented.") all_invar = tuple(self.odesys.all_invariants()) - ninvar = len(all_invar) jac = self.odesys.get_jac() nnz = self.odesys.nnz - all_exprs = self.odesys.exprs + all_invar + all_exprs = dict( + rhs=self.odesys.exprs, + invar=all_invar + ) if jac is not False and nnz < 0: jac_dfdx = list(reduce(add, jac.tolist() + self.odesys.get_dfdx().tolist())) - all_exprs += tuple(jac_dfdx) - nj = len(jac_dfdx) + all_exprs["jac_dfdt"] = jac_dfdx elif jac is not False and nnz >= 0: jac_dfdx = list(reduce(add, jac.tolist())) - all_exprs += tuple(jac_dfdx) - nj = len(jac_dfdx) - else: - nj = 0 + all_exprs["jac_dfdt"] = jac_dfdx jtimes = self.odesys.get_jtimes() if jtimes is not False: v, jtimes_exprs = jtimes - all_exprs += tuple(jtimes_exprs) - njtimes = len(jtimes_exprs) + all_exprs["jtimes"] = jtimes_exprs else: v = () jtimes_exprs = () - njtimes = 0 + + first_step = self.odesys.first_step_expr + if first_step is not None: + all_exprs["first_step"] = [first_step] + + if self.odesys.roots is not None: + all_exprs["roots"] = self.odesys.roots subsd = {k: self.odesys.be.Symbol('y[%d]' % idx) for idx, k in enumerate(self.odesys.dep)} - subsd[self.odesys.indep] = self.odesys.be.Symbol('x') - if jtimes is not False: - subsd.update({k: self.odesys.be.Symbol('v[%d]' % idx) for - idx, k in enumerate(v)}) + if self.odesys.indep is not None: + subsd[self.odesys.indep] = self.odesys.be.Symbol('x') + subsd.update({k: self.odesys.be.Symbol('m_p[%d]' % idx) for idx, k in enumerate(self.odesys.params)}) - def common_cse_symbols(): - idx = 0 - while True: - yield self.odesys.be.Symbol('m_p_cse[%d]' % idx) - idx += 1 - - def _ccode(expr): - return self.odesys.be.ccode(expr.xreplace(subsd)) - - if os.getenv('PYODESYS_NATIVE_CSE', '1') == '1': - cse_cb = self.odesys.be.cse - else: - logger.info("Not using common subexpression elimination (disabled by PYODESYS_NATIVE_CSE)") - cse_cb = lambda exprs, **kwargs: ([], exprs) - - common_cses, common_exprs = cse_cb( - all_exprs, symbols=self.odesys.be.numbered_symbols('cse_temporary'), - ignore=(self.odesys.indep,) + self.odesys.dep + v) - - common_cse_subs = {} - comm_cse_symbs = common_cse_symbols() - - for symb, subexpr in common_cses: - for expr in common_exprs: - if symb in expr.free_symbols: - common_cse_subs[symb] = next(comm_cse_symbs) - break - common_cses = [(x.xreplace(common_cse_subs), expr.xreplace(common_cse_subs)) - for x, expr in common_cses] - common_exprs = [expr.xreplace(common_cse_subs) for expr in common_exprs] - - rhs_cses, rhs_exprs = cse_cb( - common_exprs[:ny], - symbols=self.odesys.be.numbered_symbols('cse')) - - if all_invar: - invar_cses, invar_exprs = cse_cb( - common_exprs[ny:(ny + ninvar)], - symbols=self.odesys.be.numbered_symbols('cse') - ) - - if jac is not False: - jac_cses, jac_exprs = cse_cb( - common_exprs[(ny + ninvar):(ny + ninvar + nj)], - symbols=self.odesys.be.numbered_symbols('cse')) - if jtimes is not False: - jtimes_cses, jtimes_exprs = cse_cb( - common_exprs[(ny + ninvar + nj):(ny + ninvar + nj + njtimes)], - symbols=self.odesys.be.numbered_symbols('cse')) + subsd.update({k: self.odesys.be.Symbol('v[%d]' % idx) for + idx, k in enumerate(v)}) - first_step = self.odesys.first_step_expr - if first_step is not None: - first_step_cses, first_step_exprs = cse_cb( - [first_step], - symbols=self.odesys.be.numbered_symbols('cse')) + common_ignore = (() if self.odesys.indep is None else (self.odesys.indep,)) + self.odesys.dep + v + gw = GroupwiseCSE( + all_exprs, + common_cse_template="m_cse[{}]", + common_ignore=common_ignore, + subsd=subsd, + # Transformer=Transformer, + # transformer_kw=transformer_kw, + **(self.groupwise_kw or {}) + ) - if self.odesys.roots is not None: - roots_cses, roots_exprs = cse_cb( - self.odesys.roots, - symbols=self.odesys.be.numbered_symbols('cse')) + def not_arr(s): + return '[' not in s.name + types = self.types or defaultdict(lambda: (lambda lhs, rhs: (float64, rhs))) + def _cses(k): + return CodeBlock(*gw.statements(k, declare=not_arr, type_=types[k])) + cses = {k: gw.render(_cses(k)) for k in gw.keys} + n_common_cses = gw.n_remapped + common_cses = gw.render(CodeBlock(*gw.common_statements(declare=not_arr, type_=types[None]))) + assigners = {k: _AssignerGW(k, gw, **self.assigner_kws[k]) for k in gw.keys} ns = dict( _message_for_rendered=[ @@ -231,52 +264,58 @@ def _ccode(expr): "This file was generated using pyodesys-%s at %s" % ( __version__, dt.now().isoformat()) ], - p_odesys=self.odesys, p_common={ - 'cses': [(symb.name, _ccode(expr)) for symb, expr in common_cses], - 'nsubs': len(common_cse_subs) + 'cses': common_cses, + 'n_cses': n_common_cses }, + p_odesys=self.odesys, p_rhs={ - 'cses': [(symb.name, _ccode(expr)) for symb, expr in rhs_cses], - 'exprs': list(map(_ccode, rhs_exprs)) + 'cses': cses["rhs"], + 'assign': assigners["rhs"] }, p_jtimes=None if jtimes is False else{ - 'cses': [(symb.name, _ccode(expr)) for symb, expr in jtimes_cses], - 'exprs': list(map(_ccode, jtimes_exprs)) + 'cses': cses["jtimes"], + 'assign': assigners["jtimes"] }, p_jac_dense=None if jac is False or nnz >= 0 else { - 'cses': [(symb.name, _ccode(expr)) for symb, expr in jac_cses], - 'exprs': {(idx//ny, idx % ny): _ccode(expr) - for idx, expr in enumerate(jac_exprs[:ny*ny])}, - 'dfdt_exprs': list(map(_ccode, jac_exprs[ny*ny:])) + 'cses': cses["jac_dfdt"], + 'assign': assigners["jac_dfdt"] }, p_jac_sparse=None if jac is False or nnz < 0 else { - 'cses': [(symb.name, _ccode(expr)) for symb, expr in jac_cses], - 'exprs': list(map(_ccode, jac_exprs[:nj])), + 'cses': cses["jac_dfdt"], + 'assign': assigners["jac_dfdt"], 'colptrs': self.odesys._colptrs, 'rowvals': self.odesys._rowvals }, p_first_step=None if first_step is None else { - 'cses': first_step_cses, - 'expr': _ccode(first_step_exprs[0]), + 'cses': cses["first_step"], + 'assign': assigners["first_step"] }, p_roots=None if self.odesys.roots is None else { - 'cses': [(symb.name, _ccode(expr)) for symb, expr in roots_cses], - 'exprs': list(map(_ccode, roots_exprs)) + 'cses': cses["roots"], + 'assign': assigners["roots"] }, p_invariants=None if all_invar == () else { - 'cses': [(symb.name, _ccode(expr)) for symb, expr in invar_cses], - 'exprs': list(map(_ccode, invar_exprs)) + 'cses': cses["invar"], + 'assign': assigners["invar"], + 'n_invar': len(all_invar) }, p_nroots=self.odesys.nroots, p_constructor=[], - p_get_dx_max=False + p_get_dx_max=False, + p_info_comment_codegen=f"{self.groupwise_kw=}", + p_compensated_summation=self.compensated_summation ) ns.update(self.namespace_default) ns.update(self.namespace) ns.update(self.namespace_override) for k, v in self.namespace_extend.items(): - ns[k].extend(v) + if isinstance(ns[k], list): + ns[k] = ns[k] + v + elif isinstance(ns[k], set): + ns[k] = ns[k] | v + else: + raise NotImplementedError(f"Cannot extend {k} of type {type(k)}") return ns @@ -286,17 +325,14 @@ class _NativeSysBase(SymbolicSys): _NativeCode = None _native_name = None - def __init__(self, *args, **kwargs): - namespace_override = kwargs.pop('namespace_override', {}) - namespace_extend = kwargs.pop('namespace_extend', {}) - save_temp = kwargs.pop('save_temp', False) + def __init__(self, *args, native_code_kw=None, **kwargs): if 'init_indep' not in kwargs: # we need to trigger append_iv for when invariants are used kwargs['init_indep'] = True kwargs['init_dep'] = True super(_NativeSysBase, self).__init__(*args, **kwargs) - self._native = self._NativeCode(self, save_temp=save_temp, - namespace_override=namespace_override, - namespace_extend=namespace_extend) + self._native = self._NativeCode( + self, + **(native_code_kw or {})) def integrate(self, *args, **kwargs): integrator = kwargs.pop('integrator', 'native') @@ -305,7 +341,13 @@ def integrate(self, *args, **kwargs): else: kwargs['integrator'] = 'native' - return super(_NativeSysBase, self).integrate(*args, **kwargs) + return super().integrate(*args, **kwargs) + + def rhs(self, intern_t, intern_y, intern_p): + return self._native.mod.rhs(intern_t, intern_y, intern_p) + + def jac(self, intern_t, intern_y, intern_p): + return self._native.mod.dense_jac_cmaj(intern_t, intern_y, intern_p) def _integrate_native(self, intern_x, intern_y0, intern_p, force_predefined=False, atol=1e-8, rtol=1e-8, nsteps=500, first_step=0.0, **kwargs): diff --git a/pyodesys/native/cvode.py b/pyodesys/native/cvode.py index 8eae211..402590e 100644 --- a/pyodesys/native/cvode.py +++ b/pyodesys/native/cvode.py @@ -6,7 +6,7 @@ import sys from ..util import import_ -from ._base import _NativeCodeBase, _NativeSysBase, _compile_kwargs +from ._base import _NativeCodeBase, _NativeSysBase, get_compile_kwargs get_include, config, _libs = import_("pycvodes", "get_include", "config", "_libs") @@ -26,7 +26,7 @@ class NativeCvodeCode(_NativeCodeBase): _indextype = '#error "indextype_failed-to-import-pycvodes-or-too-old-version"' namespace = { - 'p_includes': ['"odesys_anyode_iterative.hpp"'], + 'p_includes': {'"odesys_anyode_iterative.hpp"'}, 'p_support_recoverable_error': True, 'p_jacobian_set_to_zero_by_solver': True, 'p_baseclass': 'OdeSysIterativeBase', @@ -36,16 +36,17 @@ class NativeCvodeCode(_NativeCodeBase): _support_roots = True def __init__(self, *args, **kwargs): - self.compile_kwargs = copy.deepcopy(_compile_kwargs) - self.compile_kwargs['define'] = ['PYCVODES_NO_KLU={}'.format("0" if config.get('KLU', True) else "1"), - 'PYCVODES_NO_LAPACK={}'.format("0" if config.get('LAPACK', True) else "1"), - 'ANYODE_NO_LAPACK={}'.format("0" if config.get('LAPACK', True) else "1")] + self.compile_kwargs = get_compile_kwargs(kwargs) + self.compile_kwargs['define'].extend([ + 'PYCVODES_NO_KLU={}'.format("0" if config.get('KLU', True) else "1"), + 'PYCVODES_NO_LAPACK={}'.format("0" if config.get('LAPACK', True) else "1"), + 'ANYODE_NO_LAPACK={}'.format("0" if config.get('LAPACK', True) else "1") + ]) self.compile_kwargs['include_dirs'].append(get_include()) - self.compile_kwargs['libraries'].extend(_libs.get_libs().split(',')) - self.compile_kwargs['libraries'].extend([l for l in os.environ.get( + self.compile_kwargs['libraries'].extend(_libs.get_libs().split(',') + [l for l in os.environ.get( 'PYODESYS_LAPACK', "lapack,blas" if config["LAPACK"] else "").split(",") if l != ""]) - self.compile_kwargs['flags'] = [f for f in os.environ.get("PYODESYS_CVODE_FLAGS", "").split() if f] - self.compile_kwargs['ldflags'] = [f for f in os.environ.get("PYODESYS_CVODE_LDFLAGS", "").split() if f] + self.compile_kwargs['flags'].extend([f for f in os.environ.get("PYODESYS_CVODE_FLAGS", "").split() if f]) + self.compile_kwargs['ldflags'].extend([f for f in os.environ.get("PYODESYS_CVODE_LDFLAGS", "").split() if f]) super(NativeCvodeCode, self).__init__(*args, **kwargs) diff --git a/pyodesys/native/gsl.py b/pyodesys/native/gsl.py index 62bbbb6..94e2a3d 100644 --- a/pyodesys/native/gsl.py +++ b/pyodesys/native/gsl.py @@ -4,7 +4,7 @@ import copy import os from ..util import import_ -from ._base import _NativeCodeBase, _NativeSysBase, _compile_kwargs +from ._base import _NativeCodeBase, _NativeSysBase, get_compile_kwargs _config, get_include = import_('pygslodeiv2', '_config', 'get_include') @@ -14,7 +14,7 @@ class NativeGSLCode(_NativeCodeBase): wrapper_name = '_gsl_wrapper' def __init__(self, *args, **kwargs): - self.compile_kwargs = copy.deepcopy(_compile_kwargs) + self.compile_kwargs = get_compile_kwargs(kwargs) self.compile_kwargs['include_dirs'].append(get_include()) self.compile_kwargs['libraries'].extend(_config.env['GSL_LIBS'].split(',')) self.compile_kwargs['libraries'].extend(os.environ.get('PYODESYS_BLAS', _config.env['BLAS']).split(',')) diff --git a/pyodesys/native/odeint.py b/pyodesys/native/odeint.py index a122f9b..9dd41ee 100644 --- a/pyodesys/native/odeint.py +++ b/pyodesys/native/odeint.py @@ -4,7 +4,7 @@ import copy from ..util import import_ -from ._base import _NativeCodeBase, _NativeSysBase, _compile_kwargs +from ._base import _NativeCodeBase, _NativeSysBase, get_compile_kwargs pyodeint = import_('pyodeint') @@ -13,7 +13,7 @@ class NativeOdeintCode(_NativeCodeBase): wrapper_name = '_odeint_wrapper' def __init__(self, *args, **kwargs): - self.compile_kwargs = copy.deepcopy(_compile_kwargs) + self.compile_kwargs = get_compile_kwargs(kwargs) self.compile_kwargs['include_dirs'].append(pyodeint.get_include()) self.compile_kwargs['libraries'].extend(['m']) super(NativeOdeintCode, self).__init__(*args, **kwargs) diff --git a/pyodesys/native/sources/_cvode_wrapper.pyx b/pyodesys/native/sources/_cvode_wrapper.pyx index e261c81..b45602d 100644 --- a/pyodesys/native/sources/_cvode_wrapper.pyx +++ b/pyodesys/native/sources/_cvode_wrapper.pyx @@ -75,6 +75,57 @@ def _as_dict(unordered_map[string, int] nfo, return dct +def rhs(double t, floating [:] y, floating [:] p): + cdef: + #cnp.ndarray[realtype, ndim=1, mode='c'] t_arr = np.asarray(t, dtype=dtype) + cnp.ndarray[realtype, ndim=1, mode='c'] y_arr = np.asarray(y, dtype=dtype) + cnp.ndarray[realtype, ndim=1, mode='c'] fout = np.zeros(y.size, dtype=dtype) + cnp.ndarray[realtype, ndim=1, mode='c'] params_arr = np.asarray(p, dtype=dtype) + realtype rtol = 1e-9 + vector[realtype] atol_vec + vector[realtype] special_settings_vec + realtype get_dx_max_factor = 1.0 + bool error_outside_bounds = True + realtype max_invariant_violation = 1.0 + atol_vec.resize(y.size, 1.0) + #assert t_arr.size == 1 + cdef CvodesOdeSys * odesys = new CvodesOdeSys( + #(NULL) if p.shape[0] == 0 else + params_arr.data, + atol_vec, + rtol, + get_dx_max_factor, error_outside_bounds, + max_invariant_violation, special_settings_vec) + odesys.rhs(t, y_arr.data, &fout[0]) + del odesys + return fout + +def dense_jac_cmaj(double t, floating [:] y, floating [:] p): + cdef: + #cnp.ndarray[realtype, ndim=1, mode='c'] t_arr = np.asarray(t, dtype=dtype) + cnp.ndarray[realtype, ndim=1, mode='c'] y_arr = np.asarray(y, dtype=dtype) + cnp.ndarray[realtype, ndim=2, mode='fortran'] j_out = np.zeros((y.size, y.size), dtype=dtype, order='F') + cnp.ndarray[realtype, ndim=1, mode='c'] params_arr = np.asarray(p, dtype=dtype) + realtype rtol = 1e-9 + vector[realtype] atol_vec + vector[realtype] special_settings_vec + realtype get_dx_max_factor = 1.0 + bool error_outside_bounds = True + realtype max_invariant_violation = 1.0 + atol_vec.resize(y.size, 1.0) + #assert t_arr.size == 1 + cdef CvodesOdeSys * odesys = new CvodesOdeSys( + #(NULL) if p.shape[0] == 0 else + params_arr.data, + atol_vec, + rtol, + get_dx_max_factor, error_outside_bounds, + max_invariant_violation, special_settings_vec) + odesys.dense_jac_cmaj(t, y_arr.data, NULL, &j_out[0,0], y.size) + del odesys + return j_out + + def integrate_adaptive(floating [:, ::1] y0, floating [::1] x0, floating [::1] xend, diff --git a/pyodesys/native/sources/_odeint_wrapper.pyx b/pyodesys/native/sources/_odeint_wrapper.pyx index b93dfa0..0fec045 100644 --- a/pyodesys/native/sources/_odeint_wrapper.pyx +++ b/pyodesys/native/sources/_odeint_wrapper.pyx @@ -87,7 +87,7 @@ def integrate_adaptive(cnp.ndarray[cnp.float64_t, ndim=2, mode='c'] y0, [atol_scalar], rtol, 1.0, False, 0.0, special_settings)) result = multi_adaptive[OdeintOdeSys]( - systems, atol, rtol, styp_from_name(_styp), y0.data, + systems, atol_scalar, rtol, styp_from_name(_styp), y0.data, x0.data, xend.data, mxsteps, &_dx0[0], &_dx_max[0], autorestart, return_on_error) diff --git a/pyodesys/native/sources/anyode.pxd b/pyodesys/native/sources/anyode.pxd index 12e3be6..b43ddfc 100644 --- a/pyodesys/native/sources/anyode.pxd +++ b/pyodesys/native/sources/anyode.pxd @@ -6,13 +6,20 @@ from libcpp.vector cimport vector from libcpp cimport bool cdef extern from "anyode/anyode.hpp" namespace "AnyODE": - cdef cppclass OdeSysBase[Real_t, Index_t]: - int nfev, njev, njvev - bool use_get_dx_max + cdef cppclass OdeSysBase[Real_t, Index_t]: + int nfev, njev, njvev + bool use_get_dx_max -cdef extern from "anyode/anyode.hpp" namespace "AnyODE": - cdef cppclass Info: + cdef cppclass Info: unordered_map[string, int] nfo_int unordered_map[string, double] nfo_dbl unordered_map[string, vector[double]] nfo_vecdbl unordered_map[string, vector[int]] nfo_vecint + + cdef cppclass Status: + pass # Status is an enum class + +cdef extern from "anyode/anyode.hpp" namespace "AnyODE::Status": + cdef Status success + cdef Status recoverable_error + cdef Status unrecoverable_error diff --git a/pyodesys/native/sources/odesys_anyode.hpp b/pyodesys/native/sources/odesys_anyode.hpp index 08a5def..521ca8e 100644 --- a/pyodesys/native/sources/odesys_anyode.hpp +++ b/pyodesys/native/sources/odesys_anyode.hpp @@ -5,10 +5,11 @@ namespace odesys_anyode { template struct OdeSys : public AnyODE::OdeSysBase { std::vector m_p; - std::vector m_p_cse; + std::vector m_cse; std::vector m_atol; std::vector m_upper_bounds; std::vector m_lower_bounds; + std::vector m_invar; std::vector m_invar0; Real_t m_rtol; Real_t m_get_dx_max_factor; diff --git a/pyodesys/native/sources/odesys_anyode_iterative.hpp b/pyodesys/native/sources/odesys_anyode_iterative.hpp index 705d7b8..2103ae4 100644 --- a/pyodesys/native/sources/odesys_anyode_iterative.hpp +++ b/pyodesys/native/sources/odesys_anyode_iterative.hpp @@ -5,10 +5,11 @@ namespace odesys_anyode { template struct OdeSys : public AnyODE::OdeSysIterativeBase { std::vector m_p; - std::vector m_p_cse; + std::vector m_cse; std::vector m_atol; std::vector m_upper_bounds; std::vector m_lower_bounds; + std::vector m_invar; std::vector m_invar0; Real_t m_rtol; Real_t m_get_dx_max_factor; diff --git a/pyodesys/native/sources/odesys_anyode_iterative.pxd b/pyodesys/native/sources/odesys_anyode_iterative.pxd index 088410a..c334f2a 100644 --- a/pyodesys/native/sources/odesys_anyode_iterative.pxd +++ b/pyodesys/native/sources/odesys_anyode_iterative.pxd @@ -6,7 +6,7 @@ from libcpp.string cimport string from libcpp.unordered_map cimport unordered_map from libcpp.vector cimport vector -from anyode cimport Info +from anyode cimport Info, Status cdef extern from "odesys_anyode_iterative.hpp" namespace "odesys_anyode": cdef cppclass OdeSys[Real_t, Index_t]: @@ -18,3 +18,9 @@ cdef extern from "odesys_anyode_iterative.hpp" namespace "odesys_anyode": bool record_jac_xvals bool record_order bool record_fpe + Status rhs(Real_t t, const Real_t * const y, Real_t * const f) except + + Status dense_jac_cmaj(Real_t t, + const Real_t * const y, + const Real_t * const fy, + Real_t * const jac, + long int ldim) except + diff --git a/pyodesys/native/sources/odesys_anyode_template.cpp b/pyodesys/native/sources/odesys_anyode_template.cpp index 8ddcab5..e7a770a 100644 --- a/pyodesys/native/sources/odesys_anyode_template.cpp +++ b/pyodesys/native/sources/odesys_anyode_template.cpp @@ -3,9 +3,25 @@ <%doc> This is file is a mako template for a C++ source file defining the ODE system. +<%! +import sympy + +%> +<% +def _inf_or_nan(s): + isinf_isnan = {'isinf(%s)': p_error_on_inf, 'isnan(%s)': p_error_on_nan} + return ' || '.join([func % s for func, pred in isinf_isnan.items()]) + +def _err_code(): + return 'AnyODE::Status::recoverable_error' if p_support_recoverable_error else 'AnyODE::Status::unrecoverable_error' + +assert p_error_on_nan in (False, True) +assert p_error_on_inf in (False, True) +%> // User provided system description: ${p_odesys.description} // Names of dependent variables: ${p_odesys.names} // Names of parameters: ${p_odesys.param_names} +// Code-generation details: ${p_info_comment_codegen} #include #include @@ -16,8 +32,15 @@ This is file is a mako template for a C++ source file defining the ODE system. #include ${inc} %endfor +%if p_compensated_summation: +#include "summation_cxx/ranged.hpp" +%endif + +// sympy wraps real around symbols which actually are real: +#define re(x) (x) + namespace { // anonymous namespace for user-defined helper functions - std::vector p_odesys_names ${'' if p_odesys.names is None else '= {"%s"}' % '", "'.join(p_odesys.names)}; + std::vector p_odesys_names ${"" if p_odesys.names is None else '= {"%s"}' % '", "'.join(p_odesys.names)}; %if p_anon is not None: ${p_anon} %endif @@ -32,10 +55,11 @@ namespace odesys_anyode { template <> struct OdeSys: public AnyODE::${p_baseclass} { std::vector m_p; - std::vector m_p_cse; + std::vector m_cse; std::vector m_atol; std::vector m_upper_bounds; std::vector m_lower_bounds; + std::vector m_invar; std::vector m_invar0; realtype m_rtol; realtype m_get_dx_max_factor; @@ -95,33 +119,31 @@ namespace odesys_anyode { realtype get_dx_max_factor, bool error_outside_bounds, realtype max_invariant_violation, - std::vector special_settings) : - m_p_cse(${p_common['nsubs']}), m_atol(atol), m_rtol(rtol), m_get_dx_max_factor(get_dx_max_factor), - m_error_outside_bounds(error_outside_bounds), m_max_invariant_violation(max_invariant_violation), - m_special_settings(special_settings) { - m_p.assign(params, params + ${len(p_odesys.params) + p_odesys.ny if p_odesys.append_iv else 0}); - <% idx = 0 %> - %for cse_token, cse_expr in p_common['cses']: - %if cse_token.startswith('m_p_cse'): - ${cse_token} = ${cse_expr}; <% assert cse_token == 'm_p_cse[{0}]'.format(idx); idx += 1 %> - %else: - const auto ${cse_token} = ${cse_expr}; - %endif - %endfor - use_get_dx_max = (m_get_dx_max_factor > 0.0) ? ${'true' if p_get_dx_max else 'false'} : false; + std::vector special_settings) + : m_p(params, params + ${len(p_odesys.params) + p_odesys.ny if p_odesys.append_iv else 0}) + , m_cse(${p_common["n_cses"]}) + , m_atol(atol) + , m_invar(${0 if p_invariants is None else p_invariants["n_invar"]}) + , m_invar0(${0 if p_invariants is None else p_invariants["n_invar"]}) + , m_rtol(rtol) + , m_get_dx_max_factor(get_dx_max_factor) + , m_error_outside_bounds(error_outside_bounds) + , m_max_invariant_violation(max_invariant_violation) + , m_special_settings(special_settings) + { + ${p_common["cses"]} + use_get_dx_max = (m_get_dx_max_factor > 0.0) ? ${"true" if p_get_dx_max else "false"} : false; %if p_invariants is not None and p_support_recoverable_error: if (m_max_invariant_violation != 0.0){ - ${'' if p_odesys.append_iv else 'throw std::runtime_error("append_iv not set to True")'} + ${"" if p_odesys.append_iv else 'throw std::runtime_error("append_iv not set to True")'} const realtype * const y = params + ${len(p_odesys.params)}; - %for cse_token, cse_expr in p_invariants['cses']: - const auto ${cse_token} = ${cse_expr}; - %endfor - %for expr in p_invariants['exprs']: - m_invar0.push_back(${expr}); - %endfor - } + m_invar0.resize(${p_invariants["n_invar"]}); + //realtype * const out = m_invar0.data(); + ${p_invariants["cses"]} + ${p_invariants["assign"].all(assign_to=lambda i: "m_invar0[%d]" % i)} + } %endif - ${'\n '.join(p_constructor)} + ${"\n ".join(p_constructor)} } indextype OdeSys::get_ny() const { @@ -146,18 +168,30 @@ namespace odesys_anyode { AnyODE::Status OdeSys::rhs(realtype x, const realtype * const ANYODE_RESTRICT y, - realtype * const ANYODE_RESTRICT f) { - %if isinstance(p_rhs, str): + realtype * const ANYODE_RESTRICT out) { + %if getattr(p_odesys, "_nonnegative", False) and p_support_recoverable_error: + for (int i=0; i<${p_odesys.ny}; ++i) if (y[i] < 0) return AnyODE::Status::recoverable_error; + %endif + %if p_error_on_inf or p_error_on_nan: + for (int i=0; i<${p_odesys.ny}; ++i) { if (${_inf_or_nan('y[i]')}) { return ${_err_code()}; } } + + %endif + %if p_compensated_summation: + summation_cxx::RangedAccumulatorNeumaier f(${p_odesys.ny}); + f.init(out); + %else: + realtype * const f = out; + %endif + %if isinstance(p_rhs, str): ${p_rhs} - %else: - ${'AnyODE::ignore(x);' if p_odesys.autonomous_exprs else ''} - %for cse_token, cse_expr in p_rhs['cses']: - const auto ${cse_token} = ${cse_expr}; - %endfor + %else: + ${"AnyODE::ignore(x);" if p_odesys.autonomous_exprs else ""} + ${p_rhs["cses"]} + ${p_rhs["assign"].all(assign_to=lambda i: sympy.Symbol("f[%d]" % i))} + %if p_compensated_summation: + f.commit(); + %endif - %for i, expr in enumerate(p_rhs['exprs']): - f[${i}] = ${expr}; - %endfor this->nfev++; %if p_support_recoverable_error: if (m_error_outside_bounds){ @@ -166,7 +200,7 @@ namespace odesys_anyode { if (y[i] < m_lower_bounds[i]) { std::cerr << "Lower bound (" << m_lower_bounds[0] << ") for " << (p_odesys_names.size() ? p_odesys_names[i] : std::to_string(i)) - << " exceeded (" << y[i] << ") at x="<< x << "\n"; + << " violated (" << y[i] << ") at x="<< x << "\n"; return AnyODE::Status::recoverable_error; } } @@ -176,28 +210,29 @@ namespace odesys_anyode { if (y[i] > m_upper_bounds[i]) { std::cerr << "Upper bound (" << m_upper_bounds[0] << ") for " << (p_odesys_names.size() ? p_odesys_names[i] : std::to_string(i)) - << " exceeded (" << y[i] << ") at x="<< x << "\n"; + << " violates (" << y[i] << ") at x="<< x << "\n"; return AnyODE::Status::recoverable_error; } } } } - %if p_invariants is not None: + %if p_invariants is not None: if (m_max_invariant_violation != 0.0){ - %for cse_token, cse_expr in p_invariants['cses']: - const auto ${cse_token} = ${cse_expr}; - %endfor - %for idx, invar_expr in enumerate(p_invariants['exprs']): - if (std::abs(${invar_expr} - m_invar0[${idx}]) > ((m_max_invariant_violation > 0) ? m_max_invariant_violation : std::abs(m_max_invariant_violation*m_invar0[${idx}]) - m_max_invariant_violation)) { - std::cerr << "Invariant (${idx}) violation at x=" << x << "\n"; - return AnyODE::Status::recoverable_error; + ${p_invariants["cses"]} + ${p_invariants["assign"].all(assign_to=lambda i: "m_invar[%d]" % i)} + for (int idx=0; idx<${p_invariants["n_invar"]}; ++idx) { + if (std::abs(m_invar[idx] - m_invar0[idx]) > ((m_max_invariant_violation > 0) + ? m_max_invariant_violation + : std::abs(m_max_invariant_violation*m_invar0[idx]) /*- m_max_invariant_violation*/)) { + std::clog << "Invariant (" << idx << ") violation at x=" << x << "\n"; + return AnyODE::Status::recoverable_error; + } } - %endfor } - %endif + %endif %endif - %if getattr(p_odesys, '_nonnegative', False) and p_support_recoverable_error: - for (int i=0; i<${p_odesys.ny}; ++i) if (y[i] < 0) return AnyODE::Status::recoverable_error; + %if p_error_on_inf or p_error_on_nan: + for (int i=0; i<${p_odesys.ny}; ++i) { if (${_inf_or_nan('out[i]')}) { return ${_err_code()}; } } %endif return AnyODE::Status::success; %endif @@ -212,38 +247,32 @@ namespace odesys_anyode { } AnyODE::Status OdeSys::jtimes( - const realtype * const ANYODE_RESTRICT v, - realtype * const ANYODE_RESTRICT Jv, - realtype x, - const realtype * const ANYODE_RESTRICT y, - const realtype * const ANYODE_RESTRICT fy) { + const realtype * const ANYODE_RESTRICT v, + realtype * const ANYODE_RESTRICT out, + realtype x, + const realtype * const ANYODE_RESTRICT y, + const realtype * const ANYODE_RESTRICT fy) + { %if p_jtimes is not None: %if isinstance(p_jtimes, str): ${p_jtimes} %else: AnyODE::ignore(fy); // Currently we are not using fy (could be done through extensive pattern matching) - ${'AnyODE::ignore(x);' if p_odesys.autonomous_exprs else ''} - - %for cse_token, cse_expr in p_jtimes['cses']: - const auto ${cse_token} = ${cse_expr}; - %endfor - - %for i in range(p_odesys.ny): - <% curr_expr = p_jtimes['exprs'][i] %> - Jv[${i}] = ${curr_expr}; - %endfor + ${"AnyODE::ignore(x);" if p_odesys.autonomous_exprs else ""} + ${p_jtimes["cses"]} + ${p_jtimes["assign"].all()} %endif this->njvev++; return AnyODE::Status::success; %else: - AnyODE::ignore(v); AnyODE::ignore(Jv); AnyODE::ignore(x); + AnyODE::ignore(v); AnyODE::ignore(out); AnyODE::ignore(x); AnyODE::ignore(y); AnyODE::ignore(fy); return AnyODE::Status::unrecoverable_error; %endif } - %for order in ('cmaj', 'rmaj'): + %for order in ("cmaj", "rmaj"): AnyODE::Status OdeSys::dense_jac_${order}(realtype x, const realtype * const ANYODE_RESTRICT y, const realtype * const ANYODE_RESTRICT fy, @@ -256,29 +285,37 @@ namespace odesys_anyode { %else: // The AnyODE::ignore(...) calls below are used to generate code free from false compiler warnings. AnyODE::ignore(fy); // Currently we are not using fy (could be done through extensive pattern matching) - ${'AnyODE::ignore(x);' if p_odesys.autonomous_exprs else ''} - ${'AnyODE::ignore(y);' if (not any([yi in p_odesys.get_jac().free_symbols for yi in p_odesys.dep]) and - not any([yi in p_odesys.get_dfdx().free_symbols for yi in p_odesys.dep])) else ''} - - %for cse_token, cse_expr in p_jac_dense['cses']: - const auto ${cse_token} = ${cse_expr}; - %endfor + ${"AnyODE::ignore(x);" if p_odesys.autonomous_exprs else ""} + ${"AnyODE::ignore(y);" if (not any([yi in p_odesys.get_jac().free_symbols for yi in p_odesys.dep]) and + not any([yi in p_odesys.get_dfdx().free_symbols for yi in p_odesys.dep])) else ""} + %if p_error_on_inf or p_error_on_nan: + for (int i=0; i<${p_odesys.ny}; ++i) { if (${_inf_or_nan('y[i]')}) { return ${_err_code()}; } } + %endif + ${p_jac_dense["cses"]} %for i_major in range(p_odesys.ny): %for i_minor in range(p_odesys.ny): + <% - curr_expr = p_jac_dense['exprs'][i_minor, i_major] if order == 'cmaj' else p_jac_dense['exprs'][i_major, i_minor] - if curr_expr == '0' and p_jacobian_set_to_zero_by_solver: - continue - %> jac[ldim*${i_major} + ${i_minor}] = ${curr_expr}; + if order == "cmaj": + i = i_minor*p_odesys.ny + i_major + else: + i = i_major*p_odesys.ny + i_minor + if p_jac_dense["assign"].expr_is_zero(i) and p_jacobian_set_to_zero_by_solver: + continue + %> + ${p_jac_dense["assign"](i, assign_to=lambda _: "jac[ldim*%d + %d]" % (i_major, i_minor))} %endfor - %endfor if (dfdt){ - %for idx, expr in enumerate(p_jac_dense['dfdt_exprs']): - dfdt[${idx}] = ${expr}; + %for idx in range(p_odesys.ny**2, p_jac_dense["assign"].n): + ${p_jac_dense["assign"](idx, assign_to=lambda _: "dfdt[%d]" % (idx - p_odesys.ny**2))} %endfor } + %if p_error_on_inf or p_error_on_nan: + for (int i = 0; i < ${p_odesys.ny}; ++i) { if (${_inf_or_nan('y[i]')}) { return ${_err_code()}; } } + %endif + this->njev++; return AnyODE::Status::success; %endif @@ -299,46 +336,40 @@ namespace odesys_anyode { %elif isinstance(p_first_step, str): ${p_first_step} %else: - %for cse_token, cse_expr in p_first_step['cses']: - const realtype ${cse_token} = ${cse_expr}; - %endfor - ${'' if p_odesys.indep in p_odesys.first_step_expr.free_symbols else 'AnyODE::ignore(x);'} - ${'' if any([yi in p_odesys.first_step_expr.free_symbols for yi in p_odesys.dep]) else 'AnyODE::ignore(y);'} - return ${p_first_step['expr']}; + ${p_first_step["cses"]} + ${"" if p_odesys.indep in p_odesys.first_step_expr.free_symbols else "AnyODE::ignore(x);"} + ${"" if any([yi in p_odesys.first_step_expr.free_symbols for yi in p_odesys.dep]) else "AnyODE::ignore(y);"} + const auto ${p_first_step["assign"](0, assign_to=lambda _: "result")} + return result; %endif } AnyODE::Status OdeSys::sparse_jac_csc(realtype x, const realtype * const ANYODE_RESTRICT y, const realtype * const ANYODE_RESTRICT fy, - realtype * const ANYODE_RESTRICT data, + realtype * const ANYODE_RESTRICT out, indextype * const ANYODE_RESTRICT colptrs, indextype * const ANYODE_RESTRICT rowvals) { %if p_jac_sparse is not None: AnyODE::ignore(fy); // Currently we are not using fy (could be done through extensive pattern matching) - ${'AnyODE::ignore(x);' if p_odesys.autonomous_exprs else ''} - ${'AnyODE::ignore(y);' if (not any([yi in p_odesys.get_jac().free_symbols for yi in p_odesys.dep]) and - not any([yi in p_odesys.get_dfdx().free_symbols for yi in p_odesys.dep])) else ''} - %for cse_token, cse_expr in p_jac_sparse['cses']: - const auto ${cse_token} = ${cse_expr}; - %endfor - - %for i in range(p_odesys.nnz): - data[${i}] = ${p_jac_sparse['exprs'][i]}; - %endfor + ${"AnyODE::ignore(x);" if p_odesys.autonomous_exprs else ""} + ${"AnyODE::ignore(y);" if (not any([yi in p_odesys.get_jac().free_symbols for yi in p_odesys.dep]) and + not any([yi in p_odesys.get_dfdx().free_symbols for yi in p_odesys.dep])) else ""} + ${p_jac_sparse["cses"]} + ${p_jac_sparse["assign"].all()} %for i in range(p_odesys.nnz): - rowvals[${i}] = ${p_jac_sparse['rowvals'][i]}; + rowvals[${i}] = ${p_jac_sparse["rowvals"][i]}; %endfor %for i in range(p_odesys.ny + 1): - colptrs[${i}] = ${p_jac_sparse['colptrs'][i]}; + colptrs[${i}] = ${p_jac_sparse["colptrs"][i]}; %endfor this->njev++; return AnyODE::Status::success; %else: AnyODE::ignore(x); AnyODE::ignore(y); AnyODE::ignore(fy); - AnyODE::ignore(data); AnyODE::ignore(colptrs); AnyODE::ignore(rowvals); + AnyODE::ignore(out); AnyODE::ignore(colptrs); AnyODE::ignore(rowvals); return AnyODE::Status::unrecoverable_error; %endif } @@ -375,23 +406,26 @@ namespace odesys_anyode { } AnyODE::Status OdeSys::roots(realtype x, const realtype * const y, realtype * const out) { + + %if p_error_on_inf or p_error_on_nan: + if (${_inf_or_nan('x')}) { return ${_err_code()}; } + for (int i=0; i<${p_odesys.ny}; ++i) { if (${_inf_or_nan('y[i]')}) { return ${_err_code()}; } } + %endif + %if p_roots is None: AnyODE::ignore(x); AnyODE::ignore(y); AnyODE::ignore(out); - return AnyODE::Status::success; %elif isinstance(p_roots, str): ${p_roots} %else: - ${'' if any(p_odesys.indep in expr.free_symbols for expr in p_odesys.roots) else 'AnyODE::ignore(x);'} - - %for cse_token, cse_expr in p_roots['cses']: - const auto ${cse_token} = ${cse_expr}; - %endfor + ${"" if any(p_odesys.indep in expr.free_symbols for expr in p_odesys.roots) else "AnyODE::ignore(x);"} + ${p_roots["cses"]} + ${p_roots["assign"].all()} + %endif - %for i, expr in enumerate(p_roots['exprs']): - out[${i}] = ${expr}; - %endfor + %if p_error_on_inf or p_error_on_nan: + for (int i=0; inrev++; return AnyODE::Status::success; - %endif } } diff --git a/pyodesys/native/sources/odesys_anyode_template.cpp.mako b/pyodesys/native/sources/odesys_anyode_template.cpp.mako new file mode 100644 index 0000000..e69de29 diff --git a/pyodesys/native/sources/summation_cxx/accumulator.hpp b/pyodesys/native/sources/summation_cxx/accumulator.hpp new file mode 100644 index 0000000..ce6200f --- /dev/null +++ b/pyodesys/native/sources/summation_cxx/accumulator.hpp @@ -0,0 +1,95 @@ +#pragma once +#include "summation_cxx/compensated.hpp" +#include +#include +#include +#include + +namespace summation_cxx { + +template +struct Accumulator : public detail::Operators> { + static constexpr Compensation compensation_scheme { scheme }; + +private: + std::array data {}; + +public: + T& accum() { return data.data()[0]; } + T& carry() { return data.data()[1]; } + const T& accum() const { return data.data()[0]; } + const T& carry() const { return data.data()[1]; } + +public: + Accumulator() = default; + // cppcheck-suppress noExplicitConstructor + Accumulator(T accum) + { + data[0] = accum; + } + explicit Accumulator(T accum, T carry) + { + data[0] = accum; + data[1] = carry; + } + void clear() { data.clear(); } + + template + static constexpr T sum(const std::array& arr) + { + Accumulator ta {}; +#if defined(__clang__) +#pragma unroll 16 +#elif defined(__GNUC__) +#pragma GCC unroll 16 +#endif + for (const auto& e : arr) { + // cppcheck-suppress useStlAlgorithm + ta += e; + } + return ta.template to(); + } + +private: + template + static constexpr void from_(Accumulator& acu, U arg, Us... args) + { + acu += arg; + if constexpr (sizeof...(args) > 0) { + from_(acu, args...); + } + } + +public: + template + static constexpr Accumulator from(U arg, Us... args) + { + Accumulator acu { arg }; + if constexpr (sizeof...(args) > 0) { + Accumulator::from_(acu, args...); + } + return acu; + } +}; +template +using AccumulatorKahan = Accumulator; +template +using AccumulatorNeumaier = Accumulator; +template +using AccumulatorNeumaierSwap = Accumulator; +template +using AccumulatorTwoSum = Accumulator; + + +template +T pow(const Accumulator& base, T exponent) +{ + return std::pow(base.accum(), exponent); +} +template +T pow(const Accumulator& base, int exponent) +{ + return std::pow(base.accum(), static_cast(exponent)); +} + +} diff --git a/pyodesys/native/sources/summation_cxx/compensated.hpp b/pyodesys/native/sources/summation_cxx/compensated.hpp new file mode 100644 index 0000000..f6e79ed --- /dev/null +++ b/pyodesys/native/sources/summation_cxx/compensated.hpp @@ -0,0 +1,87 @@ +#pragma once +#ifdef __FAST_MATH__ +#error fast math enabled (/fp:fast, -ffast-math), this would negate compensation. +#endif +#include "summation_cxx/macros.hpp" +#include // std::size_t + +namespace summation_cxx { +enum class Compensation { + NONE, + KAHAN, // should be equivalent to FAST_TWO_SUM + NEUMAIER, + NEUMAIER_SWAP, + TWO_SUM, + FAST_TWO_SUM +}; + +namespace /* anonymous */ { + template + SXX_PREFER_INLINE void accum_kahan_destructive( + T& SXX_RESTRICT accu, + T& SXX_RESTRICT carry, + T& SXX_RESTRICT elem) + { + elem -= carry; + const T tmp = accu + elem; + carry = T { tmp - accu } - elem; + accu = tmp; + } + template + SXX_PREFER_INLINE void accum_kahan( + T& SXX_RESTRICT accu, + T& SXX_RESTRICT carry, + const T& SXX_RESTRICT elem) + { + T y = elem; + accum_kahan_destructive(accu, carry, y); + } + + template + SXX_PREFER_INLINE void accum_neumaier( + T& SXX_RESTRICT acm, + T& SXX_RESTRICT carry, + const T& SXX_RESTRICT elem) + { + SXX_NEUMAIER_ADD(acm, carry, elem, T, tmp, false); + } + + template + SXX_PREFER_INLINE void accum_neumaier_swap( + T& SXX_RESTRICT acm, + T& SXX_RESTRICT carry, + const T& SXX_RESTRICT elem) + { + // cppcheck-suppress redundantAssignment + SXX_NEUMAIER_ADD(acm, carry, elem, T, tmp, true); + } + + template + SXX_PREFER_INLINE void accum_two_sum( + T& SXX_RESTRICT accu, + T& SXX_RESTRICT carry, + const T& SXX_RESTRICT elem) + { + const T s = accu + elem; + const T ap = s - elem; + const T bp = s - ap; + const T da = accu - ap; + const T db = elem - bp; + carry += da + db; + accu = s; + } + + template + SXX_PREFER_INLINE void accum_fast_two_sum( + T& SXX_RESTRICT accu, + T& SXX_RESTRICT carry, + const T& SXX_RESTRICT elem) + { + const T s = accu + elem; + const T z = s - accu; + const T t = elem - z; + carry += t; + accu = s; + } +} +} diff --git a/pyodesys/native/sources/summation_cxx/impl.hpp b/pyodesys/native/sources/summation_cxx/impl.hpp new file mode 100644 index 0000000..063c1ea --- /dev/null +++ b/pyodesys/native/sources/summation_cxx/impl.hpp @@ -0,0 +1,177 @@ +#pragma once +#include +#include +#include + +namespace summation_cxx { +template +struct Accumulator; +template +struct AccuView; +namespace detail { + template + struct Operators; +} +} +namespace summation_cxx::detail { +template +struct Operators { + typedef T underlying_type; + typedef Accumulator accumulator_type; + typedef AccuView view_type; + +#define ACCUM(cv_qual) static_cast(this)->accum() +#define CARRY(cv_qual) static_cast(this)->carry() + template + U to() const + { + if constexpr (Derived::compensation_scheme == Compensation::KAHAN) { + return ACCUM(const); + } else if constexpr (Derived::compensation_scheme == Compensation::NEUMAIER + || Derived::compensation_scheme == Compensation::NEUMAIER_SWAP + || Derived::compensation_scheme == Compensation::TWO_SUM + || Derived::compensation_scheme == Compensation::FAST_TWO_SUM + ) { + if constexpr (sizeof(T) > sizeof(U)) { + return ACCUM(const) + CARRY(const); + } else { + return static_cast(ACCUM(const)) + static_cast(CARRY(const)); + } + } else { + assert(false); + } + return U { 0 } / U { 0 }; /* unreachable code, but would return NaN */ + } + Derived& operator+=(T arg) + { + if constexpr (scheme == Compensation::KAHAN) { + accum_kahan_destructive(ACCUM(), CARRY(), arg); + } else if constexpr (scheme == Compensation::NEUMAIER) { + accum_neumaier(ACCUM(), CARRY(), arg); + } else if constexpr (scheme == Compensation::NEUMAIER_SWAP) { + accum_neumaier_swap(ACCUM(), CARRY(), arg); + } else if constexpr (scheme == Compensation::TWO_SUM) { + accum_two_sum(ACCUM(), CARRY(), arg); + } else if constexpr (scheme == Compensation::FAST_TWO_SUM) { + accum_fast_two_sum(ACCUM(), CARRY(), arg); + } else { + assert(false); + } + return *(static_cast(this)); + } + Derived& operator-=(T arg) + { + Derived& self = *(static_cast(this)); + self += -arg; + return self; + } + + Derived& operator=(const T arg) + { + Derived& self = *(static_cast(this)); + ACCUM() = arg; + CARRY() = 0; + return self; + } + void operator/=(const T& arg) + { + ACCUM() /= arg; + CARRY() /= arg; + } + void operator*=(const T& arg) + { + const T ori {ACCUM()}; + ACCUM() *= arg; + CARRY() *= arg; + CARRY() += fma(ori, arg, -ACCUM()); // 2product + } + Derived& operator+=(const accumulator_type& other) + { + Derived& self = *(static_cast(this)); + self += other.accum(); + self /*CARRY()*/ += other.carry(); + return self; + } + Derived& operator-=(const accumulator_type& other) + { + Derived& self = *(static_cast(this)); + self -= other.accum(); + self /*CARRY()*/ -= other.carry(); + return self; + } + accumulator_type operator*(const T& arg) const + { + Derived cpy = *(static_cast(this)); + cpy *= arg; + return cpy; + } + accumulator_type operator*(const accumulator_type& other) const + { + return accumulator_type(ACCUM(const) * other.accum(), + CARRY(const) * other.accum() + ACCUM(const) * other.carry() + CARRY(const) * other.carry()); + } + accumulator_type operator/(const accumulator_type& other) const + { + const T denom = other.template to(); + return accumulator_type { ACCUM(const) / denom, CARRY(const) / denom }; + } + accumulator_type operator+(const accumulator_type& other) const + { + Derived cpy = *(static_cast(this)); + cpy += other; + return cpy; + } + accumulator_type operator+(const T& arg) const + { + Derived cpy = *(static_cast(this)); + cpy += arg; + return cpy; + } + accumulator_type operator-(const accumulator_type& other) const + { + Derived cpy = *(static_cast(this)); + cpy -= other; + return cpy; + } + accumulator_type operator+() const + { + return accumulator_type(ACCUM(const), CARRY(const)); + } + accumulator_type operator-() const + { + return accumulator_type(-ACCUM(const), -CARRY(const)); + } +#define SXX_COMP(OPER_) \ + bool operator OPER_(const T& arg) const \ + { \ + const Derived& self = *(static_cast(this)); \ + return self.template to() OPER_ arg; \ + } + SXX_COMP(<) + SXX_COMP(>) +}; +#undef SXX_COMP +#define SXX_COMMUTATIVE_OP(OP) \ + template \ + typename Derived::accumulator_type operator OP( \ + const typename Derived::underlying_type& arg_a, const Derived& arg_b) \ + { \ + return arg_b OP arg_a; /* multiplication is commutative */ \ + } +SXX_COMMUTATIVE_OP(*) +SXX_COMMUTATIVE_OP(+) +#define SXX_PROMOTING_OP(OP) \ + template \ + typename Derived::accumulator_type operator OP( \ + const typename Derived::underlying_type& arg_a, const Derived& arg_b) \ + { \ + return Derived { arg_a } OP arg_b; /* multiplication is commutative */ \ + } +SXX_PROMOTING_OP(-) + +#undef SXX_COMMUTATIVE_OP +#undef SXX_PROMOTING_OP +#undef ACCUM +#undef CARRY + +} diff --git a/pyodesys/native/sources/summation_cxx/macros.hpp b/pyodesys/native/sources/summation_cxx/macros.hpp new file mode 100644 index 0000000..af6c973 --- /dev/null +++ b/pyodesys/native/sources/summation_cxx/macros.hpp @@ -0,0 +1,84 @@ +#pragma once +#ifndef SXX_RESTRICT +#if defined(__GNUC__) +#define SXX_RESTRICT __restrict__ +#elif defined(_MSC_VER) && _MSC_VER >= 1400 +#define SXX_RESTRICT __restrict +// #elif defined (__STDC_VERSION__) && __STDC_VERSION__ >= 199901L +// #define SXX_RESTRICT restrict +#else +#define SXX_RESTRICT +#endif +#endif + +#ifndef SXX_PREFER_INLINE +#if defined(__GNUC__) +#define SXX_PREFER_INLINE __attribute__((flatten)) +#elif defined(_MSC_VER) && _MSC_VER >= 1400 +#define SXX_PREFER_INLINE __forceinline +#else +#define SXX_PREFER_INLINE +#endif +#endif + +#if defined(__cpp_concepts) && __cpp_concepts >= 201907L +#include +#define SXX_FWD_IT_CONCEPT std::forward_iterator +#define SXX_RND_IT_CONCEPT std::random_access_iterator +#define SXX_UNSIGNED_INTEGRAL std::unsigned_integral +#else +#define SXX_FWD_IT_CONCEPT typename +#define SXX_RND_IT_CONCEPT typename +#define SXX_UNSIGNED_INTEGRAL typename +#endif + +// Math macros to support e.g. __float128 without std lib support: +#ifndef SXX_ABS +#define SXX_ABS(x) (((x) < 0) ? -(x) : (x)) +#endif + +#ifndef SXX_NEUMAIER_BRANCH +// see test/bench.cpp +#define SXX_NEUMAIER_BRANCH 1 +#endif + +#if defined(SXX_NEUMAIER_SWAP) +#error "API has changed, update your compilation flags accordingly" +#endif + +#define SXX_SWP_TMP_(ACCUM, CARRY, ELEM, TYPE, TMP, DO_SWAP) \ + TYPE TMP; \ + if (DO_SWAP) { \ + if ((CARRY) == 0 && (ACCUM) != 0 && SXX_ABS(ELEM) > (1u<<20)*SXX_ABS(ACCUM)) { \ + TMP = ACCUM; \ + ACCUM = CARRY; \ + CARRY = TMP; \ + } \ + } \ + TMP = (ACCUM) + (ELEM); + +#if SXX_NEUMAIER_BRANCH == 1 +#define SXX_NEUMAIER_ADD(ACCUM, CARRY, ELEM, TYPE, TMP, DO_SWAP) \ + do { \ + SXX_SWP_TMP_(ACCUM, CARRY, ELEM, TYPE, TMP, DO_SWAP) \ + if (SXX_ABS(TMP) > SXX_ABS(ELEM)) { \ + CARRY += TYPE { (ACCUM) - (TMP) } + (ELEM); \ + } else { \ + CARRY += TYPE { (ELEM) - (TMP) } + (ACCUM); \ + } \ + ACCUM = (TMP); \ + } while (0) +#else +#define SXX_NEUMAIER_ADD(ACCUM, CARRY, ELEM, TYPE, TMP, DO_SWAP) \ + do { \ + SXX_SWP_TMP_(ACCUM, CARRY, ELEM, TYPE, TMP, DO_SWAP) \ + T SXX_dat_[2] = { \ + T { (ELEM) - (TMP) } + (ACCUM), \ + T { (ACCUM) - (TMP) } + (ELEM) \ + }; \ + CARRY += SXX_dat_[SXX_ABS(TMP) > SXX_ABS(ELEM)]; \ + ACCUM = (TMP); \ + } while (0) +#endif +#undef SXX_CXX_SWP_TMP_ +#define SXX_NEUMAIER_FINALIZE(ACCUM, CARRY) ((ACCUM) + (CARRY)) diff --git a/pyodesys/native/sources/summation_cxx/ranged.hpp b/pyodesys/native/sources/summation_cxx/ranged.hpp new file mode 100644 index 0000000..464842e --- /dev/null +++ b/pyodesys/native/sources/summation_cxx/ranged.hpp @@ -0,0 +1,126 @@ +#pragma once +#include +#include // std::memset +#include // std::make_unique +#include +#include +#include + +namespace summation_cxx { +template +struct RangedAccumulator { + typedef T underlying_type; + typedef std::conditional_t, T, U> target_type; + typedef Accumulator accumulator_type; + typedef AccuView view_type; + +private: + target_type* tgt {}; + std::unique_ptr storage {}; + std::size_t sz {}; + bool cumulative {}; + +public: + RangedAccumulator() = default; + RangedAccumulator(std::size_t sz) + : storage(std::make_unique(sz * 2)) + , sz(sz) + { + } + void init(target_type* target, bool cumulative = false) + { + tgt = target; + this->cumulative = cumulative; + if (sz > 0 /* UB to call memset over zero bytes. */) { + // doing this only makes sense if commit() is not always called. + std::memset(storage.get(), 0x00, sizeof(underlying_type) * sz * 2); + } + } + view_type operator[](std::size_t idx) + { + return view_type { &storage[idx * 2] }; + } + const view_type operator[](std::size_t idx) const + { + return view_type { &storage[idx * 2] }; + } + + void commit() const + { +#if defined(SUMMTNCXX_DISTRUST_OPTIMIZING_COMPILERS) +#define SXX_OUTPUT(OP) \ + if constexpr (scheme == Compensation::KAHAN) { \ + this->tgt[i] OP this->storage[i * 2]; \ + } else if constexpr (scheme == Compensation::NEUMAIER || + Derived::compensation_scheme == Compensation::NEUMAIER_SWAP) + { + this->tgt[i] OP this->storage[i * 2] + this->storage[i * 2 + 1]; + } + else + { + assert(false); + } +#else +#define SXX_OUTPUT(OP) this->tgt[i] OP(*this)[i].template to(); +#endif +#define SXX_LOOP \ + std::size_t i = 0; \ + i < this->sz; \ + ++i + if (cumulative) { + for (SXX_LOOP) { + SXX_OUTPUT(+=) + } + } else { + for (SXX_LOOP) { + SXX_OUTPUT(=) + } + } +#undef SXX_LOOP +#undef SXX_OUTPUT + } +}; +template +using RangedAccumulatorKahan = RangedAccumulator; +template +using RangedAccumulatorNeumaier = RangedAccumulator; +template +using RangedAccumulatorNeumaierSwap = RangedAccumulator; +template +using RangedAccumulatorTwoSum = RangedAccumulator; +template +using RangedAccumulatorFastTwoSum = RangedAccumulator; + +/// Simplifies writing generic code against ranged.hpp, no compensation: +template +struct RangedUncompensatedView { + typedef T underlying_type; + typedef std::conditional_t, T, U> target_type; + typedef T accumulator_type; + typedef target_type& view_type; + +protected: + target_type* tgt {}; + std::size_t sz {}; + +public: + RangedUncompensatedView() = default; + RangedUncompensatedView(std::size_t sz) + : sz(sz) + { + } + void init(target_type* target, bool cumulative = false) + { + tgt = target; + if (!cumulative && (sz > 0 /* UB to call memset over zero bytes. */)) { + // may be skipped e.g. if we know target is already zero-initialized + std::memset(target, 0x00, sizeof(target_type) * sz); + } + } + view_type operator[](std::size_t idx) + { + return tgt[idx]; + } + void commit() const { } // no-op +}; +} diff --git a/pyodesys/native/sources/summation_cxx/view.hpp b/pyodesys/native/sources/summation_cxx/view.hpp new file mode 100644 index 0000000..3a0228c --- /dev/null +++ b/pyodesys/native/sources/summation_cxx/view.hpp @@ -0,0 +1,45 @@ +#pragma once +#include + +namespace summation_cxx { + +template +struct AccuView : public detail::Operators> { + static constexpr Compensation compensation_scheme { scheme }; + +private: + T* ptr; + +public: + T& accum() { return ptr[0]; } + T& carry() { return ptr[1]; } + const T& accum() const { return ptr[0]; } + const T& carry() const { return ptr[1]; } + +public: + AccuView() = delete; + AccuView & operator=(const AccuView&) = delete; + using detail::Operators>::operator=; + // cppcheck-suppress noExplicitConstructor + AccuView(T* data) + : ptr(data) + { + assert(data); + } + Accumulator deepcopy() + { + return Accumulator { this->accum(), this->carry() }; + } +}; + +template +using AccuViewKahan = AccuView; +template +using AccuViewNeumaier = AccuView; +template +using AccuViewNeumaierSwap = AccuView; +template +using AccuViewTwoSum = AccuView; +template +using AccuViewFastTwoSum = AccuView; +} diff --git a/pyodesys/native/symcse/__init__.py b/pyodesys/native/symcse/__init__.py new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/pyodesys/native/symcse/__init__.py @@ -0,0 +1 @@ + diff --git a/pyodesys/native/symcse/compensated.py b/pyodesys/native/symcse/compensated.py new file mode 100644 index 0000000..10f0749 --- /dev/null +++ b/pyodesys/native/symcse/compensated.py @@ -0,0 +1,396 @@ +#!/usr/bin/env python3 +"""Use this utility for symbolic treatment. + +Example usage: + +$ python3 -m symcse demo-compensated-py +$ python3 -m symcse demo-compensated-c --index 3 | clang-format --style=Google | batcat -pl C + +""" +from collections import defaultdict +from functools import reduce +from operator import add, mul +from sympy import ( + Abs, + Add, + And, + Eq, + Expr, + Lt, + Mul, + Ne, + numbered_symbols, + Piecewise, + postorder_traversal, + Symbol, + Tuple, +) +from sympy.codegen import Assignment, aug_assign, CodeBlock +from sympy.codegen.ast import AssignmentBase, Token, While, break_ + +from .core import NullTransformer +from .ordered_add import ordered_add + + +def If(cond, body): + return While(cond, CodeBlock(*body, break_)) + +class _CompensatedAdd(Token, Expr): + """Represents KBN compensated summation.""" + + _fields = __slots__ = ("terms", "accum", "carry", "temp") + _construct_terms = staticmethod(lambda args: Tuple(*args)) + + def _ccode(self, printer): + terms = ", ".join(map(printer._print, self.terms)) + return f"NA({terms} /*{str(self.accum)[:-1]}*/)" + + def to_statements(self, existing, expanded, *, transients, **kwargs): + """Transform into statements.""" + neum, ordinary = [], [] + for term in self.terms: + ex, other = [], [] + for factor in term.as_ordered_factors(): + if factor in existing: + ex.append(factor) + else: + other.append(factor) + if len(ex) == 1: # TODO: handle >1 + neum.append((reduce(mul, other) if other else 1, existing[ex[0]])) + else: + ordinary.append(term) + st = [] + if neum: + if True: + st.append(Assignment(self.accum, neum[0][0] * neum[0][1].accum)) + ordinary.extend([oth * na.accum for oth, na in neum[1:]]) + else: + st.append( + Assignment(self.accum, sum(oth * na.accum for oth, na in neum)) + ) + + st.append(Assignment(self.carry, sum(oth * na.carry for oth, na in neum))) + else: + st.append(Assignment(self.accum, ordinary.pop(0))) + st.append(Assignment(self.carry, 0)) + + for elem in ordinary: + nops = elem.count_ops() + if nops == 1 and isinstance(elem, Mul): + nops = max(nops, (-elem).count_ops()) + if nops: + # for fs in elem.free_symbols: + # if fs not in existing: + # continue + # ex = existing[fs] + # if ex not in expanded: + # continue + # e = elem.subs(fs, ex.accum + ex.carry) + # _accu, _pure_carry = [], [] + # for term in e.as_ordered_terms(): + # if ex.carry in term and ex.accum not in term: + # _pure_carry.append(term) + # else: + # _accum.append(term) + # if _pure_carry: + # ... + # ea = elem.subs(fs, ex.accum) + # ec = elem.sbus(fs, ex.carry) + tr = next(transients) + st.append(Assignment(tr, elem)) + elem = tr + st.extend( + self._impl_add(self.accum, self.carry, elem, self.temp, **kwargs) + ) + expanded.add(self) + return st + + def finalize(self): + """Close the summation.""" + return self._impl_finalize(self.accum, self.carry) + + @staticmethod + def _impl_finalize(accum, carry): + return Add(accum, carry) + +class _NeumaierAdd(_CompensatedAdd): + _impl_add_kws = ('do_swap',) + + @staticmethod + def _impl_add(accum, carry, elem, temp, do_swap=False): + """Perform Kahan-Babuska-Neumaier addition.""" + big_temp = ordered_add(ordered_add(accum, -temp), elem) + big_elem = ordered_add(ordered_add(elem, -temp), accum) + abs_elem = Abs(elem) + pw = Piecewise((big_temp, Abs(temp) > abs_elem), (big_elem, True)) + statements = [ + Assignment(temp, accum + elem), + aug_assign(carry, "+", pw), + Assignment(accum, temp), + ] + if do_swap: + return [ + If( + And(Eq(carry, 0), Ne(accum, 0), Lt(Abs(accum), abs_elem)), + [ + Assignment(temp, accum), + Assignment(accum, carry), + Assignment(carry, temp), + ], + ) + ] + statements + else: + return statements + + +class _TwoSumAdd(_CompensatedAdd): + _impl_add_kws = () + + @staticmethod + def _impl_add(accum, carry, elem, temp): + """Perform 2Sum addition.""" + a_prim = ordered_add(temp, -elem) + b_prim = ordered_add(temp, -a_prim) + delta_a = ordered_add(accum, -a_prim) + delta_b = ordered_add(elem, -b_prim) + statements = [ + Assignment(temp, accum + elem), + aug_assign(carry, "+", ordered_add(delta_a, delta_b)), + Assignment(accum, temp) + ] + return statements + + +class _FastTwoSumAdd(_CompensatedAdd): + _impl_add_kws = () + + @staticmethod + def _impl_add(accum, carry, elem, temp): + """Perform Fast2Sum addition.""" + statements = [ + Assignment(temp, ordered_add(accum, elem)), + aug_assign(carry, "+", ordered_add(elem, -ordered_add(temp, -accum))), + Assignment(accum, temp), + ] + return statements + + +class _CompensationTransformer(NullTransformer): + """Transform Add instances in CSEs to use compensated sum. + + Parameters + ---------- + up_to_debug: int, [0-100] + Code is guaranteed to compile at levels 0 (no passes, no compensation) + and 100 (all passes). + """ + + _CompAdd = None # Token + + def __init__( + self, + repl, + red, + *, + tmp_pfx="t", + trs_pfx="r", + csum_pfx="n", + up_to_debug=100, + limit=3, + parent=None, + ignore=None, + kw_to_statements=None + ): + self.repl = repl + self.red = red + self.limit = limit + self.parent = parent + self.ignore = ignore + + self.created = {} + self._all_accum = {} + self._all_carry = {} + self._all_tempv = {} + self.expanded = set() + self._analysis = defaultdict(int) + self._tmp_var = numbered_symbols(tmp_pfx) + self._trs_var = numbered_symbols(trs_pfx) + self._csum_var = numbered_symbols(csum_pfx) + self.passes = [] + for p in filter(lambda n: n.startswith("_pass_"), dir(self)): + null, rest = p.split("_pass_") + assert null == "" + num, *_ = rest.split("_") + if len(_) == 0: + continue + if int(num) <= up_to_debug: + self.passes.append(getattr(self, p)) + + self.kw_to_statements = kw_to_statements or {} + self.statements, self.final_exprs = self._pipeline() + + def remapping_for_arrayification(self, template="m_glob[{0}]"): + remapping = {} + i = 0 + for st in self.statements: + if st.lhs in remapping or st.lhs in self._all_tempv: + continue + # if st.lhs in self._all_accum or st.lhs in self._all_carry: + remapping[st.lhs] = Symbol(template.format(i), real=True) + i = i + 1 + return remapping + + def _mk_Comp(self, terms, lhs): + pfx = str(next(self._tmp_var)) if lhs is None else str(lhs) + accum = Symbol(pfx + "a", real=True) + carry = Symbol(pfx + "c", real=True) + tempv = Symbol(pfx + "t", real=True) + na = self._CompAdd(terms, accum, carry, tempv) + self._all_accum[accum] = na + self._all_carry[carry] = na + self._all_tempv[tempv] = na + return na + + @classmethod + def _is_Comp(cls, x): + return isinstance(x, cls._CompAdd) + + def _single_pass(self, statements, pass_): + new_stmts = [] + for st in statements: + if isinstance(st, AssignmentBase): + new_rhs = pass_(st.lhs, st.rhs, statements=new_stmts) + if st.lhs not in self.created: + new_stmts.append(st.__class__(st.lhs, new_rhs)) + elif hasattr(st, "body"): + assert isinstance(st.body, CodeBlock) + new_body = CodeBlock(*self._single_pass(st.body.args, pass_)) + new_args = ( + new_body if attr == "body" else getattr(st, attr) + for attr in st.__slots__ + ) + new_stmts.append(st.__class__(*new_args)) + else: + new_stmts.append(st) # no-op (e.g. BreakToken instance) + return new_stmts + + def _pipeline(self): + statements = [Assignment(*lr) for lr in self.repl] + final_exprs = self.red + for pass_ in self.passes: + statements = self._single_pass(statements, pass_) + final_exprs = [pass_(None, e, statements=statements) for e in final_exprs] + return statements, final_exprs + + def _pass_05_analysis(self, lhs, rhs, *, statements): + if rhs.is_Add: + if lhs is not None: + self._analysis[lhs] += len(rhs.args) + for term in rhs.args: + if term in self._analysis: + self._analysis[term] += len(rhs.args) - 1 + return rhs + + def _pass_10_create_nodes(self, lhs, rhs, *, statements, debug=False): + new_rhs = rhs + if self.ignore: + if lhs in self.ignore or any(rhs.has(ig) for ig in self.ignore): + return new_rhs + + while True: + for _add in filter(lambda x: x.is_Add, postorder_traversal(new_rhs)): + score = self._analysis.get(lhs, 0) + reduce( + add, [self._analysis.get(k, 1) for k in _add.args] + ) + if score >= self.limit or any(self._is_Comp(arg) for arg in _add.args): + na = self._mk_Comp(_add.args, lhs) + if _add is rhs and lhs is not None: + key = lhs + + else: + key = next(self._csum_var) + self.created[key] = na + new_rhs = new_rhs.xreplace({_add: key}) + break + else: + return new_rhs + + assert False + + def _pass_50_to_stmnts(self, lhs, rhs, *, statements): + for csum in map(self.created.get, postorder_traversal(rhs)): + if csum is None: + continue + self._pass_50_to_stmnts(lhs, csum.terms, statements=statements) + if csum not in self.expanded: + statements.extend( + csum.to_statements( + self.created, + self.expanded, + transients=self._trs_var, + **self.kw_to_statements + ) + ) + return rhs + + def _pass_60_xrepl(self, lhs, rhs, *, statements): + return rhs.xreplace(self.created) + + def _group(self, x): + all_accum, all_carry, rest = [], [], [] + for term in x.args: + if term in self.created: + all_accum.append(self.created[term].accum) + all_carry.append(self.created[term].carry) + elif term in self._all_accum: + all_accum.append(term) + elif term in self._all_carry: + all_carry.append(term) + else: + rest.append(term) + result = [] + if all_accum: + result.append(ordered_add(*all_accum)) + if all_carry: + result.append(ordered_add(*all_carry)) + if rest: + result += rest + result = reduce(add, result) + return result + + def _has(self, term): + return ( + term in self.created or term in self._all_accum or term in self._all_carry + ) + + def _pass_95_group(self, lhs, rhs, *, statements): + new_rhs = rhs.replace( + lambda s: s.is_Add and any(self._has(t) for t in s.args), self._group + ) + return new_rhs + + def _pass_90_fin(self, lhs, rhs, *, statements): + return rhs.replace(lambda x: self._is_Comp(x), lambda x: x.finalize()) + + +class NeumaierTransformer(_CompensationTransformer): + _CompAdd = _NeumaierAdd + + def __init__(self, *args, do_swap=False, **kwargs): + #kw_to_statements = kwargs.pop(kw_to_statements, {}) + # ... + super().__init__(*args, kw_to_statements=dict(do_swap=do_swap), **kwargs) + +class TwoSumTransformer(_CompensationTransformer): + _CompAdd = _TwoSumAdd + +class FastTwoSumTransformer(_CompensationTransformer): + _CompAdd = _FastTwoSumAdd + + +compensated_transformers = { + 'nocomp': NullTransformer, + 'kbn': NeumaierTransformer, + '2sum': TwoSumTransformer, + 'fast2sum': FastTwoSumTransformer, +} diff --git a/pyodesys/native/symcse/core.py b/pyodesys/native/symcse/core.py new file mode 100644 index 0000000..99776af --- /dev/null +++ b/pyodesys/native/symcse/core.py @@ -0,0 +1,85 @@ +from collections import defaultdict +from sympy import Symbol +from sympy.codegen.ast import Assignment, AssignmentBase, float64, value_const, Variable + + +class NullTransformer: + """Perform no transformation.""" + + def __init__(self, repl, red, parent=None, ignore=None): + """From CSE result.""" + self.statements = [Assignment(lhs, rhs) for lhs, rhs in repl] + self.final_exprs = red + + def statements_with_declarations( + self, + pred=lambda s: "[" not in s.name, + type_=(lambda lhs, rhs: (float64, rhs)) + ): + """Conditionally add declarations.""" + return self._declare(self.statements, pred=pred, type_=type_) + + @staticmethod + def _is_const(ctx_is_const, statements): + for st in statements: # filter(lambda s: isinstance(s, AssignmentBase), stmts): + if isinstance(st, AssignmentBase): + ctx_is_const[st.lhs] += 1 if type(st) is Assignment else 2 + elif hasattr(st, "body"): + NullTransformer._is_const(ctx_is_const, st.body.args) + else: + pass # no-op + + @staticmethod + def _as_decl(ctx_seen, statements, is_const, *, pred, type_, lift=False): + new_stmnts, lifted = [], [] + for st in statements: + if isinstance(st, Assignment) and st.lhs not in ctx_seen and pred(st.lhs): + ctx_seen.add(st.lhs) + t, new_rhs = type_(st.lhs, st.rhs) + var = Variable( + st.lhs, + type=t, + attrs=[value_const] if is_const[st.lhs] == 1 else [], + ) + if lift: + lifted.append(var.as_Declaration()) + else: + st = var.as_Declaration(value=new_rhs) + elif hasattr(st, "body"): + body_stmnts, body_lifted = NullTransformer._as_decl( + ctx_seen, st.body.args, is_const, pred=pred, type_=type_, lift=True + ) + if lift: + lifted.extend(body_lifted) + else: + new_stmnts.extend(body_lifted) + else: + pass # no-op, e.g. BreakToken() + new_stmnts.append(st) + return new_stmnts, lifted + + def _declare(self, stmts, *, pred=lambda x: True, type_=(lambda lhs, rhs: (float64, rhs))): + seen = set() + is_const = defaultdict(int) + self._is_const(is_const, stmts) + new_stmnts, lifted = self._as_decl( + seen, stmts, is_const, pred=pred, lift=False, type_=type_ + ) + assert len(lifted) == 0 + return new_stmnts + + def apply_remapping(self, remapping): + """Replace in statements & final_exprs.""" + self.statements = [st.xreplace(remapping) for st in self.statements] + self.final_exprs = [e.xreplace(remapping) for e in self.final_exprs] + + def remapping_for_arrayification(self, template="m_glob[{0}]"): + """Create a replacement dictionary.""" + remapping = {} + i = 0 + for st in self.statements: + if st.lhs in remapping: + continue + remapping[st.lhs] = Symbol(template.format(i), real=True) + i = i + 1 + return remapping diff --git a/pyodesys/native/symcse/groupwise.py b/pyodesys/native/symcse/groupwise.py new file mode 100644 index 0000000..0487afa --- /dev/null +++ b/pyodesys/native/symcse/groupwise.py @@ -0,0 +1,188 @@ +"""Handle common cses among groups of code.""" + +from collections import defaultdict +from functools import reduce +from operator import add +import numpy as np + +import sympy +from sympy import numbered_symbols +from sympy.codegen.ast import Assignment, Variable, float64 +from sympy.codegen.rewriting import create_expand_pow_optimization + +from .core import NullTransformer +from .util import CPrinter, SymbolicBackend, ccode + +expand_up_to_3 = create_expand_pow_optimization(3) + + +def pre_process(expr): + """Simplify, expand & factor. + + Example: expr.simplify(rational=True).expand().factor()""" + return expr.factor() + + +def post_process(expr): + """Expand low integer powers.""" + return expand_up_to_3(expr) + + +class GroupwiseCSE: + """Eliminate common sub-expressions from groups of expressions.""" + + def __init__( + self, + groups, + *, + common_cse_template="common_cse{0}", + common_ignore=(), + code_printer=None, + subsd=None, + transformers=None, + type_=(lambda lhs, rhs: (float64, rhs)), + pre_process=None, + post_process=None, + backend=None, + transformer_kws=None + ): + """ + Parameters + ---------- + groups : list like + \\*\\*kwargs : see code for now. + + """ + if transformers is None: + transformers = defaultdict(lambda: NullTransformer) + if code_printer is None: + code_printer = CPrinter() + if transformer_kws is None: + transformer_kws = defaultdict(dict) + self._code_printer = code_printer + self._type = type_ + self._keys, _values = zip(*groups.items()) + self._spans = np.cumsum([0] + list(map(len, _values))) + if backend is None: + backend = SymbolicBackend() + self.backend = backend + _all_values = reduce(add, map(list, _values)) + _all_exprs = list(map(pre_process, _all_values) if pre_process else _all_values) + _all_exprs = [ + e.replace(lambda s: s.is_Symbol, lambda s: sympy.Symbol(s.name, real=True)) + for e in _all_exprs + ] + common_ignore = [sympy.Symbol(ig.name, real=True) for ig in common_ignore] + for e in _all_exprs: + for fs in e.free_symbols: + if not fs.is_real: + # Switching between symengine/sympy is tricky, consistently using + # real=True with Symbols allows us to assume real=True for SymEngine + # symbols (which lack support for assumptions at the time of writing). + raise NotImplementedError( + "Only use explicitly real valued symbols." + ) + repls, reds = self._common_cse( + _all_exprs, + ignore=common_ignore, + symbols=numbered_symbols("cse_comm_locl", real=True), + ) + if post_process: + repls = [(s, post_process(e)) for s, e in repls] + reds = [post_process(e) for e in reds] + self._comm_tformr = transformers[None]( + repls, reds, ignore=common_ignore, **transformer_kws[None] + ) + remap = self._comm_tformr.remapping_for_arrayification( + template=common_cse_template + ) + self._comm_tformr.apply_remapping(remap) + _subsd = {sympy.Symbol(k.name, real=True): v for k, v in (subsd or {}).items()} + self._comm_tformr.apply_remapping(_subsd) + self.n_remapped = len(remap) + + assert len(self._comm_tformr.final_exprs) == len(reds) + del reds + self._per_g_tformrs = self._get_g_tformrs( + self._comm_tformr, + transformers=transformers, + transformer_kws=transformer_kws, + post_process=post_process, + subsd=_subsd, + ) + + @property + def keys(self): + """Retrieve the keys of the groups.""" + return self._keys + + def render(self, x): + """Generate a code string.""" + return self._code_printer.doprint(x) + + def _common_cse(self, all_exprs, **kwargs): + repls, reds = self.backend.cse(all_exprs, **kwargs) + cse_symbols = numbered_symbols("cse_t", real=True) # local temporaries + comm_subs = {} + for lhs, rhs in repls: + for expr in reds: + if lhs in expr.free_symbols: + comm_subs[lhs] = next(cse_symbols) + break + return ( + [(lhs.xreplace(comm_subs), rhs.xreplace(comm_subs)) for lhs, rhs in repls], + [r.xreplace(comm_subs) for r in reds], + ) + + def _get_g_tformrs( + self, comm_tformr, *, transformers, transformer_kws, post_process, subsd + ): + per_g = {} + for i, k in enumerate(self._keys): + g_repls, g_exprs = self.backend.cse( + comm_tformr.final_exprs[slice(*self._spans[i : i + 2])], + symbols=numbered_symbols("cse", real=True), + ) + if post_process: + g_repls = [(s, post_process(e)) for s, e in g_repls] + g_exprs = [post_process(e) for e in g_exprs] + g_tformr = transformers[k]( + g_repls, g_exprs, parent=comm_tformr, **transformer_kws[k] + ) + g_tformr.apply_remapping(subsd) + per_g[k] = g_tformr + + return per_g + + @staticmethod + def _declare(stmts, *, pred, type_=float64): + seen = set() + result = [] + for st in stmts: + if isinstance(st, Assignment) and st.lhs not in seen and pred(st.lhs): + seen.add(st.lhs) + st = Variable(st.lhs, type=type_).as_Declaration(value=st.rhs) + result.append(st) + return result + + def common_statements(self, declare=False, type_=None): + """Initialize the common sub-expressions among the groups.""" + if declare: + return self._comm_tformr.statements_with_declarations( + pred=declare, type_=type_ or self._type + ) + else: + return self._comm_tformr.statements + + def statements(self, gk, declare=False, type_=None): + """Initialize the group specific sub-expressions.""" + if declare: + return self._per_g_tformrs[gk].statements_with_declarations( + pred=declare, type_=type_ or self._type + ) + else: + return self._per_g_tformrs[gk].statements + + def exprs(self, gk): + """Retrieve the resulting expressions of the group named ``key``.""" + return self._per_g_tformrs[gk].final_exprs diff --git a/pyodesys/native/symcse/ordered_add.py b/pyodesys/native/symcse/ordered_add.py new file mode 100644 index 0000000..221478a --- /dev/null +++ b/pyodesys/native/symcse/ordered_add.py @@ -0,0 +1,17 @@ +"""Utilities used in the package.""" + +from .util import UnevaluatedRealPropagatingExpr + + +def ordered_add(*args): + """Preserve the order of terms.""" + if len(args) == 0: + return 0 + elif (len(args)) == 1: + return args[0] + a, b, *rest = args + ue = UnevaluatedRealPropagatingExpr(a + b) + if len(rest) == 0: + return ue + else: + return ordered_add(ue, *rest) diff --git a/pyodesys/native/symcse/util.py b/pyodesys/native/symcse/util.py new file mode 100644 index 0000000..5f9c9e0 --- /dev/null +++ b/pyodesys/native/symcse/util.py @@ -0,0 +1,294 @@ +"""Utilities used in the package.""" +import os +import sys +from functools import reduce +from operator import mul, add + +import sympy +from sympy.codegen.ast import float80, String, Token +from sympy.printing.c import C99CodePrinter +from sympy.printing.cxx import CXX17CodePrinter +from sympy.printing.pycode import PythonCodePrinter + +try: + import symengine as se +except ModuleNotFoundError: + se = None +else: + from symengine.lib.symengine_wrapper import sympy2symengine + + +def partn(it, pred): + """Partition an iterable into False/True groups based on predicate.""" + result = ([], []) + for elem in it: + result[pred(elem)].append(elem) + return result + + +def prod(it): + return reduce(mul, it) + + +def map_multi(funcs, its): + return (f(it) for f, it in zip(funcs, its)) + + +def idty(x): + """Identity operattion (idempotent).""" + return x + + +class UnevaluatedRealPropagatingExpr(sympy.UnevaluatedExpr): + """Propagate .is_real, but nothings else from wrapped expression.""" + + def _eval_is_real(self): + return self.args[0].is_real + +class OperatorDot(Token): + _fields = __slots__ = ('lhs', 'rhs') + #_construct_lhs = String + #_construct_rhs = String + +class _UnevaluatedExprPrinterMixin: + def _print_UnevaluatedExpr(self, arg): + return "(%s)" % super()._print_UnevaluatedExpr(arg) + + def _print_Integer(self, arg): + if abs(arg) > 2 ** 53: + return self._print(sympy.Float(arg)) + else: + return super()._print_Integer(arg) + + def _print_OperatorDot(self, arg): + return self._print(arg.lhs) + '.' + self._print(arg.rhs) + + @staticmethod + def _replace_re(arg): + if isinstance(arg, sympy.UnevaluatedExpr) and arg.args[0].is_real: + return UnevaluatedRealPropagatingExpr(arg.args[0]) + else: + return sympy.re(arg) + + def doprint(self, x): + result = super().doprint(x.replace(sympy.re, self._replace_re)) + return result + + +class CXXPrinter(_UnevaluatedExprPrinterMixin, CXX17CodePrinter): + """Patched C++-printer (https://github.com/sympy/sympy/issues/21955).""" + + +class CPrinter(_UnevaluatedExprPrinterMixin, C99CodePrinter): + """Patched C-printer (https://github.com/sympy/sympy/issues/21955).""" + + +class PyPrinter(_UnevaluatedExprPrinterMixin, PythonCodePrinter): + """Patched Python-printer (https://github.com/sympy/sympy/issues/21955).""" + + +def _cse_symengine(exprs, *, se2sympy, ignore=(), symbols=None, **kwargs): + exprs = list(exprs) + se_exprs = [sympy2symengine(e, raise_error=True) for e in exprs] + repl, red = se.cse(se_exprs, **kwargs) + repl = [(se2sympy(lhs), se2sympy(rhs)) for lhs, rhs in repl] + red = [se2sympy(e) for e in red] + if ignore: + ignore = [se2sympy(sympy2symengine(ig)) for ig in ignore] + + def has_ig(e): + return any(e.has(ig) for ig in ignore) + + def is_only_ig(e): + return all(s in ignore for s in e.free_symbols) + + keep, reintro = {}, {} + for lhs, rhs in repl: + rem = rhs.xreplace(reintro) # remaining + if not has_ig(rem): + for k, v in keep.items(): + if v == rem: + reintro[lhs] = k + break + else: + keep[lhs] = rem + continue + + candidates = [] + for part in sympy.preorder_traversal(rem): + if not has_ig(part): + candidates.append(part) + elif part.is_Add: + good, bad = partn(part.args, has_ig) + if good: + candidates.append(sum(good)) + elif part.is_Mul: + good, bad = partn(part.args, has_ig) + if good: + candidates.append(prod(good)) + if candidates: + part = sorted(candidates, key=sympy.count_ops)[-1] + if part.count_ops() == 0 or len(part.free_symbols) == 0: + reintro[lhs] = rem + else: + for k, v in keep.items(): + if part == k or part == v: + reintro[lhs] = rem + break + else: + keep[lhs] = part + reintro[lhs] = rem.subs(part, lhs) + else: + reintro[lhs] = rem + repl = list(keep.items()) + red = [e.xreplace(reintro) for e in red] + # recreated = [e.subs(reversed(repl)) for e in red] + # for e1, e2 in zip(recreated, exprs): + # assert (e1-e2).simplify() == 0 + + if symbols is not None: + remap = {} + updated = [] + for lhs, rhs in repl: + remap[lhs] = next(symbols) + updated.append((remap[lhs], rhs.xreplace(remap))) + repl = updated + red = [e.xreplace(remap) for e in red] + return repl, red + + +class SymbolicBackend: + """Allow optional use of SymEngine.""" + + def __init__(self, use_symengine=None, kw_cse=None, assume_real=True): + """Initialize a symbolic backend. + + Parameters + ---------- + use_symengine: bool + kw_cse: keywords to use in CSE + + """ + if use_symengine is None: + _req_backend = os.environ.get("SYMCSE_SYMBOLIC_BACKEND", "").lower() + if _req_backend == "symengine": + use_symengine = True + elif _req_backend == "sympy": + use_symengine = False + elif _req_backend == "": + use_symengine = False # se is not None + else: + raise ValueError("Unknown SYMCSE_SYMBOLIC_BACKEND: %s" % _req_backend) + if use_symengine and se is None: + raise ValueError("symengine missing (pip install symengine)") + self.use_symengine = use_symengine + self.kw_cse = kw_cse + self.assume_real = assume_real + + def Symbol(self, name): + if self.use_symengine: + return se.Symbol( + name + ) # https://github.com/symengine/symengine.py/issues/286 + else: + return sympy.Symbol(name, real=self.assume_real) + + def sympy2se(self, x): + if hasattr(x, "_sympy_"): + return x # looks like that's already a SymEngine object + return sympy2symengine(x, raise_error=True) + + def se2sympy(self, x): + if not hasattr(x, "_sympy_"): + assert isinstance(x, sympy.Basic) + return x # looks like that's already a SymPy object + tmp = x._sympy_() + return tmp.xreplace( + {s: sympy.Symbol(s.name, real=self.assume_real) for s in tmp.free_symbols} + ) + + def cse(self, exprs, **kwargs): + """Perform common sub-expression elimination.""" + exprs = list(exprs) + new_kw = dict(**(self.kw_cse or {}), **kwargs) + if self.use_symengine: + return _cse_symengine(exprs, se2sympy=self.se2sympy, **new_kw) + else: + return sympy.cse(exprs, **new_kw) + + def _elems_colmaj(self, mtx): + if not mtx: + return [] + result = reduce(add, mtx.T.tolist()) + if self.use_symengine: + return [self.se2sympy(e) for e in result] + else: + return result + + def _col_vec(self, v): + args = len(v), 1, tuple(v) + if self.use_symengine: + return se.Matrix(*args) + else: + return sympy.Matrix(*args) + + def jacobian(self, exprs, wrt): + args = len(exprs), 1, exprs + if self.use_symengine: + return se.Matrix(*args).jacobian(self._col_vec(wrt)) + else: + return sympy.Matrix(*args).jacobian(wrt) + + def _LU(self, mat): + if self.use_symengine: + _L, _U = self.sympy2se(mat).LU() + L, U = map(self.se2sympy, [_L, _U]) + else: + L, U, piv = mat.LUdecomposition() + if piv: + raise NotImplementedError("TODO") + for i in range(L.cols): + for j in range(i, L.cols): + L[i, j] = U[i, j] + return L + + def matmul(self, A, b): + if self.use_symengine: + return A.mul_matrix(b) + else: + return A @ b + + def lambdify(self, args, exprs): + if self.use_symengine: + return se.Lambdify(args, exprs) + else: + return sympy.lambdify(args, exprs) + + +class SymbolicBackendWithDisabledCSE(SymbolicBackend): + def cse(self, exprs, **kwargs): + return [], exprs + + +default_settings = dict(math_macros={}, type_mappings={float80: "long double"}) + + +def ccode(arg, **kwargs): + for k, v in default_settings.items(): + if k not in kwargs: + kwargs[k] = v + p = CPrinter(settings=kwargs) + return p.doprint(arg) + + +def cxxcode(arg, **kwargs): + for k, v in default_settings.items(): + if k not in kwargs: + kwargs[k] = v + p = CXXPrinter(settings=kwargs) + return p.doprint(arg) + + +def pycode(arg): + p = PyPrinter() + return p.doprint(arg) diff --git a/pyodesys/native/sympy_interface.py b/pyodesys/native/sympy_interface.py new file mode 100644 index 0000000..a548d32 --- /dev/null +++ b/pyodesys/native/sympy_interface.py @@ -0,0 +1,441 @@ +#!/usr/bin/env python3 +"""Use this utility for symbolic treatment. + +Example usage: + +$ python3 sympy_interface.py demo1 | clang-format --style=Google | batcat -pl C + +""" +import pprint +from collections import defaultdict +from functools import reduce +from operator import add +from sympy import ( + Abs, Add, Basic, ccode, cse, exp, Expr, numbered_symbols, Piecewise, pi, + postorder_traversal, preorder_traversal, pycode, Symbol, Tuple +) +from sympy.codegen import Assignment, aug_assign, CodeBlock +from sympy.codegen.ast import Token, Variable, float64 +from sympy.printing.precedence import PRECEDENCE, precedence + + +class OrderedAdd(Expr): + """Printed with parenthesis, useful for floating point math.""" + + nargs = -1 + precedence = PRECEDENCE["Add"] - 1 + + def _ccode(self, printer): + return " + ".join(printer.parenthesize(arg, precedence(self)) for arg in self.args) + + _pythoncode = _ccode + + +class _NeumaierAdd(Token, Expr): + """Represents KBN compensated summation.""" + + __slots__ = ('terms', 'accum', 'carry', 'temp') + _construct_terms = staticmethod(lambda args: Tuple(*args)) + + def _ccode(self, printer): + terms = ", ".join(map(printer._print, self.terms)) + return f"NA({terms} /*{str(self.accum)[:-1]}*/)" + + def to_statements(self, existing, expanded): + """Transform into statements.""" + neum, ordinary = [], [] + for term in self.terms: + # print(f"/*{term, term in existing}*/") + if term in existing: #.values(): + # neum.append(term) + neum.append(existing[term]) + else: + ordinary.append(term) + # print(neum, ordinary, existing) + st = [] + if neum: + st.append(Assignment(self.accum, sum(na.accum for na in neum))) + st.append(Assignment(self.carry, sum(na.carry for na in neum))) + for na in neum: + expanded.add(na) + else: + st.append(Assignment(self.accum, ordinary.pop(0))) + st.append(Assignment(self.carry, 0)) + + for elem in ordinary: + st.extend(_NeumaierAdd._impl_add(self.accum, self.carry, elem, self.temp)) + expanded.add(self) + return st + + def finalize(self): + """Close the summation.""" + return self._impl_finalize(self.accum, self.carry) + + @staticmethod + def _impl_add(accum, carry, elem, temp): + """Perform Kahan-Babuska-Neumaier addition.""" + big_temp = OrderedAdd(OrderedAdd(accum, -temp), elem) + big_elem = OrderedAdd(OrderedAdd(elem, -temp), accum) + pw = Piecewise((big_temp, Abs(temp) > Abs(elem)), (big_elem, True)) + return [ + Assignment(temp, accum + elem), + aug_assign(carry, '+', pw), + Assignment(accum, temp) + ] + + @staticmethod + def _impl_finalize(accum, carry): + return Add(accum, carry) + + +class _NeumaierTransformer: + """Transform Add instances in CSEs to use compensated sum. + + Parameters + ---------- + up_to: int, [0-100] + Code is guaranteed to compile at levels 0 (no passes, no compensation) + and 100 (all passes). + """ + + def __init__(self, repl, red, *, tmp_pfx="t", neu_pfx="n", up_to=100, limit=3): + self.repl = repl + self.red = red + self.limit = limit + + self.created = {} + self._all_accum = set() + self._all_carry = set() + self.expanded = set() + self._analysis = defaultdict(int) + self._tmp_var = numbered_symbols(tmp_pfx) + self._neu_var = numbered_symbols(neu_pfx) + self.passes = [] + for p in filter(lambda n: n.startswith('_pass_'), dir(self)): + null, rest = p.split("_pass_") + assert null == "" + num, *_ = rest.split("_") + if len(_) == 0: + continue + if int(num) <= up_to: + self.passes.append(getattr(self, p)) + + self.statements, self.final_exprs = self._pipeline() + # print("\ncreated: ", self.created) + # print("") + # print(self._analysis) + # print(self.statements) + + def statements_and_expressions(self): + return self._declare(self.statements), self.final_exprs + + def _mk_Neu(self, terms, lhs): + pfx = str(next(self._tmp_var)) if lhs is None else str(lhs) + accum = Symbol(pfx+'a') + carry = Symbol(pfx+'c') + tempv = Symbol(pfx+'t') + self._all_accum.add(accum) + self._all_carry.add(carry) + return _NeumaierAdd(terms, accum, carry, tempv) + + @staticmethod + def _is_Neu(x): + return isinstance(x, _NeumaierAdd) + + @staticmethod + def _declare(stmts, *, type_=float64): + seen = set() + result = [] + for st in stmts: + if isinstance(st, Assignment) and st.lhs not in seen: + seen.add(st.lhs) + st = Variable(st.lhs, type=type_).as_Declaration(value=st.rhs) + result.append(st) + return result + + def _pipeline(self): + statements = [Assignment(*lr) for lr in self.repl] + final_exprs = self.red + for pass_ in self.passes: + new_stmts, new_exprs = [], [] + for st in statements: + new_rhs = pass_(st.lhs, st.rhs, statements=new_stmts) + if st.lhs not in self.created: # not Neumaier, (aug)assign: + new_stmts.append(st.__class__(st.lhs, new_rhs)) + for expr in final_exprs: + new_exprs.append(pass_(None, expr, statements=new_stmts)) + statements = new_stmts + final_exprs = new_exprs + # print(f"/*{self._analysis}*/") + return statements, final_exprs + + def _pass_05_analysis(self, lhs, rhs, *, statements): + if rhs.is_Add: + if lhs is not None: + self._analysis[lhs] += len(rhs.args) + for term in rhs.args: + if term in self._analysis: + self._analysis[term] += len(rhs.args) - 1 + return rhs + + def _pass_10_create_nodes(self, lhs, rhs, *, statements, debug=False): + #new_rhs = rhs.xreplace(self.created) + new_rhs = rhs + while True: + for _add in filter(lambda x: x.is_Add, postorder_traversal(new_rhs)): + score = self._analysis.get(lhs, 0) + reduce(add, [self._analysis.get(k, 1) for k in _add.args]) + if score >= self.limit or any(self._is_Neu(arg) for arg in _add.args): + na = self._mk_Neu(_add.args, lhs) + if _add is rhs and lhs is not None: + key = lhs + + else: + key = next(self._neu_var) + #statements.append(Assignment(key, na)) + self.created[key] = na + # print(new_rhs, _add) + new_rhs = new_rhs.xreplace({_add: key}) + break + else: + return new_rhs + + assert False + # if debug and lhs is None: + # print("new_rhs: ", new_rhs) + #return + # if debug and lhs is None: + # print("") + # print("_add: ", _add) + # print("") + # print("na: ", na) + # print("") + # print("new_rhs: ", new_rhs) + + #return new_rhs + + # def _pass_11_(self, lhs, rhs, *, statements): + # return self._pass_10_create_nodes(lhs, rhs, statements=statements, debug=True) + + def _pass_50_to_stmnts(self, lhs, rhs, *, statements): + for arg in postorder_traversal(rhs): + # print(arg, self._is_Neu(arg)) + if arg in self.created: + for t in list(self.created[arg].terms)+[arg]: + if t in self.created: + if arg not in self.expanded: + statements.extend(self.created[t].to_statements(self.created, self.expanded)) + # elif self._is_Neu(arg): + # print(".") + # statements.extend(arg.to_statements(self.created, self.expanded)) + return rhs + + def _pass_60_xrepl(self, lhs, rhs, *, statements): + return rhs.xreplace(self.created) + + def _group(self, x): + # print(self.created) + # print(x)### + all_accum, all_carry = [], [] + for term in x.args: + if term in self.created: + all_accum.append(self.created[term].accum) + all_carry.append(self.created[term].carry) + elif term in self._all_accum: + all_accum.append(term) + elif term in self._all_carry: + all_carry.append(term) + else: + all_accum.append(term) + # print('//', all_accum, all_carry) + if all_accum and all_carry: + return OrderedAdd(OrderedAdd(*all_accum), OrderedAdd(*all_carry)) + else: + return OrderedAdd(*(all_accum+all_carry)) + + def _pass_95_group(self, lhs, rhs, *, statements): + # Not sure if this pass is still needed + return rhs.replace(lambda s: s.is_Add and any( + t in self.created or t in self._all_accum or t in self._all_carry for t in s.args + ), self._group) + + def _pass_90_fin(self, lhs, rhs, *, statements): + return rhs.replace(lambda x: self._is_Neu(x), lambda x: x.finalize()) + + +def _compensated_code(case, **kwargs): + repl, red = cse(case.exprs) + nm = _NeumaierTransformer(repl, red, **kwargs) + statements, new_exprs = nm.statements_and_expressions() + for i, new_expr in enumerate(new_exprs): + statements.append(Assignment(Symbol("out[%d]" % i), new_expr)) + return CodeBlock(*statements) + + +class Case: + """test case.""" + + vals_ref = None + + @classmethod + def check(cls): + """Perform self-check.""" + return _NeumaierTransformer(*cse(cls.exprs)) + + +class _Case0(Case): + x = list(map(Symbol, ["y[%d]" % i for i in range(3)])) + a = reduce(add, x[:3]) + exprs = [ + a, + a+1, + ] + + +class _Case1(Case): + x = list(map(Symbol, ["y[%d]" % i for i in range(6)])) + a = reduce(add, x[:3]) + b = reduce(add, x[3:6]) + exprs = [ + a, + b, + a + b, + a + b + 1, + ] + + +class _Case2(Case): + x = list(map(Symbol, ["y[%d]" % i for i in range(6)])) + a = reduce(add, x[:3]) + b = reduce(add, x[3:6]) + exprs = [ + a, + b, + a + b + 1, + a + b + 1 + pi, + ] + + +class _Case3(Case): + x = list(map(Symbol, ["y[%d]" % i for i in range(10)])) + a = reduce(add, x[:3]) + b = reduce(add, x[3:6]) + c = reduce(add, x[6:9]) + exprs = [ + a + 1, + b + 2, + c + 3, + a + b + x[9], + a + b + c + x[8] + ] + + +class _Case4(Case): + x = list(map(Symbol, ["y[%d]" % i for i in range(3)])) + a = 1*x[0] + b = 2*x[1] + c = 3*x[2] + exprs = [ + a + b, + a + b + c, + a + b + c + 1 + ] + + +def real(name): + return Symbol(name, real=True) + + +class _Case5(Case): + x = list(map(real, ["y[%d]" % i for i in range(5)])) + a = x[0]*x[1]**2 + x[2]**3 + x[3]**4 + b = 3*x[3] + 4*x[4] + c = exp(1*x[1] + 2*x[2] + b) + d = (1*x[1] + 2*x[2] + 3*x[3])/b + exprs = [a, a/42, b+c, c*d, b+exp(x[0])] + + +class _Case6(Case): + x = list(map(Symbol, ["y[%d]" % i for i in range(2)])) + a = reduce(add, x[:]) + exprs = [ + a, + a+1, + ] + + +class _Case7(Case): # no Neumaier + x = list(map(Symbol, ["y[%d]" % i for i in range(2)])) + a = reduce(add, x[:]) + exprs = [ + a, + a/3, + ] + + @classmethod + def check(cls): + tr = super().check() + assert(len(tr.statements) <= 1) + + +class _Case8(Case): + x = list(map(Symbol, ["y[%d]" % i for i in range(6)])) + a = reduce(add, x[:3]) + b = reduce(add, x[3:]) + exprs = [ + a, + b, + a+b + ] + +class _Case9(Case): + x = list(map(Symbol, ["y[%d]" % i for i in range(6)])) + a = reduce(add, x[:3]) + b = reduce(add, x[3:]) + exprs = [ + exp(a) + exp(b), + a+b + ] + +class _Case10(Case): + x = list(map(Symbol, ["y[%d]" % i for i in range(6)])) + a = reduce(add, x[:3]) + b = reduce(add, x[3:]) + exprs = [ + a/b, + 2*a + 3*b + ] + +class _Case11(Case): + x = list(map(Symbol, ["y[%d]" % i for i in range(4)])) + exprs = [sum(x)] + vals_ref = ([1e14, 1.0, -1e14, 1.0], [2.0]) + + +cases = [v for k, v in globals().items() + if k.startswith("_Case") and Case in v.__mro__] + +template_c = ( + "#include \n" + "void f(double * const out, const double * const y) {%s}" +) +template_py = "def f(out, y):\n %s" + +if __name__ == '__main__': + + def demo_c(index=0, up_to=100, limit=3): + """Demo C-code.""" + case = cases[index] + print('/* ' + str(case.exprs)) + code_block = _compensated_code(case, up_to=up_to, limit=limit) + print('*/') + # print(code_block) + src = ccode(code_block) + print(template_c % src) + + def demo_py(index=0, up_to=100, limit=3): + """Demo Python-code.""" + src = pycode(_compensated_code(cases[index], up_to=up_to, limit=limit)) + print(template_py % src.replace('\n', '\n ')) + + import argh + argh.dispatch_commands([demo_c, demo_py]) diff --git a/pyodesys/native/tests/_tests.py b/pyodesys/native/tests/_tests.py index d3e980c..1ca5372 100644 --- a/pyodesys/native/tests/_tests.py +++ b/pyodesys/native/tests/_tests.py @@ -13,12 +13,15 @@ ) from pyodesys.tests.bateman import bateman_full # analytic, never mind the details from pyodesys.tests.test_symbolic import decay_rhs, decay_dydt_factory, _get_decay3 +from pyodesys.native.symcse.util import SymbolicBackendWithDisabledCSE sp = import_('sympy') def _test_NativeSys(NativeSys, **kwargs): - native = NativeSys.from_callback(vdp_f, 2, 1) + native = NativeSys.from_callback( + vdp_f, 2, 1#, native_code_kw=dict(save_temp=True) + ) assert native.ny == 2 assert len(native.params) == 1 xout, yout, info = native.integrate([0, 1, 2], [1, 0], params=[2.0], **kwargs) @@ -95,8 +98,11 @@ class TransformedNativeSys(TransformedSys, NativeSys): assert np.allclose(yout, ref, rtol=rtol*forgive, atol=atol*forgive) -def _test_Decay_nonnegative(NativeSys): - odesys = NativeSys.from_other(_get_decay3(lower_bounds=[0]*3)) +def _test_Decay_nonnegative(NativeSys, use_cse, compensated): + odesys = NativeSys.from_other( + _get_decay3(lower_bounds=[0]*3), + native_code_kw=dict(groupwise_kw=None if use_cse else {'backend': SymbolicBackendWithDisabledCSE()}, compensated_summation=compensated), + ) y0, k = [3., 2., 1.], [3.5, 2.5, 0] xout, yout, info = odesys.integrate([1e-10, 1], y0, k, integrator='native') ref = np.array(bateman_full(y0, k, xout - xout[0], exp=np.exp)).T @@ -242,10 +248,19 @@ def _test_NativeSys__first_step_cb(NativeSys, forgive=20): def _test_NativeSys__first_step_cb_source_code(NativeSys, log10myconst, should_succeed, forgive=20, **kwargs): dec3 = _get_decay3() - odesys = NativeSys.from_other(dec3, namespace_override={ - 'p_first_step': 'return good_const()*y[0];', - 'p_anon': 'double good_const(){ return std::pow(10, %.5g); }' % log10myconst - }, namespace_extend={'p_includes': ['']}) + odesys = NativeSys.from_other( + dec3, + native_code_kw=dict( + save_temp=True, + namespace_override={ + 'p_first_step': 'AnyODE::ignore(x); return good_const()*y[0];', + 'p_anon': 'double good_const(){ return std::pow(10, %.5g); }' % log10myconst + }, + namespace_extend={ + 'p_includes': {''} + } + ), + ) y0, k = [.7, 0, 0], [1e23, 2, 3.] xout, yout, info = odesys.integrate(5, y0, k, integrator='native', **kwargs) ref = np.array(bateman_full(y0, k, xout - xout[0], exp=np.exp)).T @@ -275,9 +290,9 @@ def roots(t, y, p, backend): def _test_NativeSys__get_dx_max_source_code(NativeSys, forgive=20, **kwargs): dec3 = _get_decay3() - odesys = NativeSys.from_other(dec3, namespace_override={ + odesys = NativeSys.from_other(dec3, native_code_kw=dict(namespace_override={ 'p_get_dx_max': """AnyODE::ignore(y); return (1.0e-4 * x + 1.0e-3);""", - }) + })) y0, k = [.7, 0, 0], [7., 2, 3.] xout, yout, info = odesys.integrate(1, y0, k, integrator='native', get_dx_max_factor=1.0, **kwargs) @@ -321,7 +336,8 @@ def _test_NativeSys__dep_by_name__single_varied(NativeSys): def _test_return_on_error_success(NativeSys): k, y0 = [4, 3], (5, 4, 2) - native = NativeSys.from_callback(decay_rhs, len(k)+1, len(k), namespace_override={ + native = NativeSys.from_callback(decay_rhs, len(k)+1, len(k), native_code_kw=dict( + namespace_override={ 'p_rhs': """ f[0] = -m_p[0]*y[0]; f[1] = m_p[0]*y[0] - m_p[1]*y[1]; @@ -330,7 +346,7 @@ def _test_return_on_error_success(NativeSys): this->nfev++; return AnyODE::Status::success; """ - }) + })) xout = [0, 0.2, 0.4, 0.6, 0.8, 1.0] result = native.integrate(xout, y0, k, atol=1e-11, rtol=1e-11, return_on_error=True, dx_max=.05) nreached = result.info['nreached'] @@ -340,7 +356,7 @@ def _test_return_on_error_success(NativeSys): assert np.allclose(result.yout[:nreached, :], ref, rtol=1e-8, atol=1e-8) -def _test_render_native_code_cse(NativeSys): +def _test_render_native_code_cse(NativeSys, compensated): # regression test taken from chempy from pyodesys.symbolic import SymbolicSys from sympy import symbols, log, exp @@ -419,7 +435,13 @@ def _solve(odesys, **kwargs): ref = _solve(symbolic, **kw) assert ref.info['success'] - native = NativeSys.from_other(symbolic) # <-- regression test, optional: save_temp=True + native = NativeSys.from_other( + symbolic, + native_code_kw=dict( + #save_temp=True, + compensated_summation=compensated + ), + ) # regression test: sol = _solve(native, **kw) assert sol.info['success'] diff --git a/pyodesys/native/tests/test_cvode.py b/pyodesys/native/tests/test_cvode.py index 71131ee..955b3a6 100644 --- a/pyodesys/native/tests/test_cvode.py +++ b/pyodesys/native/tests/test_cvode.py @@ -1,6 +1,7 @@ # -*- coding: utf-8 -*- from __future__ import (absolute_import, division, print_function) +from itertools import product import numpy as np import pytest @@ -88,8 +89,9 @@ def test_chained_multi_native_nonnegative(): @requires('pycvodes') -def test_Decay_nonnegative(): - _test_Decay_nonnegative(NativeSys) +@pytest.mark.parametrize('use_cse,compensated', product([False, True], [False, True])) +def test_Decay_nonnegative(use_cse, compensated): + _test_Decay_nonnegative(NativeSys, use_cse=use_cse, compensated=compensated) @requires('pycvodes') @@ -148,6 +150,9 @@ def roots(t, y, p, backend): assert result.info['success'] == True # noqa assert np.min(np.abs(result.xout - 1)) < 1e-11 + f_out = odesys.rhs(np.array(0.0, dtype=np.float64), np.array([2.0]), np.array([], dtype=np.float64)) + assert f_out.size == 1 and np.all(f_out == 2.0) + @requires('sym', 'pycvodes') @pytest.mark.parametrize('idx', [0, 1, 2]) @@ -238,8 +243,9 @@ def f(t, y, p): @requires('pycvodes', 'sympy') -def test_render_native_cse_regression(): - _test_render_native_code_cse(NativeSys) +@pytest.mark.parametrize('compensated', [False, True]) +def test_render_native_cse_regression(compensated): + _test_render_native_code_cse(NativeSys, compensated=compensated) @requires('sym', 'pycvodes') diff --git a/pyodesys/native/tests/test_gsl.py b/pyodesys/native/tests/test_gsl.py index 1e69285..ac23972 100644 --- a/pyodesys/native/tests/test_gsl.py +++ b/pyodesys/native/tests/test_gsl.py @@ -1,6 +1,7 @@ # -*- coding: utf-8 -*- from __future__ import (absolute_import, division, print_function) +from itertools import product import pytest from pyodesys.util import requires @@ -74,8 +75,9 @@ def test_chained_multi_native(reduced): @requires('pygslodeiv2') -def test_Decay_nonnegative(): - _test_Decay_nonnegative(NativeSys) +@pytest.mark.parametrize('use_cse,compensated', product([False, True], [False, True])) +def test_Decay_nonnegative(use_cse, compensated): + _test_Decay_nonnegative(NativeSys, use_cse=use_cse, compensated=compensated) @requires('pygslodeiv2') diff --git a/pyodesys/native/tests/test_odeint.py b/pyodesys/native/tests/test_odeint.py index 3339e4c..050dfb7 100644 --- a/pyodesys/native/tests/test_odeint.py +++ b/pyodesys/native/tests/test_odeint.py @@ -1,6 +1,7 @@ # -*- coding: utf-8 -*- from __future__ import (absolute_import, division, print_function) +from itertools import product from pyodesys.util import requires import pytest @@ -63,8 +64,9 @@ def test_PartiallySolved_symmetric_native_multi(): @requires('pyodeint') -def test_Decay_nonnegative(): - _test_Decay_nonnegative(NativeSys) +@pytest.mark.parametrize('use_cse,compensated', product([False, True], [False, True])) +def test_Decay_nonnegative(use_cse, compensated): + _test_Decay_nonnegative(NativeSys, use_cse=use_cse, compensated=compensated) @requires('pyodeint') diff --git a/pyodesys/plotting.py b/pyodesys/plotting.py index c8dfbf1..07edddb 100644 --- a/pyodesys/plotting.py +++ b/pyodesys/plotting.py @@ -6,6 +6,12 @@ import numpy as np +DEFAULT_C = tuple( + f'#{c}' for c in ( + '2ca02c', 'c49c94', 'f7b6d2', '7f7f7f', '9edae5', '98df8a', 'bcbd22', '9467bd', 'c5b0d5', 'd62728', + '8c564b', 'c7c7c7', '1f77b4', 'e377c2', 'dbdb8d', 'ff9896', '17becf', 'aec7e8', 'ffbb78', 'ff7f0e' + ) +) def _set_scale(cb, argstr): if argstr.count(';') == 0: @@ -23,8 +29,7 @@ def _latex_from_dimensionality(dim): def plot_result(x, y, indices=None, plot_kwargs_cb=None, ax=None, ls=('-', '--', ':', '-.'), - c=('tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple', - 'tab:brown', 'tab:pink', 'tab:gray', 'tab:olive', 'tab:cyan', 'black'), + c=DEFAULT_C, m=('o', 'v', '8', 's', 'p', 'x', '+', 'd', 's'), m_lim=-1, lines=None, interpolate=None, interp_from_deriv=None, names=None, latex_names=None, xlabel=None, ylabel=None, @@ -119,6 +124,7 @@ def plot_kwargs_cb(idx, lines=False, markers=False, labels=None): ax.fill_between(x, _y[:, idx] - yerr[:, idx], _y[:, idx] + yerr[:, idx], facecolor=clr, alpha=.3) if isinstance(yscale, str) and 'linthresh' in yscale: + yscale = yscale.replace('linthreshy', 'linthresh') arg, kw = yscale.split(';') thresh = eval('dict(%s)' % kw)['linthresh'] ax.axhline(thresh, linewidth=.5, linestyle='--', color='k', alpha=.5) diff --git a/pyodesys/results.py b/pyodesys/results.py index c121077..ee37a25 100644 --- a/pyodesys/results.py +++ b/pyodesys/results.py @@ -18,7 +18,8 @@ def __init__(self, xout, yout, params, info, odesys): self.params = params self.info = info self.odesys = odesys - self.names = self.odesys.names + # names not always present: + self.names = getattr(self.odesys, 'names', None) def copy(self): return Result(self.xout.copy(), self.yout.copy(), self.params.copy(), @@ -179,7 +180,7 @@ def stiffness(self, xyp=None, eigenvals_cb=None): return (np.abs(singular_values).max(axis=-1) / np.abs(singular_values).min(axis=-1)) - def _plot(self, cb, x=None, y=None, legend=None, **kwargs): + def _plot(self, cb, x=None, y=None, legend=None, sort_legend=False, **kwargs): if x is None: x = self.xout if y is None: @@ -199,9 +200,26 @@ def _plot(self, cb, x=None, y=None, legend=None, **kwargs): if legend is None: if (kwargs.get('latex_names') or None) is not None or (kwargs['names'] or None) is not None: legend = True - return cb(x, y, legend=legend, **kwargs) + ax = cb(x, y, legend=legend, **kwargs) + if sort_legend: + handles, labels = ax.get_legend_handles_labels() + # sort both labels and handles by labels + # if kwargs.get('names', None) is None: + # sel = slice(None) + # else: + # sel = [n in kwargs['names'] for n in self.odesys.names] + reorder = np.argsort(y[-1, kwargs.get('indices', slice(None))])[::-1] + if isinstance(legend, dict): + legend_kw = legend + elif not legend: + return ax + else: + legend_kw = dict() + ax.legend([handles[i] for i in reorder], [labels[i] for i in reorder], **legend_kw) + return ax - def plot(self, info_vlines_kw=None, between=None, deriv=False, title_info=0, **kwargs): + def plot(self, info_vlines_kw=None, between=None, deriv=False, title_info=0, + **kwargs): """ Plots the integrated dependent variables from last integration. Parameters @@ -217,6 +235,7 @@ def plot(self, info_vlines_kw=None, between=None, deriv=False, title_info=0, **k \\*\\*kwargs: See :func:`pyodesys.plotting.plot_result` """ + if between is not None: if 'x' in kwargs or 'y' in kwargs: raise ValueError("x/y & between given.") @@ -231,10 +250,12 @@ def plot(self, info_vlines_kw=None, between=None, deriv=False, title_info=0, **k if 'y' in kwargs: raise ValueError("Cannot give both deriv=True and y.") kwargs['y'] = self.odesys.f_cb(*self._internals()) + ax = self._plot(plot_result, **kwargs) + if title_info: ax.set_title( - (self.odesys.description or '') + + (getattr(self.odesys, 'description', None) or '') + ', '.join( (['%d steps' % self.info['n_steps']] if self.info.get('n_steps', -1) >= 0 else []) + [ @@ -267,11 +288,12 @@ def calc_invariant_violations(self, xyp=None): val = invar(*(xyp or self._internals())) return val - val[0, :] - def plot_invariant_violations(self, **kwargs): + def plot_invariant_violations(self, *, apply_abs: bool=True, **kwargs): viol = self.calc_invariant_violations() - abs_viol = np.abs(viol) + if apply_abs: + viol = np.abs(viol) invar_names = self.odesys.all_invariant_names() - return self._plot(plot_result, x=self._internal('xout'), y=abs_viol, names=invar_names, + return self._plot(plot_result, x=self._internal('xout'), y=viol, names=invar_names, latex_names=kwargs.pop('latex_names', invar_names), indices=None, **kwargs) def extend_by_integration(self, xend, params=None, odesys=None, autonomous=None, npoints=1, **kwargs): @@ -309,6 +331,9 @@ def extend_by_integration(self, xend, params=None, odesys=None, autonomous=None, new_info[k] = [new_info[k]] new_info[k].append(v) else: - new_info[k] += v + try: + new_info[k] += v + except TypeError: + new_info[k] = v self.info = new_info return self diff --git a/pyodesys/symbolic.py b/pyodesys/symbolic.py index 105d69d..710a0c0 100644 --- a/pyodesys/symbolic.py +++ b/pyodesys/symbolic.py @@ -232,7 +232,7 @@ def __init__(self, dep_exprs, indep=None, params=None, jac=True, dfdx=True, jtimes=False, first_step_expr=None, roots=None, backend=None, lower_bounds=None, upper_bounds=None, linear_invariants=None, nonlinear_invariants=None, linear_invariant_names=None, nonlinear_invariant_names=None, steady_state_root=False, - init_indep=None, init_dep=None, sparse=False, **kwargs): + init_indep=None, init_dep=None, sparse=False, Lambdify_kw=None, **kwargs): self.dep, self.exprs = zip(*dep_exprs.items()) if isinstance(dep_exprs, dict) else zip(*dep_exprs) self.indep = indep if params is True or params is None: @@ -248,6 +248,7 @@ def __init__(self, dep_exprs, indep=None, params=None, jac=True, dfdx=True, self._dfdx = dfdx self.first_step_expr = first_step_expr self.be = Backend(backend) + self.Lambdify_kw = Lambdify_kw if steady_state_root: if steady_state_root is True: @@ -541,7 +542,7 @@ def from_other_new_params(cls, ori, par_subs, new_pars, new_par_names=None, drop_idxs = [ori.params.index(par) for par in par_subs] params = _skip(drop_idxs, ori.params, False) + list(new_pars) back_substitute = _Callback(ori.indep, ori.dep, params, list(par_subs.values()), - Lambdify=ori.be.Lambdify) + Lambdify=ori.be.Lambdify, Lambdify_kw=ori.Lambdify_kw) def recalc_params(t, y, p): rev = back_substitute(t, y, p) @@ -684,7 +685,8 @@ def get_dfdx(self): return self._dfdx def _callback_factory(self, exprs): - return _Callback(self.indep, self.dep, self.params, exprs, Lambdify=self.be.Lambdify) + return _Callback(self.indep, self.dep, self.params, exprs, + Lambdify=self.be.Lambdify, Lambdify_kw=self.Lambdify_kw) def get_f_ty_callback(self): """ Generates a callback for evaluating ``self.exprs``. """ @@ -738,8 +740,8 @@ def get_jtimes_callback(self): if jtimes is False: return None v, jtimes_exprs = jtimes - return _Callback(self.indep, tuple(self.dep) + tuple(v), self.params, - jtimes_exprs, Lambdify=self.be.Lambdify) + return _Callback(self.indep, tuple(self.dep) + tuple(v), self.params, jtimes_exprs, + Lambdify=self.be.Lambdify, Lambdify_kw=self.Lambdify_kw) def get_first_step_callback(self): if self.first_step_expr is None: @@ -1336,7 +1338,8 @@ def analytic_factory(x0, y0, p0, be): @staticmethod def _get_analytic_callback(ori_sys, analytic_exprs, new_dep, new_params): - return _Callback(ori_sys.indep, new_dep, new_params, analytic_exprs, Lambdify=ori_sys.be.Lambdify) + return _Callback(ori_sys.indep, new_dep, new_params, analytic_exprs, + Lambdify=ori_sys.be.Lambdify, Lambdify_kw=ori_sys.Lambdify_kw) def __getitem__(self, key): ori_dep = self.original_dep[self.names.index(key)] diff --git a/pyodesys/tests/test_robertson.py b/pyodesys/tests/test_robertson.py index 9bda692..1f6cf15 100644 --- a/pyodesys/tests/test_robertson.py +++ b/pyodesys/tests/test_robertson.py @@ -16,12 +16,8 @@ @requires('sym', 'sympy', 'pyodeint') -@pytest.mark.xfail(reason="internal numerical check in Boost-ublas lu routine triggers here.") def test_run_integration(): - try: - xout, yout, info = run_integration(integrator='odeint')[:3] - except Exception: # flaky test 'internal logic' from boost's ublas lib (?) - xout, yout, info = run_integration(integrator='odeint', atol=1e-4, rtol=1e-4)[:3] + xout, yout, info = run_integration(integrator='odeint')[:3] assert info['success'] is True @@ -29,7 +25,7 @@ def test_run_integration(): @pycvodes_double def test_run_integration__atol_dict(): xout, yout, info = run_integration( - integrator='cvode', atol={'A': 1e-9, 'B': 1e-10, 'C': 1e-5}, nsteps=1173*3)[:3] + integrator='cvode', atol={'A': 1e-10, 'B': 1e-11, 'C': 1e-6}, nsteps=5000)[:3] assert info['success'] is True @@ -37,7 +33,7 @@ def test_run_integration__atol_dict(): @pycvodes_double def test_run_integration__atol_list(): xout, yout, info = run_integration( - integrator='cvode', atol=[1e-9, 1e-10, 1e-5], nsteps=2500*3)[:3] + integrator='cvode', atol=[1e-10, 1e-11, 1e-6], nsteps=5000)[:3] assert info['success'] is True @@ -140,48 +136,29 @@ def test_get_ode_exprs_ODESys(): atol=1e-8, rtol=1e-12, extra_forgive=2) _test_goe(symbolic=False, logc=False, logt=True, zero_conc=0, zero_time=1e-12, atol=1e-8, rtol=1e-12, extra_forgive=0.4) + for reduced in range(4): + _test_goe(symbolic=False, reduced=reduced, + extra_forgive=5 if reduced == 2 else 3) + if reduced != 2: + _test_goe(symbolic=False, reduced=reduced, logc=True, logt=False, zero_conc=1e-18, + atol=1e-10, rtol=1e-10, extra_forgive=20, first_step=1e-14, nsteps=17000) + if reduced == 3: + _test_goe(symbolic=False, reduced=reduced, logc=True, logt=True, zero_conc=1e-18, zero_time=1e-12, + atol=1e-12, rtol=5e-13, extra_forgive=1e-3, first_step=1e-13) # note extra_forgive -@requires('sym', 'sympy', 'pycvodes') -@pycvodes_double -@pytest.mark.parametrize("reduced", [0,1,2,3]) -def test_get_ode_exprs_ODESys__reduced__1(reduced): - _test_goe(symbolic=False, reduced=reduced, extra_forgive=6) - -@requires('sym', 'sympy', 'pycvodes') -@pycvodes_double -@pytest.mark.parametrize("reduced", [0,1,3]) -def test_get_ode_exprs_ODESys__reduced__2(reduced): - _test_goe(symbolic=False, reduced=reduced, logc=True, logt=False, zero_conc=1e-18, - atol=1e-10, rtol=1e-10, extra_forgive=20, first_step=1e-14, nsteps=17000) - - -@requires('sym', 'sympy', 'pycvodes') -@pycvodes_double -def test_get_ode_exprs_ODESys__reduced__3(): - reduced = 3 - _test_goe(symbolic=False, reduced=reduced, logc=True, logt=True, zero_conc=1e-18, zero_time=1e-12, - atol=1e-12, rtol=5e-13, extra_forgive=1e-3, first_step=1e-13) # note extra_forgive - + _test_goe(symbolic=False, reduced=reduced, logc=False, logt=True, zero_time=1e-9, + atol=1e-13, rtol=1e-14, first_step=1e-14, extra_forgive=3) -@requires('sym', 'sympy', 'pycvodes') -@pycvodes_double -@pytest.mark.parametrize("reduced", [0,1,2,3]) -def test_get_ode_exprs_ODESys__reduced__4(reduced): - _test_goe(symbolic=False, reduced=reduced, logc=False, logt=True, zero_time=1e-9, - atol=1e-13, rtol=1e-14, first_step=1e-14, extra_forgive=3) +@pytest.mark.xfail(reason="flaky due to rounding errors") @requires('sym', 'sympy', 'pycvodes') @pycvodes_double -@pytest.mark.parametrize("reduced", [ - 0, - pytest.param(1, marks=pytest.mark.xfail(reason="flaky, CV_REPTD_RHSFUNC_ERR on CI-server")), - 2, - 3 -]) -def test_get_ode_exprs_ODESys__reduced__5(reduced): - _test_goe(symbolic=False, reduced=reduced, logc=False, logt=True, zero_time=1e-12, - atol=1e-9, rtol=1e-10, extra_forgive=8, nonnegative=True) # tests RecoverableError - +def test_get_ode_exprs_ODESys__flaky(): + for reduced in range(4): + _test_goe(symbolic=False, reduced=reduced, logc=False, logt=True, zero_time=1e-12, + atol=1e-8, rtol=1e-10, extra_forgive=1, nonnegative=True, + first_step=1e-14) # tests RecoverableError + @requires('sym', 'sympy', 'pycvodes') @pycvodes_double diff --git a/pyodesys/tests/test_symbolic.py b/pyodesys/tests/test_symbolic.py index 2bc7dc9..6f23fd7 100644 --- a/pyodesys/tests/test_symbolic.py +++ b/pyodesys/tests/test_symbolic.py @@ -542,7 +542,7 @@ def callback(*args, **kwargs): name='vode', method='bdf', first_step=1e-10) assert info['njev'] > 0 min_time_band = min(min_time_band, time_band) - check(yout_dens[-1, :], n, p, a, atol, rtol, 1.5) + check(yout_dens[-1, :], n, p, a, atol, rtol, 4) check(yout_band[-1, :], n, p, a, atol, rtol, 1.5) assert min_time_dens*2 > min_time_band # (2x: fails sometimes due to load) diff --git a/pyodesys/util.py b/pyodesys/util.py index 6e087f0..b66fa4f 100644 --- a/pyodesys/util.py +++ b/pyodesys/util.py @@ -37,7 +37,7 @@ def _parse_requirement(req): return req.strip(), None, None -class requires: +class requires(object): """Conditional skipping (on requirements) of tests in pytest Examples @@ -184,7 +184,7 @@ def transform_exprs_indep(fw, bw, dep_exprs, indep, check=True): return [(e/fw.diff(indep)).subs(indep, bw) for e in exprs] -class _Blessed: +class _Blessed(object): pass @@ -281,7 +281,7 @@ def pycvodes_klu(cb): reason="Sparse jacobian tests require pycvodes and sundials with KLU enabled.")(cb) -class MissingImport: +class MissingImport(object): def __init__(self, modname, exc): self._modname = modname diff --git a/scripts/render_notebooks.sh b/scripts/render_notebooks.sh index 681c5ea..06f9f5e 100755 --- a/scripts/render_notebooks.sh +++ b/scripts/render_notebooks.sh @@ -17,8 +17,9 @@ function quiet_unless_fail { cd examples/ -PREC=$(python -c "from pycvodes import config; print(config.get('SUNDIALS_PRECISION', 'double'))") +PREC=`python3 -c "from pycvodes import config; print(config.get('SUNDIALS_PRECISION', 'double'))"` set -x +set +e for ipynb in *.ipynb; do if [[ $ipynb == "_native_standalone.ipynb" ]]; then continue # issue with boost's program options @@ -29,15 +30,14 @@ for ipynb in *.ipynb; do if [[ $ipynb == "_bench_native_odesys_multi.ipynb" ]]; then continue fi - #quiet_unless_fail - python -m jupyter nbconvert \ - --log-level=INFO \ - --to=html \ - --ExecutePreprocessor.enabled=True \ - --ExecutePreprocessor.timeout=900 \ - "${ipynb}" + if [[ $ipynb == "_const_step_size.ipynb" ]]; then + continue + fi + quiet_unless_fail jupyter nbconvert --log-level=INFO --to=html --ExecutePreprocessor.enabled=True --ExecutePreprocessor.timeout=900 "${ipynb}" \ + | grep -v -e "^\[NbConvertApp\] content: {'data':.*'image/png'" #if [ ${QUIET_EXIT_CODE} -ne 0 ]; then # exit ${QUIET_EXIT_CODE} #fi done +set -e #../scripts/render_index.sh *.html diff --git a/setup.py b/setup.py index 1f82102..15fb097 100755 --- a/setup.py +++ b/setup.py @@ -60,6 +60,7 @@ def _path_under_setup(*args): submodules = [ 'pyodesys.native', + 'pyodesys.native.symcse', ] tests = [ @@ -77,8 +78,8 @@ def _path_under_setup(*args): _author, _author_email = open(_path_under_setup('AUTHORS'), 'rt').readline().split('<') extras_req = { - 'integrators': ['pyodeint>=0.10.3', 'pycvodes>=0.13.1', 'pygslodeiv2>=0.9.3'], - 'native': ['pycompilation>=0.4.9', 'pycodeexport>=0.1.2'], + 'integrators': ['pyodeint>=0.10.10', 'pycvodes>=0.15.0', 'pygslodeiv2>=0.9.8'], + 'native': ['pycompilation>=0.4.9', 'pycodeexport>=0.1.2', 'appdirs'], 'docs': ['Sphinx', 'sphinx_rtd_theme', 'numpydoc'], 'testing': ['pytest', 'pytest-cov', 'pytest-flakes', 'pytest-pep8', 'rstcheck'] } @@ -96,7 +97,7 @@ def _path_under_setup(*args): license=license, packages=[pkg_name] + submodules + tests, include_package_data=True, - install_requires=['numpy>=1.16.4', 'scipy>=1.2.3', 'sym>=0.3.7', 'appdirs', + install_requires=['numpy>=1.16.4', 'scipy>=1.2.3', 'sym>=0.3.4', 'sympy>=1.5.1', 'matplotlib>=2.2.5', 'notebook' + ('>=5.7.8' if 'CI' in os.environ else '')], tests_require=['pytest>=5.4.1'], extras_require=extras_req,