# SINDy
> Discovering governing equations using the SINDy algorithm

- toc: true 
- hide: false
- branch: master
- search_exclude: false
- badges: true
- comments: true
- categories: [differential equations, machine learning]


```python
# hide
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

from scipy.integrate import odeint

from sklearn.metrics import r2_score, mean_absolute_error
from sklearn.metrics import mean_squared_error, median_absolute_error
```


```python
# hide
def mean_absolute_percentage_error(y_true, y_pred): 
    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100 

def symmetric_mape(y_true, y_pred, eps = 1e-8):
    summ = ((np.abs(y_true) + np.abs(y_pred)) + eps)
    return np.mean(np.abs(y_pred - y_true) / summ) * 100
```


```python
# hide
def print_scores(y_test, y_pred):
    print(f"R2 score: {r2_score(y_test, y_pred)}")
    print(f"MSE score: {mean_squared_error(y_test, y_pred)}")
    print(f"MAE score: {mean_absolute_error(y_test, y_pred)}")
    print(f"Median AE score: {median_absolute_error(y_test, y_pred)}")
    print(f"MAPE score: {mean_absolute_percentage_error(y_test, y_pred)}")
    print(f"SMAPE score: {symmetric_mape(y_test, y_pred)}")
```

Sparse Identification of Nonlinear Dynamical Systems (SINDy) is an algorithm to discover governing dynamical equations for time series ${\bf x}(t)$. The key idea is to construct a differential equation:

$$ \frac{d{\bf x}}{d t} = \Theta({\bf x}^T) {\bf x}$$, 

where the time derivative is computed from the time series data and $\Theta({\bf x})$ is a non-linear basis. 

