Add Gaussian quadrature integration for line shape models by vsnever · Pull Request #378 · cherab/core

vsnever · 2022-08-31T18:29:36Z

This PR fixes #366 and makes the following changes:

adds Integrator1D class for calculating definite integrals of Function1D,
adds GaussianQuadrature integrator,
adds NumericallyIntegrableLineShapeModel class for lineshapes that cannot be analytically integrated over a spectral bin,
replaces half sums in StarkBroadenedLine.add_line() with integration over spectral bins.

…ce value.

…ine to use numerical integration over the spectral bin.

Mateasek

Hi @vsnever , thank you for bringing this PR. I really like it and I think it is a great idea. I added some minor comments to the code, but I have one conceptual idea.

For every lineshape, we have to have the integral value of the power within the bin. For some line shapes there exists an analytical function, some have to be done numerically. Also definite integrals are functions of 2 variables. What if you made Function2D base class for the Integrator1D and then instead of having LineShapeModel and NumericallyIntegrableLineShapeModel classes, there could be again only one class, which would take a Function2D instance as the integrator parameter? In the code it does not matter really, if you are calling an analytically defined Function2D (created with the function framework) of if you call a numerical integrator, also being a Function2D.

This would also need to store the func parameter from integrate method as a attribute of the Integrator1D class, but I think that it is only a minor change.

I hope I managed to express myself clearly. What do you think?

Mateasek · 2022-09-09T14:41:32Z

+
+        i = 0
+        for order in range(self._min_order, self._max_order + 1):
+            self._roots[i:i + order], self._weights[i:i + order] = roots_legendre(order)


in the integrate method you do operation (self._roots_mv[i] + 1.), storing root + 1 value in the array could bring some speedups and it seems that you don't use the value of root elsewhere

I rewrote the integration, so x now taken as x = c + d * self._roots_mv[i], where c = 0.5 * (a + b) and d = 0.5 * (b - a). Adding 1 is no longer required.

Mateasek · 2022-09-09T14:48:55Z

+        for order in range(self._min_order, self._max_order + 1):
+            newval = 0
+            for i in range(ibegin, ibegin + order):
+                x = a + 0.5 * (b - a) * (self._roots_mv[i] + 1.)
+                newval += self._weights_mv[i] * func.evaluate(x)
+            newval *= 0.5 * (b - a)
+


a + 0.5 * (b - a) and 0.5 * (b - a) are constants over the iteration, please consider calculating the values outside of the cycles. If you also apply the comment suggesting to store root + 1 value in self._roots the code would speed up a bit. Considering how many times the integrate method can be called in an observation, any speedup would be I think welcomed (regarding the fact that it is for minimum amount of work)

Mateasek · 2022-09-09T14:51:09Z

