Hello,
I’ve been using the new scipy.integrate.cubature function, which is very cool! I wish to use it in my projects, but I am required to re-use points, weights and function evaluations (nodes, weights, f_nodes) used in the integration process. Glancing over the documentation, it seems the function does not support this quadrature data extraction.
Looking through the source code on scipy.integrate._rules._base.py, I see that apply_fixed_rule function has all the data there, so it should be possible to hack together something that returns all of the quadrature data.
Perhaps anyone has done this before in some branch/fork? If not, I would also be happy to hear any advice on what’s the best way to go about this 
Thank you!
Hi,
Glad you’re enjoying using cubature! I was its main author earlier this year. If you’re interested in a very hacky solution, then something like the following should work.
cubature accepts an argument called rule which is ordinarily a string like "gk21", "gk15", "genz-malik", and this argument tells cubature the underlying cubature rule to use. However, it’s possible to pass a scipy.integrate._rules.Rule class to this argument (see here and here) which can be used to implement a custom integration rule, although this currently isn’t part of the public API.
A Rule is just a class which implements a method estimate(f, a, b, args) and estimate_error(f, a, b, args), corresponding to calculating an estimate of the integral and an estimate of the error over the region (a, b). Since these two methods are arbitrary, you could implement a custom rule which calculated the estimate of the integral over that region and but also saved the corresponding (nodes, weights, f_nodes) data it used to make the calculation.
I’ve written a small example here which modifies the default 21-node Gauss-Kronrod cubature rule to also save a list of (nodes, weights, f_nodes) for each set of evaluations.
I hope that’s helpful – apologies that it’s not the most elegant solution to the problem. If this were to be supported directly by cubature, I think the best approach would be to make the Rule API public and implement integration rules that save/cache the points, weights and function evaluation data directly.
Let me know if you have any more questions or if I’ve misunderstood the question, I’m happy to help.
1 Like
Thanks for the reply and the example @ollybritton!
I think the implementation through the Rule API makes a lot of sense.
However, before asking the question, I didn’t realize that there could be two ways to interpret my query. You could either extract all the quadrature data:
- generated through the quadrature process.
- needed to reproduce the final integral result.
I am interested in both of these cases, and your example deals with the former unless I am missing something.
For the latter, one needs to take care of heappop of the largest error regions here. I managed to hack together a fork which takes care of the latter case and outputs data into the CubatureRegion object.
In any case, your reply answers my original question
P.S. Since I have this opportunity to chat with the main author, I would like to ask if you think that it would be possible to use this cached quadrature data to restart the integration?
Happy to hear that! With regards to the last question about restarting the integration, unless I am misunderstanding, you could add an additional argument to cubature such as previous_result which took in the returned CubatureResult and would initialise regions, est, err here using the previous results rather than subdividing the initial region.