Finer steps sizes for Gauss-Legendre integration rules

[pkienzle]

The current set of Gauss-Legendre rules for adaptive integration uses npoints in [20, 76, 500, 5000].

This is okay for a single integration loop, but that step between 500 and 5000 makes the model run 100x slower when you have nested integrals.

However, more steps means more tuning of the cutoff values for qr_max. There is also a concern about the amount of memory used for GPU processes. Adding something around npoints=1000 will use another 16kB of constant memory.

The algorithm is too slow to generate points on the fly, so they need to be cached.

See also #717

[andyfaff]

@lru_cache(maxsize=128)
def gauss_legendre(n):
    """
    Calculate gaussian quadrature abscissae and weights
    Parameters
    ----------
    n : int
        Gaussian quadrature order.
    Returns
    -------
    (x, w) : tuple
        The abscissae and weights for Gauss Legendre integration.
    """
    return scipy.special.p_roots(n)

?

[andyfaff]

I think the underlying machinery to calculate them is also available in C.

[pkienzle]

Much faster than the numpy version!

In [2]: from scipy.special import p_roots

In [3]: from numpy.polynomial.legendre import leggauss

In [4]: leggauss(5)
Out[4]: 
(array([-0.90617985, -0.53846931,  0.        ,  0.53846931,  0.90617985]),
 array([0.23692689, 0.47862867, 0.56888889, 0.47862867, 0.23692689]))

In [5]: p_roots(5)
Out[5]: 
(array([-0.90617985, -0.53846931,  0.        ,  0.53846931,  0.90617985]),
 array([0.23692689, 0.47862867, 0.56888889, 0.47862867, 0.23692689]))

In [6]: %timeit leggauss(1000)
70.7 ms ± 346 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [7]: %timeit p_roots(1000)
20.1 ms ± 8.23 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [8]: %timeit p_roots(5000)
488 ms ± 174 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [9]: %timeit leggauss(5000)
14.2 s ± 1.37 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

On cuda we will need to put the cache in global memory since it will be too big for constant memory. Need to see what that does to performance.

[pkienzle]

I adjusted adaptive integration in #658 so that it limits the θ-φ mesh size to max 200000 points per q. Models are now fast enough and accurate enough for the ranges covered by SANS/USANS.

Still need to tune the qr cutoff values for each integration order, and consider adding orders between 76 and 500 and between 500 and 5000.

Use python explore/check_adaptive.py from #658 after tuning. It systematically tests every model for small and large shapes at different aspect ratios (rod, disk, cube).

[pkienzle]

To keep the overall runtime down I've limited the θ-φ grid to 100 000 points, with the outer loop limited to 500 points. We could instead check if the shape is a long rod (i.e., qr_max for the a-b cross section is small enough that a 20 point integration will suffice) and use 5000 points in this case. Handing both inner and outer qr max to the gauss_weights function would allow it to balance the integration order for the inner and outer loops.