Replace tabulated Gaussian quadrature with rule

[rjmoerland]

Rationale

Gaussian quadrature is an efficient scheme for numerical integration. The current GaussianQuadrature distribution in Optiland, for integration over the unit disk, uses tabulated numbers for orders 1 to 6. Higher order integration is therefore limited to other distribution types, that may be suboptimal, such as uniform or random sampling.

This PR implements a new GaussianQuadrature class, a subclass of BaseDistribution, where the number of rings and spokes can be adjusted at will. The correct locations are calculated on the fly, together with the integration weights.

Implementation

The implementation is based on the paper by William H. Peirce, "Numerical Integration Over the Planar Annulus,", Journal of the Society for Industrial and Applied Mathematics, Vol. 5, No. 2 (Jun., 1957), pp. 66-73.

I believe it is a drop-in replacement for many cases, but I still have open questions about the weights:

1. The typical normalization of the weights is such that the sum is equal to 1.0 or π. However, the current Optiland implementation uses a value of 1.0 if is_symmetric is False and 3.0 if is_symmetric is True. Why the difference?
2. The description of is_symmetric and the implementation seem to differ - according to the text, is_symmetric is to indicate symmetry in y, but activating it seems to indicate that the distribution is circularly symmetric. Not activating it seems to imply mirror symmetry in x. Is this correct, and should it be maintained?
3. It is not clear to me what exactly is the purpose of the symmetric/asymmetric distribution toggle in ray.py, when Hx == Hy == 0. I tried to maintain the functionality of OPD_difference, but I will need a little help because I don't understand the weights 😶

As a consequence of my lack of understanding of these details, there are a few pytests failing.

[HarrisonKramer]

Thanks for adding this, @rjmoerland!

I recognize the need to move away from a direct lookup here.

Using your numbering:

1. I believe this was a simplified attempt to apply a scale factor based on the number of points used. It's currently hardcoded to use 1 spoke for the symmetric case and 3 for the asymmetric case. This is the source of that scale factor.
2. You are right. It's currently written incorrectly in the docs.
3. This is an attempt to enforce that a rotationally symmetric field (Hx=Hy=0) only uses a single spoke for computational efficiency. I am okay with abandoning the manual weights approach used and moving towards your more mathematically rigorous approach. Perhaps we can default to one spoke for the symmetric case and 3 for the asymmetric case, but it's also fine to let users configure that themselves.

What do you think? I'd say you can update the OPD difference function to make it mathematically consistent, then update the corresponding unit tests.

Thanks,
Kramer

[rjmoerland]

Hi @HarrisonKramer , thank you for the explanations. Keeping the numbering, but swapping them around a little:

• 3: I think the implied boundary conditions for the current implementation of OPD_difference are that:
  a. Hx == Hy == 0: The lens system is rotationally symmetric
  b. else: The lens system has mirror symmetry along the x axis. Implicitly, that means the lens system is rotationally symmetric and Hx == 0.
  I think it's dangerous to always assume these to be true, as a cylindrical lens or an off-axis element is sufficient to break the symmetry. I would strongly recommend to return distributions that "always" work, and leave it to a the user to make an informed choice to use symmetry.

• 2: Pending on what we decide on 3, I'll update the documentation accordingly.

• 1: I do see the use for having distributions that use some kind of symmetry. I could look at some switch to reduce the number of points in x, y, or both. I recommend returning weights that are normalized such that 'integrating' the constant function 1.0 simply returns area/pi (that is, 1), regardless of the used symmetry. This means users do not need to branch or adjust the outcome of a numerical integration when they use symmetry.

I'll gladly update the OPD_difference function, but is there a standard model that is typically used in geometrical optics, with known analytical results? Using the Hubble model seems to test more for consistency than accuracy (but then again, I don't know much about the Hubble's optical system).

[HarrisonKramer]

Hi @rjmoerland,

Same numbering:

• 3: Very good point. It's indeed risky to assume symmetry. We could check symmetry of the system manually (e.g., via the is_symmetric attribute of geometries), but this is likely overkill. I agree with your proposal. Users can manually change this, if they wish.

• 1: It would be great to exploit symmetry where possible. More rays for a single raytrace of course results in much more computation when run over many optimization iterations. As you mentioned, a fallback could be to simply let users choose this. That's likely the best option. And I like your recommendation for integration.

AFAIK, there's no standard model used here. A Gaussian quadrature with N radial points should perfectly integrate a polynomial up to order 2N-1, based on the Forbes paper below. To check the accuracy of the integration approach, we can either compare to a mathematicall constructed pupil (e.g., some high-dimensional Zernike representation), or compare to a very high density raytraced grid, treating it as "exact". I think Forbes used the latter. If we wanted to be thorough, we could check for low-NA, high-NA, low-FOV, high-FOV, and refractive/reflective systems.

G. W. Forbes, "Optical system assessment for design: numerical ray tracing in the Gaussian pupil," J. Opt. Soc. Am. A, Vol. 5, No. 11, pp. 1943-1956 (November 1988).

