riemann.pyx gives lots of invalid value in divide warnings

[jason]

I keep getting Warning: invalid value encountered in divide, and the error happens in these lines in the _generate_theta_array function in calculus/riemann.pyx:

        K = np.array(
            [-TWOPI / N * sadp * sadp[t] * 1 / (TWOPI * I) *
              ((dp / adp) / (cp - cp[t]) -
               ((dp[t] / adp[t]) / (cp - cp[t])).conjugate())
              for t in np.arange(NB)], dtype=np.complex128)

[jason]

Ethan, I noticed several things that could be pulled out of the inner loop, and so rewrote this line in the attached patch. However, I don't know this code or the algorithms very well---can you look at the patch?

It appears that the warnings happen when there is a division by zero, as in this test case:

sage: import numpy as np
sage: a=np.array([1,2],dtype=np.complex128)
sage: b=np.array([0,1],dtype=np.complex128)
sage: a/b
Warning: invalid value encountered in divide
array([ nan nanj,   2. +0.j])
sage: a[0]/b[0]
(nan+nan*j)

So should we turn off the error checking for this error when doing this operation since at least one entry in cp-cp[t] will be zero? Or should we check for zero denominators and not calculate the fraction with those?

[jason]

Unfortunately, it appears that I really changed something in my patch, or introduced a significant amount of numerical error, as now the doctests give results that are way different than documented.

[fbissey]

I am not sure what to do with arrays quoted with indices and without, as in:

cp -cp[t]

but I am fairly sure you attached .conjugate() to the wrong expression. I am attaching a corrected patch. I cannot try it right now.

[fbissey]

my own version of the patch

[evanandel]

This code generates a square array of dimension NB. The 't'th row is generated by (among other things) computing the difference between every element of cp and cp[t]. That's the meaning of

cp - cp[t]

fbissey: It looks like your patch is correct, I'll try it shortly.

Jason: If it interests you, this section of code is preparing a square matrix for nystrom numerical integration. The cp array holds the collocation points, that is the points around the boundary of the figure where the integral is being evaluated. dp contains the derivatives at those points.

computing the square array will of course end up with illegal divisions for the 't'th element of the 't'th row. However, the next line of code:

for i in xrange(NB):
    K[i, i] = 1

overwrites the bad elements with 1's. Thus the divisions by zero don't affect the algorithm. (If we have any other divisions by 0, that means that the user has tried to evaluate a self intersecting figure (and gotten very unlucky) for which the algorithm is meaningless anyway.)

We can deal with the zero divides whatever way fits sage style best. If there's an efficient way to not perform those divisions that's fine. Otherwise I have no problems with simply ignoring the warnings.

• Ethan

[fbissey]

Replying to evanandel:

We can deal with the zero divides whatever way fits sage style best. If there's an efficient way to not perform those divisions that's fine. Otherwise I have no problems with simply ignoring the warnings.

I will think about that a little bit but given the context we can indeed ignore these warnings. I think they show up in stderr and don't affect the results of the tests. They are obvious when the tests fail but I am not sure they are otherwise.

[fbissey]

Do I have the syntax right in thinking that K off-diagonal elements are defined by

K[i,j]= -TWOPI / N * sadp[i] * sadp[j] * 1 / (TWOPI * I) *
              ((dp[i] / adp[i]) / (cp[i] - cp[j]) -
               ((dp[j] / adp[j]) / (cp[i] - cp[j])).conjugate())

or do I have i and j mixed up?

[kcrisman]

Seeing this made me look through riemann.pyx, and it turns out that file could use some general TLC. Esp. with regard to a few doc things and some fairly massive redundancy with complex plots. See #10945.

[evanandel]

slight fix

[evanandel]

Alright, I wrapped the offending code in calls that tell numpy to ignore those invalid divide warnings then restore the original settings. It seems to fix things well.

[kcrisman]

Okay, the code seems right, (and slightly refactored but fine), conforms with Numpy specs, passes the tests, and Ethan says it's correct for his algorithm.