Opened 9 years ago

Closed 9 years ago

#10821 closed defect (fixed)

riemann.pyx gives lots of invalid value in divide warnings

Reported by: jason Owned by: burcin
Priority: major Milestone: sage-4.7.1
Component: calculus Keywords:
Cc: evanandel, fbissey Merged in: sage-4.7.1.alpha0
Authors: François Bissey, Jason Grout Reviewers: Ethan Van Andel, Karl-Dieter Crisman
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description (last modified by evanandel)

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)

Attachments (3)

trac-10821-riemann-invalid-divide.patch (1.1 KB) - added by jason 9 years ago.
trac-10821-riemann-invalid-divide.2.patch (1.1 KB) - added by fbissey 9 years ago.
my own version of the patch
trac-10821-riemann-invalid-divide.3.patch (1.3 KB) - added by evanandel 9 years ago.
slight fix

Download all attachments as: .zip

Change History (16)

comment:1 Changed 9 years ago by jason

  • Cc evanandel added
  • Component changed from PLEASE CHANGE to calculus
  • Description modified (diff)
  • Milestone set to sage-4.6.2
  • Owner changed from tbd to burcin

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?

comment:2 Changed 9 years ago by jason

  • Type changed from PLEASE CHANGE to defect

Changed 9 years ago by jason

comment:3 Changed 9 years ago by 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.

comment:4 Changed 9 years ago by 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.

Changed 9 years ago by fbissey

my own version of the patch

comment:5 follow-up: Changed 9 years ago by evanandel

  • Description modified (diff)

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

comment:6 in reply to: ↑ 5 Changed 9 years ago by fbissey

  • Milestone changed from sage-4.6.2 to sage-4.7
  • Status changed from new to needs_work

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.

comment:7 Changed 9 years ago by 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?

comment:8 Changed 9 years ago by 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.

Changed 9 years ago by evanandel

slight fix

comment:9 Changed 9 years ago by evanandel

  • Status changed from needs_work to needs_review

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.

comment:10 Changed 9 years ago by fbissey

  • Cc fbissey added

comment:11 Changed 9 years ago by kcrisman

  • Authors set to François Bissey, Jason Grout
  • Reviewers set to Ethan Van Andel, Karl-Dieter Crisman
  • Status changed from needs_review to positive_review

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.

comment:12 Changed 9 years ago by jdemeyer

  • Milestone changed from sage-4.7 to sage-4.7.1

comment:13 Changed 9 years ago by jdemeyer

  • Merged in set to sage-4.7.1.alpha0
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.