Factor Gaussian to improve stability

[iangrooms]

This PR gives users the option of splitting a Gaussian filter into a sequence of Gaussian filters each of which has a smaller length scale. Each smaller-scale Gaussian filter is more stable to roundoff. For the true/target filter, and in the absence of roundoff errors, there is no difference between applying one Gaussian filter or applying several small-scale filters. Even without roundoff errors, it is not exactly equivalent to filter once versus more than once because of the polynomial approximation to the Gaussian filter. But because our defaults give a very high accuracy in the polynomial approximation, the two methods give extremely similar results, and the factored method can be more stable to roundoff.

This PR is marked as Draft because @NoraLoose is going to this method it to the tutorial on numerical stability.

[codecov-commenter]

Codecov Report

Merging #127 (f339f06) into master (75f172c) will increase coverage by 0.08%. The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #127      +/-   ##
==========================================
+ Coverage   98.41%   98.49%   +0.08%     
==========================================
  Files           9        9              
  Lines         944      997      +53     
==========================================
+ Hits          929      982      +53     
  Misses         15       15              

Flag	Coverage Δ
unittests	98.49% <100.00%> (+0.08%)	:arrow_up:

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
gcm_filters/filter.py	99.03% <100.00%> (+0.04%)	:arrow_up:
tests/test_filter.py	99.55% <100.00%> (+0.10%)	:arrow_up:

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 75f172c...f339f06. Read the comment docs.

[NoraLoose]

Thanks @iangrooms for implementing this. I think it may be better if I put in a separate PR for showing this feature in the numerical instability notebook. I would like to avoid conflicts with PR #115, which also changes the same notebook in the documentation. The latter PR is essentially ready to merge but needs PR #123 to be reviewed and merged first. So let's not complicate things further.

If you change this PR from "Draft" to "Ready for review", I can have a look.

[NoraLoose]

Thanks @iangrooms for this new feature. The option to break large filters into a sequence of smaller filters is really exciting because this may bail gcm-filters out of the numerical instability issue. 🎉

I left some minor comments above. A major comment is that we should add a test in filter.py that checks on some test data that specifying filter_scale=L and n_iterations=1 gives indeed (approximately) the same as specifying filter_scale=L/sqrt(2) and n_iterations=2.

[rabernat]

This is awesome! Thanks so much Ian. Sounds like this approach could really help out with numerical stability. I agree with Nora that tests are a must.

I wonder if there is a generalization of this splitting approach that might also work with the taper filter.

[iangrooms]

Replying to @rabernat's comment about doing something like this for the Taper filter. I wonder too, but haven't had any good ideas yet. We could apply a sequence of Taper filters that successively zero out more and more small scales. This might improve the stability, because we know that filtering smooth fields tends to be more stable than filtering rough fields. But it also might not help, because the last filter you apply still has a big filter "factor" and could be unstable. The nice thing about the Gaussian version is that you get away with only ever using filters with small filter scales, and I can't figure out how to do that with the Taper.

[iangrooms]

This is great!! The new documentation section is written very clearly, and the new code and test looks good to me.

Final comment / question is whether we want to implement the same test for the viscosity filter, too (for completeness).

The latest commit adds a test of the iterated Gaussian for the vector filter. It uses n_iterations = 4 so it satisfies @NoraLoose's other suggestion too, as well as the minor changes to comments and docs.

[NoraLoose]

The PR looks great to me. I don't have any further comments.

It would be nice to have another reviewer look over this PR before we merge, since this PR makes some major changes to the code. Any volunteers?