Add Gauss-Legendre integrator to documentation

[slel]

A method for Gauss-Legendre integration of fast-callable polynomials was added in #23140.

It lives in src/sage/numerical/gauss_legendre.pyx.

This ticket is to make it appear in the documentation.

Ideally add examples and do a bit of pep8 formatting.

Related:

• Ask Sage question 44545: Gauss-Legendre quadrature

CC: @slel

Component: numerical

Keywords: gauss_legendre

Author: Linden Disney-Hogg

Branch/Commit: e74a32b

Reviewer: Samuel Lelièvre, Nils Bruin

Issue created by migration from https://trac.sagemath.org/ticket/30698

[slel]

Description changed:

--- 
+++ 
@@ -1,7 +1,12 @@
 A method for Gauss-Legendre integration of fast-callable
 polynomials was added in #23140.

+It lives in `src/sage/numerical/gauss_legendre.pyx`.
+
 This ticket is to make it appear in the documentation.

-Ideally add examples and do a bit of pep8 reformatting.
+Ideally add examples and do a bit of pep8 formatting.

+Related:
+
+- [Ask Sage question 44545: Gauss-Legendre quadrature](https://ask.sagemath.org/question/44545)

[mkoeppe]

comment:3

Setting new milestone based on a cursory review of ticket status, priority, and last modification date.

[DisneyHogg]

Branch: u/gh-DisneyHogg/add_gauss_legendre_integrator_to_documentation

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

1f65c2d

Documentation edited and integrate_vector changed

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Commit: 1f65c2d

[nbruin]

comment:6

Looks good! Some comments based on an initial reading of the code:

• I don't think "REC" algorithm is universally accepted terminology. If you want to call it that, I think it needs a literature reference with it. This could be something like:

Compute the integration nodes and weights for the Gauss-Legendre quadrature scheme.

We use the recurrence relations for Legendre polynomials to compute their values. This is a version of the algorithm that in [Neurohr ...] is called the REC algorithm.

NOTE:

It may be worth testing if using the Arb algorithm for finding the nodes and weights in `arb/acb_calc/integrate_gl_auto_deg.c` has better performance.

• Concerning range in cython code: See Cython doc. This gets compiled to an actual c-for loop (it's important here that cython has type-information on what len returns (an actual c-integer-type) so that it can provide the loop with sufficient type hints so that the C compiler can actually compile it. The fact that xrange still works indicates that in these situations, cython considers range and xrange as synonyms and that this behaviour has been maintained (at least with the settings that sage uses for its cython library). Moving to range if you're touching the code anyway makes sense.

[slel]

comment:7

Remember to set the ticket to "needs review" when ready for review.

Some comments on the branch you pushed, in addition to those by @nbruin.

Split into two tickets?

I like your proposed change to integrate_vector allowing to specify a number of nodes instead of specifying epsilon, but for clarity I would separate that change out to a different ticket, and stick to documentation and reformatting here.

Now some comments on the changes in the existing branch.

Addition to the documentation

Well done adding a line to src/doc/en/reference/numerical/index.rst to make Gauss-Legendre integration appear in the documentation.

Reformatting

Good job on the reformatting overall.

Good call adding a double-colon changing NOTE: into NOTE::, and indenting its content.

To finish up, I think it needs to become .. NOTE::.

Some ways I would push the reformatting a bit further:

-        sage: L1 = nodes(24,53)
+        sage: L1 = nodes(24, 53)
-        sage: P = RR['x'](sage.functions.orthogonal_polys.legendre_P(24,x))
+        sage: P = RR['x'](sage.functions.orthogonal_polys.legendre_P(24, x))
         sage: Pdif = P.diff()
-        sage: L2 = [((r+1)/2,1/(1-r^2)/Pdif(r)^2) for r,_ in RR['x'](P).roots()]
+        sage: L2 = [((r + 1)/2, 1/(1 - r^2)/Pdif(r)^2)
+        ....:       for r, _ in RR['x'](P).roots()]
-        sage: all((a[0]-b[0]).abs() < 10^-15 and (a[1]-b[1]).abs() < 10^-9 for a,b in zip(L1,L2))
+        sage: all((a[0] - b[0]).abs() < 1e-15 and (a[1] - b[1]).abs() < 1e-9
+        ....:     for a, b in zip(L1, L2))
         True

         sage: prec = 200
         sage: K = RealField(prec)
-        sage: V = VectorSpace(K,2)
+        sage: V = VectorSpace(K, 2)
-        sage: epsilon = K(2^(-prec+4))
+        sage: epsilon = K(2^(-prec + 4))
-        sage: f = lambda t:V((1+t^2,1/(1+t^2)))
+        sage: f = lambda t: V((1 + t^2, 1/(1 + t^2)))
-        sage: I = integrate_vector(f, prec, epsilon)
+        sage: I = integrate_vector(f, prec, epsilon=epsilon)
-        sage: J = V((4/3, pi/4))
+        sage: J = V((4/3, pi/4))
-        sage: max(c.abs() for c in (I-J)) < epsilon
+        sage: max(c.abs() for c in I - J) < epsilon

         sage: prec = 200
         sage: Kreal = RealField(prec)
         sage: K = ComplexField(prec)
-        sage: V = VectorSpace(K,2)
+        sage: V = VectorSpace(K, 2)
-        sage: epsilon = Kreal(2^(-prec+4))
+        sage: epsilon = Kreal(2^(-prec + 4))
-        sage: f = lambda t: V((t,K(exp(2*pi*t*K.0))))
+        sage: f = lambda t: V((t, K(exp(2*pi*t*K.0))))
-        sage: I = integrate_vector(f, prec, epsilon)
+        sage: I = integrate_vector(f, prec, epsilon=epsilon)
-        sage: J = V((1/2,0))
+        sage: J = V((1/2, 0))
-        sage: max(c.abs() for c in (I-J)) < epsilon
+        sage: max(c.abs() for c in I - J) < epsilon

-        for i in range(1,len(nodelist)):
+        for i in range(1, len(nodelist)):

Changes to integrate_vector

Perhaps the new N=None keyword should come after the existing N=None keyword.

Otherwise, this line, which appears twice in the doctests:

        sage: I = integrate_vector(f, prec, epsilon)

is really executing

        sage: I = integrate_vector(f, prec, N=epsilon)

which is likely not intended, and this change might similarly break existing user code.

In any case, the examples should be explicit about which keyword is passed, using one of:

        sage: I = integrate_vector(f, prec, epsilon=epsilon)
        sage: I = integrate_vector(f, prec, N=N)

and possibly new examples should be added using

        sage: I = integrate_vector(f, prec, epsilon=epsilon, N=N)
        sage: I = integrate_vector(f, prec, N=N, epsilon=epsilon)

to illustrate that specifying N overrides specifying epsilon.

There is a typo "specificed" for "specified" in the new line describing N. Suggested change:

-     - ``N`` -- integer. Number of nodes to use. If specificed, the target error ``epsilon`` is ignored. 
+     - ``N`` -- integer. Number of nodes to use.
+       If specified, the target error ``epsilon`` is ignored. 

[nbruin]

comment:8

Replying to @slel:

Perhaps the new N=None keyword should come after the existing N=None keyword.

after the existing epsilon=None keyword, is probably what you meant. Good point! It actually suggests that override is not safe, since it's too easy to mix the two parameters up. So probably a ValueError (or something like that) should be raised if both N and epsilon have a non-None value. It would still not guard against accidental mix-ups due to positional parameters being confused with keyword parameters.

What this is actually indicative of is that the API of the routine is getting crowded and that there really are two routines here: one that integrates with a give number of nodes and one that uses heuristic error estimation to try and get a result within a given accuracy. The second routine may call the first for code reuse, but inlining the code would be acceptable for performance reasons.

So ... splitting is the way to go! It also leads to the fastest code because you select the branch by the name of the routine you call (and with cython that means the choice can be locked in compile-time!) rather than by the type of parameters that are passed in (which need runtime testing to determine in python and cython)

[mwageringel]

comment:9

Replying to @nbruin:

It would still not guard against accidental mix-ups due to positional parameters being confused with keyword parameters.

This could also be solved by using keyword-only arguments. In order to maintain backward compatibility, it could be implemented like this:

def integrate_vector(f, prec, epsilon_deprecated=None, *, epsilon=None, N=None):
    if epsilon_deprecated is not None:
        from sage.misc.superseded import deprecation
        deprecation(30698, "positional argument epsilon should be passed as keyword argument")
        if epsilon is None:
            epsilon = epsilon_deprecated
        else:
            raise ValueError("...")

[DisneyHogg]

comment:10

Replying to @slel:

Thanks for all the comments, will get cracking on them!

Split into two tickets?

I like your proposed change to integrate_vector allowing to specify a number of nodes instead of specifying epsilon, but for clarity I would separate that change out to a different ticket, and stick to documentation and reformatting here.

I have opened a new ticket #32014 for the changes to integrate_vector itself.

[nbruin]

comment:11

Replying to @mwageringel:

Replying to @nbruin:

It would still not guard against accidental mix-ups due to positional parameters being confused with keyword parameters.

This could also be solved by using keyword-only arguments. In order to maintain backward compatibility, it could be implemented like this:

def integrate_vector(f, prec, epsilon_deprecated=None, *, epsilon=None, N=None):
    if epsilon_deprecated is not None:
        from sage.misc.superseded import deprecation
        deprecation(30698, "positional argument epsilon should be passed as keyword argument")
        if epsilon is None:
            epsilon = epsilon_deprecated
        else:
            raise ValueError("...")

It doesn't remove the Code smell of boolean parameters, and these code-path-determining signatures have a similar smell. There is also the complexity: testing the function arguments is forcing us an extra runtime branch that is almost always resolvable at compile-time (because the branch is determined by essentially the call signature that is used, and writing code where the call signature is only determined at runtime is, while possible in python, difficult to do and, thankfully, rare). So distinguishing between the two uses with different callables is in this case definitely the right design. These routines are in fact used in rather tight loops. While in most applications any call overhead will be dwarfed by the required evaluations, it's near a point where you might not want to introduce extra overhead unnecessarily. I suspect the contribution will be very small, though, so the main argument here is really just clean code design.

(it's also an internal routine; at least it was up to this ticket, so deprecation wouldn't be necessary. I think in this case there's a cleaner solution that doesn't require API changes to existing routines at all! -- for those reading this and not familiar with the term "bikeshedding": this is a prime example of it, and it illustrates the limited productivity that ensues from it)

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

b55bec3

Format changes

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 1f65c2d to b55bec3

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

87886cc

Added reference for REC algorithm

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from b55bec3 to 87886cc

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 87886cc to f8c19a2

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

f8c19a2

Added TODO for nodes method

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from f8c19a2 to a56b3c1

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

a56b3c1

Undid changes to integrate_vector

[nbruin]

Reviewer: Nils Bruin

[nbruin]

comment:17

Looks good to me; checked the resulting HTML typesetting as well.

[slel]

comment:18

Missing continuation prompts:

         sage: L2 = [((r + 1)/2, 1/(1 - r^2)/Pdif(r)^2)
-                    for r, _ in RR['x'](P).roots()]
+        ....:       for r, _ in RR['x'](P).roots()]
         sage: all((a[0] - b[0]).abs() < 1e-15 and (a[1] - b[1]).abs() < 1e-9
-                  for a, b in zip(L1, L2))
+        ....:     for a, b in zip(L1, L2))

If you fix that, set back to positive review on my behalf.

[slel]

Changed reviewer from Nils Bruin to Samuel Lelièvre, Nils Bruin

[vbraun]

comment:19

Author name is missing

[slel]

Author: Linden Disney-Hogg

[nbruin]

comment:21

Thanks, Samuel.

The continuation prompts indeed should be fixed. I'm setting the ticket back to "needs work" because it does, and otherwise it might get missed.

Once Linden has added continuation prompts to those two lines, he can set it back to "positive"; or alert me to do that.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

e74a32b

Added continuation prompts

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from a56b3c1 to e74a32b

[vbraun]

Changed branch from u/gh-DisneyHogg/add_gauss_legendre_integrator_to_documentation to e74a32b