Gaussian quadrature #141
Gaussian quadrature #141
Conversation
|
I changed PR base as new features go into develop first 🙏 |
|
Just request a review when I should review this! |
|
Hi, I can't seem to be able to request a review through the usual method. However, please feel free to take a look at your convenience. All results for 1-dimensional integrals are correct now. You can see the tests that I did in this Colab. Probably the integration does not work correctly for higher dimensions. When I first looked at this I was interested in vectorized integration - integrating over vectors of limits in one dimension - which is just the result of treating dims as the number of limits and not summing over all dimensions in the end. This is what quadpy does as well, and since it turns out that doing fixed-point and iterative integration in Torch is very slow even with large numbers of points, I won't be looking into it much more, or into high-dimension integration at all. (Also with the particular integral I wanted to do quickly, it looks like tensor-based backends won't help either since it is bogged down by a high-order cross-section calculation.) Hopefully someone who has a better knowledge of how high-dimensional integration should work can quickly get it working based on the start that's been made here. |
|
Yay! |
|
I do not know what is going on with that test. |
I had this problem in another repo as well, it's only the posting of the coverage comment, not sure if the thing I implemented is quite as robust as I had hoped 🥲 Will investigate soon if it continues. In any event not related to this PR, I think. |
ilan-gold
force-pushed
the
gaussian-quadrature
branch
from
5b93185
to
c31ccae
Compare
| return f | ||
|
|
||
| def _weights(self, N, dim, backend, requires_grad=False): | ||
| return None |
There was a problem hiding this comment.
This should never be called, right? If so, then we might as well raise an exception in case it somehow does?
There was a problem hiding this comment.
I made this an implemented method in the base class because the Gaussian integral inherits from this and I've implemented everything here around generating the grid etc. which requires this method.
There was a problem hiding this comment.
Sure, I am just thinking practically speaking this function should always be overwritten by any derived class, right? And if not, then it at least shouldn't be called. So instead of returning None and potentially failing silently if either it wasn't overwritten or called by accident we should raise an exception?
There was a problem hiding this comment.
The idea here is that some integration methods have weights and some do not, but the GridIntegrator should be able to handle both. So if you prefer I can make the base class method abstract here instead of having this default behavior and then add this implementation to the NewtonCotes class instead tor return None. Let me know :)
There was a problem hiding this comment.
The idea here is that some integration methods have weights and some do not, but the
GridIntegratorshould be able to handle both. So if you prefer I can make the base class method abstract here instead of having this default behavior and then add this implementation to theNewtonCotesclass instead tor returnNone. Let me know :)
Yes, let's do that then :)
|
@gomezzz I generally cleaned up the tests. I also addressed the other comments and left the ones unresolved where we might have a discussion. |
|
@ilan-gold I am trying to fix the CI in #169 but on that PR it works. Could you try to merge the branch into this one to see if I fixed it? Otherwise, my suspicion is that it maybe because the PR comes from a fork? There seems to be some permission problems as described here MishaKav/pytest-coverage-comment#68 (the mentioned repo setting is, however, correct, but maybe the one in the fork also affects this? Not quite sure yet what is going on :) ) |
…ussian_quadrature
|
A class MC failure @gomezzz - I will re-run but posting for posterity: https://github.com/esa/torchquad/actions/runs/4530799642/jobs/7980177519?pr=141 Not really sure what to do about these MC failures in the single-dimension integral domain case. It seems like they are generally against the concept of unit-testing (since they are random). Maybe we could set a seed? |
I think they are already seeded, but not all the frameworks provide deterministic seeding on different machines, iirc? I am afraid we may just need to increase the bounds instead 🤷 |
|
@gomezzz I increased the bounds. I think everything is passing now. I have one more small PR based on https://github.com/ilan-gold/torchquad/tree/special_grid to make - basically allowing the user to disable a check and then generalizing the usage of |
|
don't know what to say here @gomezzz about this coverage business... |
|
So we are good to merge @gomezzz? |
Add Gaussian quadrature methods
Base Gaussian quadrature class created
Several Gaussian quadrature methods implemented, for integration over [a,b], [0,inf] and [-inf,inf]
Minor changes to BaseIntegrator to allow for args for input functions
Resolved Issues
How Has This Been Tested?
How Has This Not Been Tested?
All except Gauss-legendre have not been tested on the GPU. Most have also not been tested in multiple dimensions.
No backends besides Torch and numpy have been tested
Known issues
Finding the roots of some polynomials (Jacobi, Laguerre) use scipy functions. I don't know if/how this is compatible with autoray.do()
Transforming points and weights for integration can be tricky with the GPU involved - need to be able to check/enforce that all tensors are on the same device
Gauss-Jacobi integration does not currently give the correct solution outside the interval [-1,1]
Not yet implemented
other Gaussian methods including Gauss-Kronrod / QUADPACK