Add GaussianQuadrature2D#476
Conversation
jacklovell
left a comment
jacklovell
left a comment
There was a problem hiding this comment.
Very nice, only minor comments from me. Some of the (code comment related) comments apply equally to the existing GaussianQuadrature class on which this is based, so if you can add some clarifying code comments there too that would be nice but not a deal breaker.
| @@ -17,5 +17,5 @@ | |||
| # under the Licence. | |||
|
|
|||
| from cherab.core.math.integrators.integrators1d cimport Integrator1D, GaussianQuadrature | |||
There was a problem hiding this comment.
For consistency, we should rename GaussianQuadrature to GaussianQuadrature1D to distinguish it from the new 2D class. But I suggest keeping GaussianQuadrature as an alias to maintain backwards compatibility for now, and deprecate it in another release or so (removed in Cherab 2.0).
| @@ -17,4 +17,4 @@ | |||
| # under the Licence. | |||
|
|
|||
| from .integrators1d import Integrator1D, GaussianQuadrature | |||
| # See the Licence for the specific language governing permissions and limitations | ||
| # under the Licence. | ||
|
|
||
| from numpy cimport ndarray |
There was a problem hiding this comment.
This is unnecessary: you don't make use of the typed _x_roots, _x_weights, _y_roots and _y_weights in any tight loops: they're only assigned to per the output of roots_legendre, which goes through the Python layer anyway. So you can cdef these as object instead of ndarray and save a numpy cimport, which will speed up compilation.
All the typed access is (correctly) done through the memoryviews ayway.
| - :math:`x_{\mathrm{lower}}`: Lower limit of integration for the x-dimension. | ||
| - :math:`x_{\mathrm{upper}}`: Upper limit of integration for the x-dimension. | ||
| - :math:`y_{\mathrm{lower}}(x)`: Lower limit of integration for the y-dimension, a function of :math:`x`. | ||
| - :math:`y_{\mathrm{upper}}(x)`: Upper limit of integration for the y-dimension, a function of :math:`x`. |
There was a problem hiding this comment.
Either double-escape \mathrm or (preferably) use a raw string for the docstring (r''') and go with single backslashes for all the Latex code.
|
|
||
| self._rtol = value | ||
|
|
||
| cdef _build_cache(self): |
There was a problem hiding this comment.
This doesn't need to be cdef, a simple def will do.
| :param int x_max_order: Maximum order on Gaussian quadrature in the x dimension. Default is `50`. | ||
| :param int x_min_order: Minimum order on Gaussian quadrature in the x dimension. Default is `1`. | ||
| :param int y_max_order: Maximum order on Gaussian quadrature in the y dimension. Default is `50`. | ||
| :param int y_min_order: Minimum order on Gaussian quadrature in the y dimension. Default is `1`. |
There was a problem hiding this comment.
I think it's worth clarifying here that order is iteratively increased from min_order to max_order until the desired tolerance is reached, and if either x or y reaches max order before tolerance is reached the other one continues to increase order while the first is pegged at its max order. It takes a lot of reading of the code to realise this.
This explanation is also missing from the 1D case so it would be nice to add it there but at least the 1D integrator's code is simple enough that a reader can quickly work it out.
| # cache roots and weights for x dimension | ||
| n = (self._x_max_order + self._x_min_order) * (self._x_max_order - self._x_min_order + 1) // 2 | ||
|
|
||
| self._x_roots = np.zeros(n, dtype=np.float64) |
There was a problem hiding this comment.
These are essentially 1D views of ragged arrays, right? Where each row would be the coefficients for a given order. They would be lower triangular matrices if min_x_order = 1. Maybe put a comment somewhere here mentioning this, as it'll save someone at a later date wasting time trying to work out why you're manually indexing instead of just using 2D - or 3D with [x, y, order] - arrays here (like I might have just done...).
| self._y_roots_mv = self._y_roots | ||
| self._y_weights_mv = self._y_weights | ||
|
|
||
| cpdef double profile_evaluate(self, int n, double x_lower, double x_upper, Function1D y_lower, Function1D y_upper): |
There was a problem hiding this comment.
What's the purpose of this function, and evaluate_overhead below? They're not called anywhere: are they left over from testing?
| int x_ibegin, y_ibegin | ||
| int i, j | ||
| double current_integral, previous_integral, y_contribution, error | ||
| double x, y, x_offset, y_offset, x_slope, y_slope, scale, rtol, y_lower_val, y_upper_val |
| self.assertEqual(quadrature.y_max_order, y_max_order) | ||
| self.assertEqual(quadrature.integrand, Arg2D("x")) | ||
|
|
||
| def test_integrate(self): |
There was a problem hiding this comment.
Add tests for x_max_order > y_max_order and x_max_order < y_max_order to complement the existing x_max_order == x_min_order. This may help avoid future regressions.
2 participants
Adds
GaussianQuadrature2Dclass mentioned in #475.Adds also integrators (1D, 2D) to our documentation.
This PR will be marked as draft until PR #473 is merged.