Skip to content

Commit

Permalink
Reformat text
Browse files Browse the repository at this point in the history
  • Loading branch information
briandesilva committed Dec 2, 2021
1 parent d694cf7 commit c40d10f
Showing 1 changed file with 30 additions and 22 deletions.
52 changes: 30 additions & 22 deletions examples/12_weakform_SINDy_examples.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -306,7 +306,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"### Performance should improve as the number of sub-domain integrations points increases (numerically approximating the integrals better and better) and number of sub-domains increases (more points for regression). Let's use some noisy Lorenz data and investigate."
"Performance should improve as the number of sub-domain integrations points increases (numerically approximating the integrals better and better) and number of sub-domains increases (more points for regression). Let's use some noisy Lorenz data and investigate."
]
},
{
Expand Down Expand Up @@ -379,7 +379,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"### Performance clearly improves as the number of subdomains and number of integrations points increase! We can also see that Lorenz is correctly identified despite ~20% noise levels\n",
"Performance clearly improves as the number of subdomains and number of integrations points increase! We can also see that Lorenz is correctly identified despite ~20% noise levels.\n",
"\n",
"\n",
"The plot belows shows that we can use the weak-formulation to build models that are robust to noise, and additionally indicates convergence as the regression becomes larger and more accurate. "
]
},
Expand Down Expand Up @@ -415,10 +417,10 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"### The default scikit-learn functionality for model.predict is to return x_dot of the same type as the training data. \n",
"So for weak form, model.predict returns a prediction of the weak form of x_dot, rather than a prediction of x_dot. \n",
"#### We can get around this with a bit of a cheat... inputting the model coefficients from the weak form into the original (not weak) model, and use this to predict!\n",
"#### Beware, this requires: \n",
"The default scikit-learn functionality for `model.predict` is to return x_dot of the same type as the training data. So for weak form, model.predict returns a prediction of the weak form of x_dot, rather than a prediction of x_dot.\n",
"\n",
"We can get around this with a bit of a cheat... inputting the model coefficients from the weak form into the original (not weak) model, and use this to predict!\n",
"Beware, this requires: \n",
"1. That the libraries and library ordering are identical in the two models!\n",
"2. For PDEs, the spatial grids must be identical. This means you need to reuse the library. If you initialize a new PDE library, a new set of subdomains is randomly chosen.\n",
"3. Note that the candidate libraries $\\Theta$ are fundamentally different in the weak and non-weak models. In the former, all the columns are integrated in time (and for PDEs, also in space)! This means if you forecast the weak model coefficients with the non-weak model, you are using a $\\Theta$ matrix that is very noisy! In other words, using the weak form fixed the issues with noise, but forecasting with the original model still has the noise in $\\Theta$.\n",
Expand Down Expand Up @@ -610,7 +612,8 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Now we reformat all the training and testing data and define the PDELibrary\n",
"Now we reformat all the training and testing data and define the PDELibrary.\n",
"\n",
"We also define two different PDELibraries and show that despite identical parameters, they have different subdomain centers, since the centers are randomly chosen within the domain when the PDELibrary is initialized."
]
},
Expand Down Expand Up @@ -736,8 +739,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"#### The takeaway here is that the weak formulation drastically improved our system identification on the Burgers' equation with added noise. \n",
"The weak formulation can handle even more than 20% noise here, although then a fairly large value for K is required to obtain a decent system identification. This is required for some of the examples below, and the downside is that this slows down the code considerably. "
"The takeaway here is that the weak formulation drastically improved our system identification on the Burgers' equation with added noise. The weak formulation can handle even more than 20% noise here, although then a fairly large value for K is required to obtain a decent system identification. This is required for some of the examples below, and the downside is that this slows down the code considerably. "
]
},
{
Expand Down Expand Up @@ -907,12 +909,12 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Okay, so all the optimizers essentially capture the model but often have some much smaller \"extra\" terms. \n",
"**Okay, so all the optimizers essentially capture the model but often have some much smaller \"extra\" terms.**\n",
"The best way to deal with these spurious terms is to use ensembling, i.e. the generation of many models by sub-sampling the data, or sub-sampling the SINDy candidate library. See notebook 13 for many examples of how to use these methods. \n",
"\n",
"You can also deal with this by scanning over the hyperparameters for each method although this is more laborious for the user. \n",
"\n",
"#### Next we try the SR3 optimizer on the same data but with added noise of varying levels to illustrate the robustness to noisy data\n",
"**Next we try the SR3 optimizer on the same data but with added noise of varying levels to illustrate the robustness to noisy data.**\n",
"Ideally, we would cross-validate over 10-20 noise instantiations, but with this high-dimensional data this can be computationally slow. We compute the coefficient model errors defined through\n",
"$$\\Delta\\xi_{u_{xx}} = \\|\\xi_{u_{xx}}^{true} - \\xi_{u_{xx}}^{pred}\\| / \\|\\xi_{u_{xx}}^{true}\\| = \\|-1 - \\xi_{u_{xx}}^{pred}\\|,$$\n",
"and similarly for the other coefficients."
Expand Down Expand Up @@ -980,7 +982,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Now that we have computed the error in the model coefficients at each noise level, we plot the results\n",
"Now that we have computed the error in the model coefficients at each noise level, we plot the results.\n",
"We show below that the weak form nicely works for even 50% Gaussian noise added to every point, showing the power of the weak-formulation for robust system identification."
]
},
Expand Down Expand Up @@ -1064,7 +1066,7 @@
"$$A^2 = u^2 + v^2.$$\n",
"The main change will be a significantly larger library... cubic terms in (u, v) and all their first and second order derivatives. We will also need to generate the data because saving a high-resolution form of the data makes a fairly large file. See the Example 10 Jupyter notebook for the non-weak-form system identification of the reaction-diffusion system.\n",
"\n",
"#### Note that Rudy PDE-FIND paper and Messenger Weak SINDy paper use 256 spatial points in each spatial direction, but Reinbold weak SINDy PRE paper uses 512 points in each direction. We will try and get away with only 64 points in each direction for speed (with normal PDE-FIND this would be a liability because the high order derivatives are very noisy), and still show robustness to ~ 10% noise levels"
"Note that the Rudy PDE-FIND paper and Messenger Weak SINDy paper use 256 spatial points in each spatial direction, but Reinbold weak SINDy PRE paper uses 512 points in each direction. We will try and get away with only 64 points in each direction for speed (with normal PDE-FIND this would be a liability because the high order derivatives are very noisy), and still show robustness to ~ 10% noise levels"
]
},
{
Expand Down Expand Up @@ -1260,7 +1262,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Okay, so we have captured the essential terms in the clean data case... can we repeat with some added noise?"
"Okay, so we have captured the essential terms in the clean data case... can we repeat with some added noise?"
]
},
{
Expand Down Expand Up @@ -1295,7 +1297,8 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Don't get the correct PDE... why is this?\n",
"We don't get the correct PDE... why is this?\n",
"\n",
"It is because the resolution of the weak formulation is not high enough. One needs to increase K (the number of subdomains = size of the regression) and num_pts_per_domain (improve the integration, especially important when data is noisy, and we are using simple trapezoid rules). Unfortunately, this is computationally intensive so we leave this improvement out of this example notebook."
]
},
Expand All @@ -1306,7 +1309,7 @@
"# Test weak form PDE functionality on 3D Reaction-Diffusion system\n",
"Can even use weak-form for 3D PDEs although this is getting very computationally intensive! \n",
"\n",
"#### We illustrate the functionality below but the resolution is comically lower here -- too low here to get a good model -- think of this as an illustration of how you could use this, but higher memory is needed for success!\n",
"We illustrate the functionality below but the resolution is comically lower here -- too low here to get a good model -- think of this as an illustration of how you could use this, but higher memory is needed for success!\n",
"\n",
"We will use a 3D reaction-diffusion equation called the Gray-Scott Equation. We are folllowing the example in Section 3.3.3 of Maddu, S., Cheeseman, B. L., Sbalzarini, I. F., & Müller, C. L. (2019). Stability selection enables robust learning of partial differential equations from limited noisy data. arXiv preprint arXiv:1907.07810., https://arxiv.org/pdf/1907.07810.pdf.\n",
"$$u_t = D_u\\nabla^2 u - uv^2 + 0.014(1-u)$$\n",
Expand Down Expand Up @@ -1514,11 +1517,16 @@
"model = ps.SINDy(feature_library=pde_lib, optimizer=optimizer)\n",
"model.fit(u_flattened, x_dot=u_dot_integral, quiet=True)\n",
"model.print()\n",
"coefs = optimizer.coef_\n",
"\n",
"# Note the fit here is really bad... this is because K and num_pts_per_domain need to be \n",
"# larger for an accurate fit, but we are constrained by the heavy memory usage. Reducing the \n",
"# memory use in the code is a good place for future improvements."
"coefs = optimizer.coef_"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note the fit here is really bad... this is because K and num_pts_per_domain need to be \n",
"larger for an accurate fit, but we are constrained by the heavy memory usage. Reducing the \n",
"memory use in the code is a good place for future improvements."
]
}
],
Expand All @@ -1538,7 +1546,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.4"
"version": "3.6.9"
},
"toc": {
"base_numbering": 1,
Expand Down

0 comments on commit c40d10f

Please sign in to comment.