+    Compute a definite integral of a one-dimensional function.
+    """
+
+    cpdef (double, double) integrate(self, Function1D func, double a, double b):


please consider renaming the integrate method to evaluate, which is used in the function framework. Also, why not to add the python __call__ method, as is done in the function framework?

Also, is returning the error necessary and what we want? The maximum value of the error can be controlled and the error value is not needed during observations, where you are throwing out the value anyway.

The Integrator1D now inherits from the Function2D. I renamed integrate to evaluate and now this method returns a single value. Now the __call__ method from the parent class can be used.

Mateasek · 2022-09-09T15:05:23Z

+    def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, Integrator1D integrator=None,
+                 dict stark_model_coefficients=None):


Please consider putting integrator parameter at the end for backwards compatibility

…peModel.

vsnever · 2022-09-11T15:45:28Z

Thank you very much for the review, @Mateasek. I've implemented all your suggestions.

Integrator1D now inherits from Function2D, which is a good idea. The integrand can be set via the integrand property in Python or directly accessed using the function attribute in the C-interface.

NumericallyIntegratedLineShapeModel is removed and the integrator attribute is added to the base LineShapeModel. The integrator is set to None by default except for the StarkBroadenedLine, which sets it to GaussianQuadrature.

Mateasek

Hi @vsnever , thanks for the changes. I think now the code is ready to be accepted.

jacklovell

I like this a lot: having more specialised integration routines to supplement trapezoidal integration is very useful indeed.

However, I don't think it makes sense to subclass the function framework classes for these integrators: just because Integrator1D.evaluate is a function of 2 arguments doesn't make the integrator a 2D function. For example, doing Integrator1D(function) + Arg1D() makes no sense.

Cherab has precedent for this: the core rate objects don't subclass from Function2D, even though they have evaluate(ne, te) methods. I think the same should apply here.

Otherwise, looks good. Couple of minor comments addressed inline.

jacklovell · 2022-10-13T15:48:38Z

+        if self.function is None:
+            raise AttributeError("Integrand is not set.")


Is there an expected use case for these integrators where the integrand is None? Ideally you would fail earlier than at the evaluation point: maybe require an integrand to be passed in when the object is instantiated and assert in the property setter that it cannot be None.

If Integrator1D is not a subclass of Function2D, then I see no point in making function an attribute of the class. I think it's better to pass it as an argument to evaluate():

cdef evaluate(self, Function1D function, double a, double b):

The wrapping can be done in the __call__ then:

def __call__(self, object function, double a, double b): return self.evaluate(autowrap_function1d(func), a, b)

Performing autowrap in every call could lead to a performance degradation if this is called many times: autowrap does a sequence of type checks and type conversions to validate the input. If these integrators are called many fewer times than the underlying function the performance hit is likely to be small though, so I don't see this as a major problem unless you can think of a use case where one might be for example sampling the output of the integral on a high resolution grid.

Ultimately I think either way is fine, but using the property you're effectively caching the result of the validations done using autowrap.

I would keep the function as a property. The integration is going to be called many times during observations

The performance of the Cython code will not suffer because evaluate is called directly in Cython, and the wrapping is done only in the Python __call__ method. But never mind, let's leave the integrand as a property.

The problem is that we can't make it mandatory to set the integrand on initialisation. This is so because an instance of Integrator1D is passed to the LineShapeModel constructor as an argument, but the integrand cannot be set prior LineShapeModel initialisation. If setting the integrand to None is a bad option, then it should have some default value, for example Constant1D(0).

…tegrand to None.

…lised line shape.

vsnever · 2022-10-17T19:27:07Z

Thanks for the review, @jacklovell. I fixed the issues you mentioned.

Also, I've added analytical line shape normalisation to StarkBroadenedLine using the expression:
$$\int_{-Ra}^{Ra}\frac{1}{|x|^{5/2}+a^{5/2}}dx=\frac{2R}{a^{3/2}}{_2}F_1(\frac{2}{5},1,\frac{7}{5},-R^{5/2}),$$
where a > 0, R > 0 and $_2F_1$ is a hypergeometric function.

vsnever added 4 commits August 30, 2022 20:35

Add Gaussian quadrature integrator and tests.

b202dbb

Fix an error in GaussianQuadrature. Increase default relative toleran…

a32a6a7

…ce value.

Add NumericallyIntegrableLineShapeModel class. Update StarkBroadenedL…

c4a604d

…ine to use numerical integration over the spectral bin.

Update changelog.

6308e0f

vsnever added the enhancement label Aug 31, 2022

Mateasek requested changes Sep 9, 2022

View reviewed changes

vsnever added 2 commits September 11, 2022 03:05

Make Integrator1D a Function2D instance.

526cc17

Remove NumericallyIntegratedLineShapeModel. Add integrator to LineSha…

86d908b

…peModel.

Mateasek approved these changes Sep 12, 2022

View reviewed changes

Merge branch 'development' into feature/quadrature_integration

803d529

jacklovell requested changes Oct 13, 2022

View reviewed changes

vsnever added 2 commits October 17, 2022 18:59

Remove Function2D inheritance for Integrator1D. Forbid setting the in…

6a915f9

…tegrand to None.

Add initialisation checks to StarkFunction and make it return a norma…

e62f0ad

…lised line shape.

vsnever and others added 2 commits October 18, 2022 23:14

Update CHANGELOG.md

308208c

Merge branch 'development' into feature/quadrature_integration

f0bc1aa

jacklovell merged commit 98223bc into cherab:development Dec 22, 2022

vsnever deleted the feature/quadrature_integration branch May 11, 2023 20:46

vsnever mentioned this pull request Dec 1, 2023

Add a numerical integration method to use with line shapes and plasma models #366

Closed

		def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, Integrator1D integrator=None,
		dict stark_model_coefficients=None):

		if self.function is None:
		raise AttributeError("Integrand is not set.")

Conversation

vsnever commented Aug 31, 2022

Uh oh!

Mateasek left a comment

Choose a reason for hiding this comment

Uh oh!

Mateasek Sep 9, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

vsnever commented Sep 11, 2022

Uh oh!

Mateasek left a comment

Choose a reason for hiding this comment

Uh oh!

jacklovell left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

vsnever commented Oct 17, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

Mateasek Sep 9, 2022 •

edited

Loading

vsnever commented Oct 17, 2022 •

edited

Loading