Steps:
 - Compute time derivative: Perhaps the trickiest part, especially for noisy time series although one can use total variation regularized derivatives for such cases as in [here](https://github.com/stur86/tvregdiff)
 - Choose basis: Some non-linear basis constructed from the time series data, for example, polynomial basis.
 - Apply regularized linear regression: Apply Lasso/Ridge with one step or sequential thresholding least squares.
 
Paper: https://www.pnas.org/content/113/15/3932 

Once the underlying dynamical equations are discovered, **forecasting** becomes a lot easier.

**Extensions**:

Knowing the right coordinates and basis functions to use is often difficult. The most interesting extension of SINDy so far has been to use latent basis / coordinates instead of physical space basis. In latent space, 

$$ {\bf x} \longrightarrow {\bf z} = Encoder({\bf x}).$$

A non-linear basis in ${\bf z}$ is used to perform SINDy. The full cost function also takes into account the physical space SINDy expression. For more details, see the paper:

Paper: https://www.pnas.org/content/116/45/22445

## Example

### Create synthetic dataset


```python
# hide

# **REMEMBER**: change `dx1dt`, `dx2dt` also if you change the ranfun below.
def fun(y, t):
    x1, x2 = y
    dxdt = [-x1**3 - x2, x1 - x2**3]
    return dxdt

t_full = np.linspace(0, 15, 1501)

# The following two functions work great, with alpha = 0.0001, they are being correctly classified.
#y0 = [0.1, 0.05] # LINEAR: small x1, x2 => Linear coupled oscillator
y0 = [0.5, 0.5] # NON-LINEAR: Getting close to triggering non-linear effects
```

**Solve** the ODE:

$$
\begin{align}
\frac{dx_1}{dt} &= -x_1^3 - x_2 \\
\frac{dx_2}{dt} &= x_1 - x_2^3 
\end{align}
$$


```python
# hide
sol = odeint(fun, y0, t_full)
sol_new = sol[:1001]
```

Plot the solutions


```python
# hide_input
t_discover = t_full[:1001]

plt.plot(t_discover, sol_new[:, 0], 'b', label='x1(t)')
plt.plot(t_discover, sol_new[:, 1], 'g', label='x2(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
# plt.savefig('time_series.png', bbox_inches='tight')
```

## Compute the time derivative. 


```python
x1 = sol_new[:,0]
x2 = sol_new[:,1]
```

**Actual** time derivatives


```python
dx1dt = -x1**3 - x2
dx2dt =  x1    - x2**3
```

**Numerically** computed derivatives from data


```python
# hide
dt = t_discover[1] - t_discover[0]
dx1dt_data = np.gradient(x1, dt)
dx2dt_data = np.gradient(x2, dt)
# or using TVRegDiff
# dx1dt = TVRegDiff(x1, 10, 0.1, dx=dt, plotflag=0)
# See https://github.com/stur86/tvregdiff 
# for information on different parameters
```


```python
# hide
print_scores(dx1dt, dx1dt_data)
```

    R2 score: 0.9999999965191242
    MSE score: 2.8269156140981004e-10
    MAE score: 4.206507155221408e-06
    Median AE score: 3.6181352647546294e-06
    MAPE score: 0.004477342875542224
    SMAPE score: 0.0022391799382288345



```python
# hide
print_scores(dx2dt, dx2dt_data)
```

    R2 score: 0.9999996026514059
    MSE score: 2.195307457531667e-08
    MAE score: 9.892101759080555e-06
    Median AE score: 4.016199533329878e-06
    MAPE score: 0.009102822948966595
    SMAPE score: 0.004535741601342113



```python
# hide
plt.plot(dx1dt[-50:], label='Actual', alpha=0.4)
plt.plot(dx1dt_data[-50:], ':', lw=2, label='from_data')

plt.xlabel('Time', fontsize=14)
plt.ylabel(r'$dx_1/dt$', fontsize=14)
plt.legend()
```


```python
# hide
plt.plot(dx2dt[-50:], label='Actual', alpha=0.4)
plt.plot(dx2dt_data[-50:], ':', lw=2, label='from_data')

plt.xlabel('Time', fontsize=14)
plt.ylabel(r'$dx_2/dt$', fontsize=14)
plt.legend()
```


```python
# hide
X = np.zeros((sol_new.shape[0], sol_new.shape[1]))

X[:, 0] = x1 # dx1/dt
X[:, 1] = x2 # dx2/dt
```

## Construct a basis: polynomial (or trig)


```python
# collapse
from sklearn.preprocessing import PolynomialFeatures

dum_data = pd.DataFrame({'x1': x1, 'x2': x2})

deg = 3 # Polynomial degree to use

p = PolynomialFeatures(degree=deg,include_bias=True).fit(dum_data)
xpoly = p.fit_transform(dum_data)
newdf = pd.DataFrame(xpoly, columns = p.get_feature_names(dum_data.columns))

print("Feature names:", list(newdf))#newdf.columns.values.tolist())
print("Feature array shape:", newdf.shape)
```

    Feature names: ['1', 'x1', 'x2', 'x1^2', 'x1 x2', 'x2^2', 'x1^3', 'x1^2 x2', 'x1 x2^2', 'x2^3']
    Feature array shape: (1001, 10)


## Regression using Lasso

(or Ridge/OLS with sequential thresholding)

Lasso does regularized linear regression with L1-norm. **alpha** is a hyperparameter.
 - low alpha -> OLS
 - high alpha -> most features zero


```python
from sklearn.linear_model import Lasso
mod = Lasso(alpha=0.0001)
mod
```




    Lasso(alpha=0.0001)




```python
# hide
# Prepare training data
newdf_train, newdf_test = newdf[:800], newdf[800:]
dx1dt_train, dx1dt_test = dx1dt[:800], dx1dt[800:]
dx2dt_train, dx2dt_test = dx2dt[:800], dx2dt[800:]
```

### $dx_1 / dt$


```python
# hide
mod.fit(newdf_train, dx1dt_train)
print(mod.coef_) # should give the 3rd (x2) + 4th (x1^2) argument non-zero
print(mod.intercept_)
mod.score(newdf_test, dx1dt_test)
```

    [ 0.         -0.05576026 -0.99852418  0.         -0.         -0.01819165
     -0.50202469 -0.          0.         -0.        ]
    0.0033302444886464527





    0.9943287924302203




```python
# hide
fit_dx1 = pd.DataFrame(columns=newdf.columns)
fit_dx1.loc[0] = mod.coef_
fit_dx1.abs().sort_values(by=0, axis=1, ascending=False)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>x2</th>
      <th>x1^3</th>
      <th>x1</th>
      <th>x2^2</th>
      <th>1</th>
      <th>x1^2</th>
      <th>x1 x2</th>
      <th>x1^2 x2</th>
      <th>x1 x2^2</th>
      <th>x2^3</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.998524</td>
      <td>0.502025</td>
      <td>0.05576</td>
      <td>0.018192</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
    </tr>
  </tbody>
</table>
</div>




```python
# hide
ypred_x1 = mod.predict(newdf_test)
print_scores(dx1dt_test, ypred_x1)
```

    R2 score: 0.9943287924302203
    MSE score: 8.443005835967624e-05
    MAE score: 0.009127412469564932
    Median AE score: 0.00942351242559783
    MAPE score: 59.78370313040213
    SMAPE score: 9.034670809972091


Identified ODE:

$dx_1/dt \sim -x_2 - x_1^3$

with minor contributions from other terms that could be gotten rid of by applying a threshold


```python
# Drop features with absolute values less than 0.1
dx1_thr = fit_dx1[fit_dx1.columns[fit_dx1.abs().max() > 0.1]]
dx1_thr
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>x2</th>
      <th>x1^3</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>-0.998524</td>
      <td>-0.502025</td>
    </tr>
  </tbody>
</table>
</div>




```python
# hide_input
# PLOT results
t_test  = np.linspace(8, 10, 201)

#plt.plot(t, sol_new[:, 1], label='Actual', alpha = 0.4)
#plt.plot(ypred, label='Actual')
plt.plot(t_test, ypred_x1, 'r--', lw = 2, label='Prediction')
plt.plot(t_discover, dx1dt, 'g', label='Actual', alpha=0.4)
plt.axvline(x = 8, color='k', linestyle='--')
plt.xlabel('Time', fontsize=14)
plt.ylabel(r'$dx_1/dt$', fontsize=14)
plt.grid()
plt.legend()
# plt.savefig('dx1dt_fit.png', bbox_inches='tight')
```

### $dx_2 / dt$


```python
# hide
mod.fit(newdf_train, dx2dt_train)
print(mod.coef_) # should give the 2nd (x1) + last (x2^3) argument non-zero
print(mod.intercept_)
mod.score(newdf_test, dx2dt_test)
```

    [ 0.          0.99491832 -0.02198579 -0.         -0.         -0.04791986
      0.          0.         -0.         -0.82689872]
    0.0022126622685510813





    0.9965028656692735




```python
# hide
fit_dx2 = pd.DataFrame(columns=newdf.columns)
fit_dx2.loc[0] = mod.coef_
fit_dx2.abs().sort_values(by=0, axis=1, ascending=False) 
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>x1</th>
      <th>x2^3</th>
      <th>x2^2</th>
      <th>x2</th>
      <th>1</th>
      <th>x1^2</th>
      <th>x1 x2</th>
      <th>x1^3</th>
      <th>x1^2 x2</th>
      <th>x1 x2^2</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.994918</td>
      <td>0.826899</td>
      <td>0.04792</td>
      <td>0.021986</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.0</td>
    </tr>
  </tbody>
</table>
</div>




```python
# hide
ypred_x2 = mod.predict(newdf_test)
print_scores(dx2dt_test, ypred_x2)
```

    R2 score: 0.9965028656692735
    MSE score: 1.3128657613492703e-05
    MAE score: 0.003284761866157678
    Median AE score: 0.0036951175861815455
    MAPE score: 1.962931213051045
    SMAPE score: 0.9952953131151899


Identified ODE:

$dx_2/dt \sim x_1 - x_2^3$

with minor contributions from other terms that could be gotten rid of by applying a threshold


```python
# Drop features with absolute values less than 0.1
dx2_thr = fit_dx2[fit_dx2.columns[fit_dx2.abs().max() > 0.1]]
dx2_thr
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>x1</th>
      <th>x2^3</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.994918</td>
      <td>-0.826899</td>
    </tr>
  </tbody>
</table>
</div>




```python
# hide_input
# PLOT results
t_test  = np.linspace(8, 10, 201)

#plt.plot(t, sol_new[:, 0], label='Actual', alpha = 0.4)
#plt.plot(ypred, label='Actual')
plt.plot(t_test, ypred_x2, 'r--', lw = 2, label='Prediction')
plt.plot(t_discover, dx2dt, 'g', label='Actual', alpha = 0.4)
plt.axvline(x = 8, color='k', linestyle='--')
plt.xlabel('Time', fontsize=14)
plt.ylabel(r'$dx_2/dt$', fontsize=14)
plt.grid()
plt.legend()
# plt.savefig('dx2dt_fit.png', bbox_inches='tight')
```

## Forecasting using SINDy

Now that we have discovered the differential equations, we can use the ODE solver to forecast the future.


```python
# Manually entering values for coefficients but this can be automated
def fun_forecast(y, t):
    x1, x2 = y
    dxdt = [-0.5*(x1**3) - 0.9985*x2, 
            0.995*x1     - 0.8269*x2**3]
    return dxdt

t_forecast = np.linspace(10, 15, 500)

y0 = [x1[-1], x2[-1]]

sol_forecast = odeint(fun_forecast, y0, t_forecast)
```


```python
# hide
plt.plot(t_forecast, sol_forecast[:, 0], 'b', label='x1(t)')
plt.plot(t_forecast, sol_forecast[:, 1], 'g', label='x2(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
```


```python
# hide_input
plt.plot(t_forecast, sol_forecast[:, 0], 'b--', label='Forecast')
plt.plot(t_full, sol[:, 0], 'g', label='Actual', alpha=0.4)
plt.axvline(x = 10, color='k', linestyle='--')
plt.legend(loc='best')
plt.xlabel('Time', fontsize=14)
plt.ylabel(r'$x_1$', fontsize=14)
plt.grid()
```


```python
# hide_input
plt.plot(t_forecast, sol_forecast[:, 1], 'b--', label='Forecast')
plt.plot(t_full, sol[:, 1], 'g', label='Actual', alpha=0.4)
plt.axvline(x = 10, color='k', linestyle='--')
plt.legend(loc='best')
plt.xlabel('Time', fontsize=14)
plt.ylabel(r'$x_2$', fontsize=14)
plt.grid()
```


```python
# hide
y_pred1, y_pred2 = sol_forecast[:, 0], sol_forecast[:, 1]
y_test1, y_test2 = sol[1001:, 0],      sol[1001:, 1]
```


```python
# hide
for (y_test, y_pred) in zip([y_test1, y_test2], [y_pred1, y_pred2]):
    print(f"R2 score: {r2_score(y_test, y_pred)}")
    print(f"MSE score: {mean_squared_error(y_test, y_pred)}")
    print(f"MAE score: {mean_absolute_error(y_test, y_pred)}")
    print(f"Median AE score: {median_absolute_error(y_test, y_pred)}")
    print(f"MAPE score: {mean_absolute_percentage_error(y_test, y_pred)}")
    print(f"SMAPE score: {symmetric_mape(y_test, y_pred)}")
```

    R2 score: 0.9986014512322046
    MSE score: 2.9194295629861532e-05
    MAE score: 0.004498499556886438
    Median AE score: 0.004095611629683124
    MAPE score: 8.356195293001417
    SMAPE score: 3.427959986309787
    R2 score: 0.9988787480251635
    MSE score: 2.8082044451838518e-05
    MAE score: 0.004387715903581068
    Median AE score: 0.0032235369101030545
    MAPE score: 8.17055503299457
    SMAPE score: 3.0643658219157226



```python

```

f569a834c02ee2dd9dfe7f64231dde771eb83c07
ipynb
Jupyter Notebook
_notebooks/2019-09-10-sindy-basic.ipynb
fnauman/ds-blog
2620141d4758b3d14e20478e30756905d0b3773e
Apache-2.0
_notebooks/2019-09-10-sindy-basic.ipynb
fnauman/ds-blog
2620141d4758b3d14e20478e30756905d0b3773e
Apache-2.0
2020-10-03T12:26:27.000Z
2021-01-16T06:41:44.000Z
_notebooks/2019-09-10-sindy-basic.ipynb
fnauman/ds-blog
2620141d4758b3d14e20478e30756905d0b3773e
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES