Merge branch 'master' into trapping-plumes

.github/workflows/draft-pdf.yml

@@ -1,4 +1,10 @@
-on: [push]
+on:
+  push:
+    branches:
+      - master
+  pull_request:
+    branches:
+      - master
 
 jobs:
   paper:

.github/workflows/main.yml

@@ -40,7 +40,7 @@ jobs:
     - name: Install dependencies
       run: |
         pip install --upgrade pip
-        pip install .[dev,miosr,cvxpy]
+        pip install .[dev,miosr,cvxpy,sbr]
     - name: Build the docs
       # Not exactly how RTD does it, but close.
       run: |

.github/workflows/notebooks.yml

@@ -0,0 +1,39 @@
+on:
+  schedule:
+    # run twice a month during times when hopefully few other jobs are scheduled
+    - cron: '0 12 6,21 * *'
+
+jobs:
+  find-notebooks:
+    runs-on: ubuntu-latest
+    outputs:
+      paths: ${{ steps.find-notebooks.outputs.paths }}
+    steps:
+      - name: List Files
+        id: find-notebooks
+        uses: mirko-felice/list-files-action@v3.0.5
+        with:
+          repo: ${{ github.repository }}
+          ref: ${{ github.ref }}
+          path: "examples"
+          ext: ".ipynb"
+
+  run-notebooks:
+    needs: find-notebooks
+    runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+      matrix:
+        files: ${{ fromJson(needs.find-notebooks.outputs.paths) }}
+    steps:
+      - uses: actions/checkout@v3
+      - name: Set up Python
+        uses: actions/setup-python@v3
+        with:
+          python-version: "3.9"
+      - name: Install dependencies
+        run: |
+          pip install .[cvxpy,miosr] sympy nbconvert jupyter matplotlib seaborn pandas dysts
+      - name: Run Notebook
+        run: |
+          jupyter nbconvert --execute --to notebook  --inplace ${{ matrix.files }}

.gitignore

@@ -16,6 +16,7 @@ venv/
 ENV/
 env.bak/
 venv.bak/
+env8
 
 # automatically generated by setuptools-scm
 pysindy/version.py

examples/13_ensembling.ipynb

