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 )
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)
Change History (16)
comment:1 Changed 9 years ago by
- 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
comment:2 Changed 9 years ago by
- Type changed from PLEASE CHANGE to defect
Changed 9 years ago by
comment:3 Changed 9 years ago by
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
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.
comment:5 follow-up: ↓ 6 Changed 9 years ago by
- 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
- 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
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
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.
comment:9 Changed 9 years ago by
- 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
- Cc fbissey added
comment:11 Changed 9 years ago by
- 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
- Milestone changed from sage-4.7 to sage-4.7.1
comment:13 Changed 9 years ago by
- Merged in set to sage-4.7.1.alpha0
- Resolution set to fixed
- Status changed from positive_review to closed
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:
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?