Opened 10 years ago

# Upgrade numpy to 1.5.1 — at Version 16

Reported by: Owned by: jason tbd major sage-4.7 packages: standard evanandel, kcrisman N/A

New version of numpy is out: 1.5.1

### comment:1 Changed 10 years ago by fbissey

• Milestone changed from sage-4.6.2 to sage-4.7

Updating should be straightforward but some doctest will break because of some extra verbosity:

```sage -t -long -force_lib "devel/sage-main/sage/calculus/riemann.pyx"
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 102:
sage: m = Riemann_Map([lambda t: e^(I*t)], [lambda t: I*e^(I*t)], 0)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 563:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 613:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 676:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 760:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 822:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 152:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginaWarning: invalid value encountered in divide
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 102:
sage: m = Riemann_Map([lambda t: e^(I*t)], [lambda t: I*e^(I*t)], 0)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 563:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 613:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 676:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 760:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 822:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 152:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 196:
sage: isinstance(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)._repr_(), str)  # lo
ng time
Expected:
True
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
True
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 208:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 308:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 376:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 417:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 487:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/riemann.pyx", line 513:
sage: m = 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)  # long time (4 sec)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
```

and

```sage -t -long "devel/sage-main/sage/calculus/interpolators.pyx"
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/interpolators.pyx", line 52:
sage: m = Riemann_Map([lambda x: ps.value(x)], [lambda x: ps.derivative(x)],0)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
File "/usr/share/sage/devel/sage-main/sage/calculus/interpolators.pyx", line 183:
sage: m = Riemann_Map([lambda x: cs.value(x)], [lambda x: cs.derivative(x)], 0)
Expected nothing
Got:
doctest:1: ComplexWarning: Casting complex values to real discards the imaginary part
**********************************************************************
```

I think that's all the failures related to numpy-1.5.1

### comment:2 Changed 10 years ago by jason

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 fbissey

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 jason

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 jason

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 fbissey

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:8 Changed 10 years ago by jason

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 jason

• Status changed from new to needs_work

### comment:10 Changed 10 years ago by fbissey

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 jason

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?

### comment:12 Changed 10 years ago by jason

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)
**********************************************************************
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 jason

• Description modified (diff)

### comment:14 Changed 10 years ago by jason

• Status changed from needs_work to needs_info

### comment:15 follow-up: ↓ 16 Changed 10 years ago by 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

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 fbissey

• Description modified (diff)

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.

Note: See TracTickets for help on using tickets.