New benchmarking examples

[akaptano]

Hi all,

Lanyue Zhang and I have been working on getting a big benchmark paper + example done, and I think the branch is pretty close to being ready for a merge.

Almost all of the changes are in the examples/16_noise_robustness/ folder. There are three jupyter notebooks there, and a utils.py file for some extra functions. One of the notebooks reproduces all the figures in the paper, another does some visualization and generates plots, and the last performs an example hyper parameter scan using the ensembling and other fancy functionality.

I also changed gitlab -> github in the .pre-commit file, since it looks like the issue was that flake8 moved from gitlab to github very recently (https://jira.mongodb.org/browse/PYTHON-3531). This appears to have resolved the issue. However, there is now a pytest issue related to sphinx and/or ipython (jupyter/nbconvert#528) and I'm not sure how to fix this.

[znicolaou]

This example looks really great to me!

I dug into details a little bit, and pushed a small change just now. These are just suggestions--feel free to revert anything back to the last version if you like. I think it is ready for merge with master when you are.

The main change I made was to speed up the Pareto sweeps by calculating x_dot_test and the library .transform() matrices, which I called mats, just once before the sweep rather than for each value of the hyperparameters. The model predictions are then calculated as x_dot_test_pred = [coef_new[j].dot(mat.T).T for mat in mats], where coef_new = np.array(optimizer.coef_list). This hack reproduces the model.predict functionality but allows us to recycle the mats, which is especially useful in the weak case, since it is relatively expensive to calculate the same weak integrals many times. It may be helpful in practice, since the Pareto front scans can take a while, and the hyperparameter_scan functions in utils.py can be used in other cases.

I also included an option strong_rmse in Pareto_scan_ensembling, which allows the Pareto-optimal AIC calculation to use the "strong" rmse value in conjunction with weak_form=True (for any optimizer algorithm), and verified that it is consistent for the STLSQ case in 16_benchmark_paper.ipynb.

For the MIOSR algorithm, I did encounter "GurobiError: Model too large for size-limited license; visit https://www.gurobi.com/free-trial for a full license." I was able to acquire an academic license following the instructions here: https://www.gurobi.com/features/academic-named-user-license/, and everything works okay after copying the license to the right gurobi directory in my anaconda site-packages. The MIOSR runs do take a long time, and I think they will hog resources if running in parallel. I didn't run the sweeps again or save any results. I'd suggest (if you have the patience) to run all the sweeps one-at-a-time on a cluster machine so the runtime can be estimated without competing for resources.

But everything look really nice to me, and I think it is a really helpful example to include. Nice work!

[znicolaou]

PS: I pushed once more fixing the order of the predicted_coefficient in the weak case. For STLSQ, the average coefficient errors are totally consistent with the previous version, as are the average RMSE errors in the weak case with strong_rmse=True. The average RMSE errors are even a couple orders of magnitude smaller for the weak case with strong_rmse=False (which is the setting that I recommend using). Anyway, I should move on to other things now--let me know if you find any new issues!