Closed iangrooms closed 2 years ago
Check out this pull request on
See visual diffs & provide feedback on Jupyter Notebooks.
Powered by ReviewNB
Thanks @iangrooms, this sounds really exciting! :rocket:
It is a bit difficult for me to review this PR because there are no notes on how the new method exactly operates. Ian, it would be helpful if you could share notes (if available) and/or provide more comments in the code. That way, I could try to match up theory and code.
In addition, I am wondering if we can leverage the verification tests that @rabernat set up in https://github.com/ocean-eddy-cpt/gcm-filters/pull/79 (where we would need verifications not only for kernels but for full filter operators). The rationale is the following: This PR swaps out the central filtering algorithm which this package is based on for another algo. Don't we want to test that the new algo leads to the same/similar filtered fields as the old one?
It is a bit difficult for me to review this PR because there are no notes on how the new method exactly operates. Ian, it would be helpful if you could share notes (if available) and/or provide more comments in the code. That way, I could try to match up theory and code.
I think I'll update the Filter Theory section to describe how it works. In exact arithmetic it should produce exactly the same answer as the old method, but the new method is more stable to roundoff errors. Sort of like re-arranging the order of the filter steps.
In addition, I am wondering if we can leverage the verification tests that @rabernat set up in https://github.com/ocean-eddy-cpt/gcm-filters/pull/79 (where we would need verifications not only for kernels but for full filter operators). The rationale is the following: This PR swaps out the central filtering algorithm which this package is based on for another algo. Don't we want to test that the new algo leads to the same/similar filtered fields as the old one?
I probably have to look more closely at this. I guess we could put in some kind of test to compare the result of the new algorithm to the old one, but eventually I think we should just drop the old one and use the new one.
I think I'll update the Filter Theory section to describe how it works.
Sounds great! I will wait with the review until the changes to the filter theory are pushed.
I guess we could put in some kind of test to compare the result of the new algorithm to the old one, but eventually I think we should just drop the old one and use the new one.
Yes, exactly. As outlined in #79, the idea is to have tests like these:
filtered = filter.apply(data)
np.testing.assert_allclose(filtered, expected)
where expected
is a field we save explicitly into the test suite based on the current code.
Merging #149 (0cbaba3) into master (e4177da) will increase coverage by
0.33%
. The diff coverage is95.12%
.
@@ Coverage Diff @@
## master #149 +/- ##
==========================================
+ Coverage 97.98% 98.32% +0.33%
==========================================
Files 9 9
Lines 1044 1014 -30
==========================================
- Hits 1023 997 -26
+ Misses 21 17 -4
Flag | Coverage Δ | |
---|---|---|
unittests | 98.32% <95.12%> (+0.33%) |
:arrow_up: |
Flags with carried forward coverage won't be shown. Click here to find out more.
Impacted Files | Coverage Δ | |
---|---|---|
gcm_filters/filter.py | 98.01% <95.00%> (+1.55%) |
:arrow_up: |
tests/test_filter.py | 99.55% <100.00%> (-0.02%) |
:arrow_down: |
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 e4177da...0cbaba3. Read the comment docs.
I added a section on the theory that underpins the new iterative algorithm for applying the filter to data. @NoraLoose can you take a look?
Looks great @iangrooms! I left some minor comments above.
Two bigger comments are:
Thanks for the code review! In response to the two big comments:
which would be unnecessary if we believe the new algorithm is correct, not to mention being onerous to implement
The point of the verification tests would be that we don't have to "believe". They would catch any bugs in the derivation / implementation of the new algorithm, which we may have overseen by simply looking at the code.
I think it needs its own issue (cf https://github.com/ocean-eddy-cpt/gcm-filters/pull/79) and its own PR.
Sure, we can make the verification tests its own issue and PR. But I think those should be resolved first, before we can merge this PR. I can help with these tests.
I ran the old and new algorithms on the numerical instability notebook. I saved the result of the old algorithm (zeta_2Deg
) then plotted a comparison with the result of the new algorithm. The plot is in the latest commit of the numerical instability notebook. The difference is about 7 orders of magnitude smaller than the field itself. So it looks like the old and new algorithms are giving the same results up to roundoff errors (single precision), at least in situations where the old algorithm is stable. Between the theory, the code implementing the theory, and the comparison of the old and new methods, I hope this is enough to convince that the new algorithm works as advertised.
Hi @iangrooms, we merged #153 so we can finalize this PR too. Here are some final minor things that have to be done / discussed:
v0.3
, it is okay to drop the feature as long as we communicate it.) If we keep it, we should add some sentences at the top of the "Factoring the Gaussian Filter" docs section, to the effect of: "Using this feature is not recommended, and only necessary if you run into numerical instability issues, which you most likely won't."three_term_recurrence
, and run the tests without overwriting the test data. (You won't overwrite anything if you simply run pytest
and don't mess with the GCM_FILTERS_OVERWRITE_TEST_DATA
variable.) If the tests pass, the filter algorithm gives the same result as the previous one.I'll merge the new tests. In the meantime, I do plan to remove the numerical instability notebook as well as any related comments in the docs, and the 'factoring the Gaussian.' It's a big update, so I think it makes sense to bump the version number. The old version will still be available if anyone needs it.
If we remove the "Factoring the Gaussian" section in the docs, we should probably also delete the parameter n_iterations
in the code (+ functionality associated to it). I'm fine with this.
I've removed the numerical instability notebook, removed the factored Gaussian from the docs and code, and merged the latest updates to testing. For some reason building the docs fails (all other tests pass); I noticed that this happened after I pulled in the latest changes from the master branch and before I made my own changes. Not sure how to fix it. Can @NoraLoose or @rabernat take a look?
I can build the docs locally for this branch. Not sure what is going on either.
I am guessing what happened is that some version has changed in the readthedocs environment, since we are not pinning specific versions. It looks a bit like this: https://github.com/pydata/pydata-sphinx-theme/issues/511
I will try to dig deeper. If you want, you can merge this PR and we will fix the docs in a separate PR.
I can also build the docs locally, so if you're both agreed then I think it's OK to merge this PR and fix the docs separately.
From the beginning we've had problems with numerical stability, i.e. roundoff error, that accumulates during the iterative application of the polynomial filter to data. See, for example, #33, #124, and #135, and the things we've tried to improve this, e.g. #67, #130, and #134.
This PR uses exactly the same polynomial approximation of the target filter, but applies the polynomial to the data in a completely different way. The old way was based on factoring the polynomial in terms of its roots and then iteratively applying each factor to the data. The new method is based on an algorithm for evaluating a polynomial using its coordinates in the Chebyshev basis, the algorithm being based on the three-term recurrence for Chebyshev polynomials.
For now I have left the "Numerical Instability" example notebook in, showing that the new code can filter 1/4 degree vorticity data to 10 degrees with no problems. In fact I have not been able to break the new code in any example I've tried including the Taper filter with very large filter factors.