[rjmoerland]

Hi @rjmoerland,

Same numbering:

3: Very good point. It's indeed risky to assume symmetry. We could check symmetry of the system manually (e.g., via the is_symmetric attribute of geometries), but this is likely overkill. I agree with your proposal. Users can manually change this, if they wish.

Cool, agreed.

1: It would be great to exploit symmetry where possible. More rays for a single raytrace of course results in much more computation when run over many optimization iterations. As you mentioned, a fallback could be to simply let users choose this. That's likely the best option. And I like your recommendation for integration.

Agreed. I'll check what I can do. I'm thinking along the lines of extra switches like symmetry_x and symmetry_y that could be True or False, or, alternatively, a string like x_symmetry, y_symmetry or xy_symmetry, which would allow other options to be added in the future. But before implementing this, I'll check if such a symmetry condition can be applied robustly.

AFAIK, there's no standard model used here. A Gaussian quadrature with N radial points should perfectly integrate a polynomial up to order 2N-1, based on the Forbes paper below. To check the accuracy of the integration approach, we can either compare to a mathematicall constructed pupil (e.g., some high-dimensional Zernike representation), or compare to a very high density raytraced grid, treating it as "exact". I think Forbes used the latter. If we wanted to be thorough, we could check for low-NA, high-NA, low-FOV, high-FOV, and refractive/reflective systems.

G. W. Forbes, "Optical system assessment for design: numerical ray tracing in the Gaussian pupil," J. Opt. Soc. Am. A, Vol. 5, No. 11, pp. 1943-1956 (November 1988).

My apologies for going completely off track here, but: I think the statement in the cited paper is possibly a mistake, or somewhat misleading. Indeed, Gaussian quadrature with N points can integrate a polynomial of degree 2N-1, but that is for the 1D case (meaning something like $$y = f(x)$$) as far as I can tell. Since the tabulated radial evaluation positions are in agreement with William H. Peirce, "Numerical Integration Over the Planar Annulus,", Journal of the Society for Industrial and Applied Mathematics, Vol. 5, No. 2 (Jun., 1957), pp. 66-73, and the weights are the same except for a constant factor, I think it's safe to assume that the underlying conditions are the same as well. According to Peirce, the recipe given in his paper is good for integrating 2D polynomials of order N in both $x$ and $y$. Implicitly, this means that the function spanning the unit disk has to be smooth and continuous. Now, if we would look at the simple radial polynomial $$f(r) = r$$, then according to Forbes, N = 1 would suffice. On the other hand, if we treat it as a 2D function, it actually has a discontinuous first order derivative at $r = 0$ (or perhaps even undefined), and this discontinuity is difficult to approximate with polynomials. In fact, as a 2D function in $x$ and $y$, it would need to be written as $$f(x, y) = \sqrt{x^2 + y^2}$$, which highlights how it is not a polynomial in $x$ and $y$.

I think I can provide and example of the above by doing the numerical integration. The analytical result for the integral of $f(r) = r$ over the unit disk is $\int_0^{2 \pi} \int_0^1 f(r) r dr d\phi = 2 \pi \int_0^1 r\ r d r = 2\pi/3$. For N = 1 and the rule from Forbes/Peirce, the evaluation point is $r = \sqrt{1/2}$. The weight is (in my implementation) 1.0. Therefore, the approximation of the integral above is $\sqrt{1/2} \pi$, which is still quite far from $2 \pi / 3$. Increasing the number of points to 15 improves the result to $0.6666943765953873 \pi$.

On the other hand, if you know that the function is circularly symmetric, then instead of approximating the 2D integral, it is more efficient to evaluate $2 \pi \int_0^1 f(r) r dr = 2\pi \int_0^1 r^2 dr$. In that case, N = 2 is the minimum required order, and shifting the domain from -1... +1 to 0 ... 1 makes that the evaluation points are r1 = 0.21132487, r2 = 0.78867513. Evaluating $f(r) r $, and summing these after weighting gives $0.6666666666666666 \pi$, which is much more accurate than the 2D evaluation.

Anyway, this is as far as I'll bore you, but I think it's useful to establish what the 2D Gaussian quadrature can do and what it cannot do.

Then, more on topic: I have implemented a few pytest to check for the accuracy of the numerical integration, but what I meant to ask is: are there any models for the Optical Path Difference that are typically used in geometrical optics? Something like the analytical result for the Optical Path Difference when a collimated beam is focused by a plano-convex lens for example? Anything that is non-trivial, but has an analytical solution would be perfect for this.

Thanks!

[rjmoerland]

Just in case you might have missed it on Discord:

I've been looking at that OPD function for this PR, but I cannot figure out what the correct value should be for the failing pytest (OPD of Hubble telescope). I also couldn't really match OPD values to a toy model I made for a plano-convex lens, so I must be misunderstanding something here. Is this something we can discuss?