@@ -42,7 +42,16 @@
      "start_time": "2020-10-22T22:20:40.308783Z"
     }
    },
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/var/folders/37/5qnc4jls0l90nkqjm2xxg1bm0000gp/T/ipykernel_80870/1554735471.py:10: DeprecationWarning: Please use `ODEintWarning` from the `scipy.integrate` namespace, the `scipy.integrate.odepack` namespace is deprecated.\n",
+      "  from scipy.integrate.odepack import ODEintWarning\n"
+     ]
+    }
+   ],
    "source": [
     "import matplotlib.pyplot as plt\n",
     "from mpl_toolkits.mplot3d import Axes3D\n",
@@ -86,12 +95,14 @@
    "metadata": {},
    "outputs": [
     {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "(x)' = 0.018 1 + -10.024 x + 10.013 y\n",
-      "(y)' = 27.808 x + -0.930 y + -0.996 x z\n",
-      "(z)' = 0.169 1 + -2.669 z + 0.998 x y\n"
+     "ename": "TypeError",
+     "evalue": "SINDy.fit() got an unexpected keyword argument 'ensemble'",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[2], line 25\u001b[0m\n\u001b[1;32m     23\u001b[0m model \u001b[38;5;241m=\u001b[39m ps\u001b[38;5;241m.\u001b[39mSINDy(feature_names\u001b[38;5;241m=\u001b[39mfeature_names, optimizer\u001b[38;5;241m=\u001b[39mensemble_optimizer)\n\u001b[1;32m     24\u001b[0m \u001b[38;5;66;03m# Ensemble with replacement (V1)\u001b[39;00m\n\u001b[0;32m---> 25\u001b[0m model\u001b[38;5;241m.\u001b[39mfit(x_train, t\u001b[38;5;241m=\u001b[39mdt, ensemble\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m, quiet\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n\u001b[1;32m     26\u001b[0m model\u001b[38;5;241m.\u001b[39mprint()\n\u001b[1;32m     27\u001b[0m ensemble_coefs \u001b[38;5;241m=\u001b[39m model\u001b[38;5;241m.\u001b[39mcoef_list\n",
+      "\u001b[0;31mTypeError\u001b[0m: SINDy.fit() got an unexpected keyword argument 'ensemble'"
      ]
     }
    ],
@@ -1079,7 +1090,7 @@
    "hash": "3ee6f1cb9fc3b265a5f24cdb7fa225f31e54d7494aa3be0e32b8f891af359708"
   },
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -1093,7 +1104,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.4"
+   "version": "3.11.5"
   },
   "toc": {
    "base_numbering": 1,

examples/19_bayesian_sindy/example.py

@@ -0,0 +1,197 @@
+#!/usr/bin/env python
+# coding: utf-8
+# # Bayesian UQ-SINDy
+# In[1]:
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.integrate import solve_ivp
+
+import pysindy as ps
+
+
+# In[2]:
+
+
+# set seed for reproducibility
+np.random.seed(987)
+
+
+# ### Lotka-Volterra Predator-Prey Model
+#
+# In this example, we generate the data using the Lotka-Volterra equations, which is a simplified model of Predator-Prey interactions. They specify a system of Ordinary Differential Equations (ODEs):
+# \begin{align}
+# \frac{dP}{dt} &= a P - b P Q\\
+# \frac{dQ}{dt} &= c P Q - d Q
+# \end{align}
+# where $P$ is the concentration of prey, $Q$ is the concentration of predators, $a$ is the birth rate of prey, $b$ is the death rate of prey, $c$ is the birth rate of predators and $d$ is the death rate of predators.
+#
+# For more details, see e.g. Rockwood L. L. and Witt J. W. (2015). *Introduction to population ecology*. Wiley Blackwell, Chichester, West Sussex, UK, 2nd edition
+
+# In[3]:
+
+
+# set up a class that defines the Lotka-Volterra equations
+class PredatorPreyModel:
+    def __init__(self, a=1.0, b=0.3, c=0.2, d=1.0):
+        # internalise the model parameters.
+        self.a = a
+        self.b = b
+        self.c = c
+        self.d = d
+
+    def dydx(self, t, y):
+        # Lotka-Volterra Model model, see e.g. https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations.}
+        return np.array(
+            [self.a * y[0] - self.b * y[0] * y[1], self.c * y[0] * y[1] - self.d * y[1]]
+        )
+
+    def solve(self, **kwargs):
+        # solve the system of ODEs.
+        return solve_ivp(self.dydx, **kwargs)
+
+
+# In[4]:
+
+
+# set some hyperparameters.
+t_span = [0, 30]
+y0 = np.array([10, 5])
+max_step = 0.1
+
+# initialise the model and solve.
+my_model = PredatorPreyModel()
+sol = my_model.solve(t_span=t_span, y0=y0, max_step=max_step)
+
+
+# In[5]:
+
+
+# the noise level.
+noise = 0.1
+
+# extract the timesteps and perturb the solution with noise.
+t = sol.t
+P = sol.y[0, :] + np.random.normal(scale=noise, size=sol.t.size)
+Q = sol.y[1, :] + np.random.normal(scale=noise, size=sol.t.size)
+
+
+# In[6]:
+
+
+# plot the solution.
+plt.figure(figsize=(12, 4))
+plt.plot(sol.t, sol.y[0, :], label="Prey")
+plt.scatter(t, P)
+plt.plot(sol.t, sol.y[1, :], label="Predators")
+plt.scatter(t, Q)
+plt.legend()
+plt.show()
+
+
+# ### Bayesian UQ-SINDy
+#
+# Here we recover the governing equations using UQ-SINDy. For more details on the theory of the method, see Hirsh, S. M., Barajas-Solano, D. A., & Kutz, J. N. (2021). *Sparsifying Priors for Bayesian Uncertainty Quantification in Model Discovery* (arXiv:2107.02107). arXiv. http://arxiv.org/abs/2107.02107
+#
+# Note that the current implementation differs from the method described in Hirsh et al. (2021) by imposing the error model directly on the derivatives, rather than on the states, circumventing the need to integrate the equation to evaluate the posterior density. One consequence of this is that the noise standard deviation "sigma" is with respect to the derivatives instead of the states and hence should not be interpreted.
+#
+# The underlying code used to find the posterior distribution of model parameters is `numpyro.infer.MCMC` using the `numpyro.infer.NUTS` kernel. Note that all keyword arguments passed to `pysindy.optimizers.SBR` are sent forward to the [MCMC sampler](https://num.pyro.ai/en/stable/mcmc.html).
+
+# In[7]:
+
+
+# set sampler hyperparameters
+sampling_seed = 123
+
+if __name__ == "testing":
+    num_warmup = 10
+    num_samples = 100
+    num_chains = 1
+else:
+    num_warmup = 500
+    num_samples = 2000
+    num_chains = 2
+
+
+# In[8]:
+
+
+# initialise the Sparse bayesian Regression optimizer.
+optimizer = ps.optimizers.SBR(
+    num_warmup=num_warmup,
+    num_samples=num_samples,
+    mcmc_kwargs={"seed": sampling_seed, "num_chains": num_chains},
+)
+
+# use the standard polynomial features.
+feature_library = ps.feature_library.polynomial_library.PolynomialLibrary(
+    include_interaction=True
+)
+
+# initialise SINDy and fit to the data.
+sindy = ps.SINDy(optimizer, feature_library, feature_names=["P", "Q"])
+sindy.fit(np.c_[P, Q], t=t)
+
+
+# In[9]:
+
+
+# set up a new differential equation that uses the Bayesian SINDy predictions.
+def surrogate_dydt(t, y):
+    _y = y[np.newaxis, :]
+    return sindy.predict(x=_y)
+
+
+# solve using the Bayesian SINDy equations.
+surrogate_sol = solve_ivp(surrogate_dydt, t_span=t_span, y0=y0, max_step=max_step)
+
+
+# In[10]:
+
+
+# plot the surrogate solution.
+plt.figure(figsize=(12, 4))
+
+plt.plot(surrogate_sol.t, surrogate_sol.y[0, :], label="Prey")
+plt.scatter(t, P)
+
+plt.plot(surrogate_sol.t, surrogate_sol.y[1, :], label="Predators")
+plt.scatter(t, Q)
+
+plt.legend()
+plt.show()
+
+
+# ### Get MCMC diagnostics
+#
+# We can inspect the posterior samples in more detail using `arviz`. Note that this is not included as a dependency of `pysindy` and must be installed separately.
+
+# In[11]:
+
+
+# import arviz.
+import arviz as az
+
+
+# In[12]:
+
+
+# convert the numpyro samples to an arviz.InferenceData object.
+samples = az.from_numpyro(sindy.optimizer.mcmc_)
+
+# have a look at the summray.
+az.summary(samples)
+
+
+# In[13]:
+
+
+# plot the traces.
+az.plot_trace(samples, divergences=False)
+plt.tight_layout()
+plt.plot()
+
+
+# In[ ]:
+
+
+# In[ ]: