Opened 10 years ago
Last modified 10 years ago
#10792 closed task
Upgrade numpy to 1.5.1 — at Version 16
Reported by: | jason | Owned by: | tbd |
---|---|---|---|
Priority: | major | Milestone: | sage-4.7 |
Component: | packages: standard | Keywords: | |
Cc: | evanandel, kcrisman | Merged in: | |
Authors: | Reviewers: | ||
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
New version of numpy is out: 1.5.1
SPKG: http://spkg-upload.googlecode.com/files/numpy-1.5.1.spkg
Change History (17)
comment:1 Changed 10 years ago by
- Milestone changed from sage-4.6.2 to sage-4.7
comment:2 Changed 10 years ago by
It sounds like there may be just a small change we should make to the Riemann map function. Do you know line it is complaining about?
comment:3 Changed 10 years ago by
I am not sure. I don't understand the code very well in Riemann map. The only thing suspect to me is this
pt1 = np.complex(pt)
on line 499 in riemann.pyx. pt1 type has never been defined so I don't know if cython automatically give it the type of the right hand side or if it defaults to something else. In which case it would have to be cast.
comment:4 Changed 10 years ago by
Okay. I'll try to do some printf debugging to see if I can find it "soon".
comment:5 follow-up: ↓ 6 Changed 10 years ago by
Do you have a 1.5.1 spkg that you could post? It sounds like you already made one.
comment:6 in reply to: ↑ 5 Changed 10 years ago by
Replying to jason:
Do you have a 1.5.1 spkg that you could post? It sounds like you already made one.
I get that a lot. I will prepare one soon. But I don't need a new spkg to test things in sage-on-gentoo because I can use the system version, which is how I am aware of this issue and a few others regarding upgrading components.
comment:7 Changed 10 years ago by
You can find a numpy-1.5.1 spkg here http://spkg-upload.googlecode.com/files/numpy-1.5.1.spkg
comment:8 Changed 10 years ago by
I don't understand the riemann code very well either, but it seems like the warning is coming from the line where temp (a real array?) is being assigned a value from self.cps (a complex array)---see the code that is deleted in the trac-10792-riemann.patch file. I simplified the code a lot to avoid having a temporary array, but this changed lots of numerical results. So one question is: are the results it returns now better or worse than before?
I also seem to get a huge number of Warning: invalid value encountered in divide
warnings when doing Riemann_Map([lambda t: e^(I*t) - 0.5*e^(-I*t)], [lambda t: I*e^(I*t) + 0.5*I*e^(-I*t)], 0)
. These come from NaN values popping up in the array. Does anyone else see those?
comment:9 Changed 10 years ago by
- Status changed from new to needs_work
comment:10 Changed 10 years ago by
The warnings have been there since cython-0.13 has landed. The code itself is numerically sensitive to various details including lapack implementation. You may have seen something about it on sage-release.
I will take a few slow days because the town where I live (Christchurch, New Zealand) has been hit by a 6.3 earthquake. Where I live is ok, where I work isn't.
comment:11 Changed 10 years ago by
- Cc evanandel added
CCing Ethan, since he was the original author on the riemann.pyx file.
Ethan: I've been looking at the code that seems to cause two problems. First is a ComplexWarnings? that now appears in numpy 1.5.1, which I believe comes from a temporary array that was assumed to be real, but maybe should be complex. I've rewritten the affected lines to not use a temporary array.
The other problem seems like it is a problem even before the upgrade---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)
I noticed several things that could be pulled out of the inner loop, and so rewrote this line in the attached patch (which I'll update shortly). 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 this entry and not calculate it?
Changed 10 years ago by
comment:12 Changed 10 years ago by
Ethan,
I've split the issue into two problems: the ComplexWarnings? problem that I hopefully fixed with the patch on this ticket and the invalid division warning, which I've made #10821. Can you check this code and the numerical differences in output? I now get these failures:
% sage -t -long riemann.pyx | grep -v "Warning: invalid value encountered in divide" E:128 Detected SAGE64 flag Building Sage on OS X in 64-bit mode sage -t -long "trees/sage-4.6.1/devel/sage-main/sage/calculus/riemann.pyx" ********************************************************************** File "/Users/grout/sage-trees/sage-4.6.1/devel/sage-main/sage/calculus/riemann.pyx", line 123: sage: print "error =", m.inverse_riemann_map(m.riemann_map(x)) - x # long time Expected: error = (-0.0001211863...+0.001606139...j) Got: error = (-0.000291742263558+0.00160874617675j) ********************************************************************** File "/Users/grout/sage-trees/sage-4.6.1/devel/sage-main/sage/calculus/riemann.pyx", line 484: sage: m.riemann_map(0.25 + sqrt(-0.5)) # long time Expected: (0.137514137885...+0.876696023004...j) Got: (0.13751414123737432+0.87669602715226935j) ********************************************************************** File "/Users/grout/sage-trees/sage-4.6.1/devel/sage-main/sage/calculus/riemann.pyx", line 486: sage: m.riemann_map(1.3*I) # long time Expected: (-1.561029396...+0.989694535737...j) Got: (-1.560925774791804e-05+0.98969453688560927j) ********************************************************************** File "/Users/grout/sage-trees/sage-4.6.1/devel/sage-main/sage/calculus/riemann.pyx", line 489: sage: m.riemann_map(0.4) # long time Expected: (0.733242677182...+3.50767714...j) Got: (0.73324267718110436+3.2076771460460016e-06j) ********************************************************************** File "/Users/grout/sage-trees/sage-4.6.1/devel/sage-main/sage/calculus/riemann.pyx", line 492: sage: m.riemann_map(np.complex(-3, 0.0001)) # long time Expected: (1.405757135...+1.05106...j) Got: (1.4057571374988058e-05+8.0616369487397992e-10j) ********************************************************************** 2 items had failures: 1 of 15 in __main__.example_1 4 of 9 in __main__.example_8 ***Test Failed*** 5 failures. For whitespace errors, see the file /Users/grout/.sage//tmp/.doctest_riemann.py [105.0 s] ---------------------------------------------------------------------- The following tests failed: sage -t -long "trees/sage-4.6.1/devel/sage-main/sage/calculus/riemann.pyx" Total time for all tests: 105.1 seconds
It looks like the errors are all slight numerical errors. Can someone confirm that?
comment:13 Changed 10 years ago by
- Description modified (diff)
comment:14 Changed 10 years ago by
- Status changed from needs_work to needs_info
comment:15 follow-up: ↓ 16 Changed 10 years ago by
Sorry for the slow response.
It looks like everything is slight numerical problems, which is expected given (as mentioned above) the numerical sensitivity. I'll happily test, but I'm not sure how to upgrade Sage's numpy. If someone can tell me how to get that running, I'll be able to confirm this patch and #10821.
Side Note 1: I see that my code is apparently quite incomprehensible to people. The underlying mathematics is quite a bear, so it'll never be nice and accessible, but is there anything I can do to improve the documentation?
Side Note 2: What's the status of fast_callable these days? Patch #8867 is still out there, and if fast_callable is ready, I think we might be able to merge that as well (otherwise, I can try amending it to work with fast_callable as is)
comment:16 in reply to: ↑ 15 Changed 10 years ago by
- Description modified (diff)
Replying to evanandel:
Sorry for the slow response.
It looks like everything is slight numerical problems, which is expected given (as mentioned above) the numerical sensitivity. I'll happily test, but I'm not sure how to upgrade Sage's numpy. If someone can tell me how to get that running, I'll be able to confirm this patch and #10821.
You can go through the following steps:
1) drop the numpy spkg in the spkg/standard folder
2) sage -f numpy-1.5.1
3) apply the patches (from here and #10821 ) the easiest way is to follow the procedure from the developer's manual http://www.sagemath.org/doc/developer/walk_through.html#reviewing-a-patch
4) sage -b
5) test
Side Note 1: I see that my code is apparently quite incomprehensible to people. The underlying mathematics is quite a bear, so it'll never be nice and accessible, but is there anything I can do to improve the documentation?
could you refer to a review of it in a journal or an online resource?
Side Note 2: What's the status of fast_callable these days? Patch #8867 is still out there, and if fast_callable is ready, I think we might be able to merge that as well (otherwise, I can try amending it to work with fast_callable as is)
No idea.
Updating should be straightforward but some doctest will break because of some extra verbosity:
and
I think that's all the failures related to numpy-1.5.1