[Numpy-discussion] Is this a bug in numpy.ma.reduce?

Bruce Southey bsouthey@gmail....
Mon Mar 8 13:35:42 CST 2010


On 03/08/2010 12:17 PM, David Goldsmith wrote:
> On Mon, Mar 8, 2010 at 6:52 AM, Bruce Southey <bsouthey@gmail.com 
> <mailto:bsouthey@gmail.com>> wrote:
>
>     On 03/08/2010 01:30 AM, David Goldsmith wrote:
>>     On Sun, Mar 7, 2010 at 4:41 AM, Friedrich Romstedt
>>     <friedrichromstedt@gmail.com
>>     <mailto:friedrichromstedt@gmail.com>> wrote:
>>
>>         I would like to stress the fact that imo this is maybe not
>>         ticket and not a bug.
>>
>>         The issue arises when calling a.max() or similar of empty
>>         arrays a, i.e., with:
>>
>>         >>> 0 in a.shape
>>         True
>>
>>         Opposed to the .prod() of an empty array, such a .max() or .min()
>>         cannot be defined, because the set is empty.  So it's fully
>>         correct to
>>         let such calls fail.  Just the failure is a bit deep in
>>         numpy, and
>>         only the traceback gives some hint what went wrong.
>>
>>         I posted something similar also on the matplotlib-users list,
>>         sorry
>>         for cross-posting thus.
>>
>>
>>     Any suggestions, then, how to go about figuring out what's
>>     happening in my code that's causing this "feature" to manifest
>>     itself?
>>
>>     DG
>     Perhaps providing the code with specific versions of Python, numpy
>     etc. would help.
>
>     I would guess that aquarius_test.py has not correctly setup the
>     necessary inputs (or has invalid inputs) required by matplotlib
>     (which I have no knowledge about). Really you have to find if the
>     _A in cmp.py used by 'self.norm.autoscale_None(self._A)' is valid.
>     You may be missing a valid initialization step because the
>     TypeError exception in autoscale_None ('You must first set_array
>     for mappable') implies something need to be done first.
>
>     Bruce
>
>
> Python 2.5.4, Numpy 1.4.0, Matplotlib 0.99.0, Windows 32bit Vista Home 
> Premium SP2
>
> # Code copyright 2010 by David Goldsmith
> # Comments and unnecessaries edited for brevity
> import numpy as N
> import matplotlib as MPL
> from matplotlib import pylab
> from matplotlib.backends.backend_agg import FigureCanvasAgg as 
> FigureCanvas
> from matplotlib.figure import Figure
> import matplotlib.cm <http://matplotlib.cm> as cm
>
> J = complex(0,1); tol = 1e-6; maxiter = 20;
> roots = (-2 + J/2, -1 + J,  -0.5 + J/2, 0.5 + J,  1 + J/2, 2 + J, 2.5 
> + J/2,
>          -2 - J,   -1 - J/2, -0.5 - J,  0.5 - J/2, 1 - J,  2 - J/2, 
> 2.5 - J)
>
> def ffp(z):
>     w, wp = (0J, 0J)
>     for root in roots:
>         z0 = z - root
>         w += N.sin(1/z0)
>         wp -= N.cos(1/z0)/(z0*z0)
>     return (w, wp)
>
> def iter(z):
>     w, wp = ffp(z)
>     return z - w/wp
>
> def find_root(z0):#, k, j):
>     count = 0
>     z1 = iter(z0)
>     if N.isnan(z1):
>         return N.complex64(N.inf)
>     while (N.abs(z1 - z0) > tol) and \
>            (count < maxiter):
>         count += 1
>         z0 = z1
>         z1 = iter(z0)
>     if N.abs(z1 - z0) > tol:
>         result = 0
>     else:
>         result = z1
>     return N.complex64(result)
>
> w, h, DPI = (3.2, 2.0, 100)
> fig = Figure(figsize=(w, h),
>              dpi=DPI,
>              frameon=False)
> ax = fig.add_subplot(1,1,1)
> canvas = FigureCanvas(fig)
> nx, xmin, xmax = (int(w*DPI), -0.5, 0.5)
> ny, ymin, ymax = (int(h*DPI),  0.6, 1.2)
>
> X, xincr = N.linspace(xmin,xmax,nx,retstep=True)
> Y, yincr = N.linspace(ymin,ymax,ny,retstep=True)
> W = N.zeros((ny,nx), dtype=N.complex64)
>
> for j in N.arange(nx):
>     if not (j%100): # report progress
>         print j
>     for k in N.arange(ny):
>         x, y = (X[j], Y[k])
>         z0 = x + J*y
>         W[k,j] = find_root(z0)#,k,j)
>
> print N.argwhere(N.logical_not(N.isfinite(W.real)))
> print N.argwhere(N.logical_not(N.isfinite(W.imag)))
> W = W.T
> argW = N.angle(W)
> print N.argwhere(N.logical_not(N.isfinite(argW)))
> cms = ("Blues",)# "Blues_r", "cool", "cool_r",
>
> def all_ticks_off(ax):
>     ax.xaxis.set_major_locator(pylab.NullLocator())
>     ax.yaxis.set_major_locator(pylab.NullLocator())
>
> for cmap_name in cms:
>     all_ticks_off(ax)
>     ax.hold(True)
>     for i in range(4):
>         for j in range(4):
>             part2plot = argW[j*ny/4:(j+1)*ny/4, i*nx/4:(i+1)*nx/4]
>             if N.any(N.logical_not(N.isfinite(part2plot))):
>                 print i, j,
>                 print N.argwhere(N.logical_not(N.isfinite(part2plot)))
>             extent = (i*nx/4, (i+1)*nx/4, (j+1)*ny/4, j*ny/4)
>             ax.imshow(part2plot, cmap_name, extent = extent)
>             ax.set_xlim(0, nx)
>             ax.set_ylim(0, ny)
>     canvas.print_figure('../../Data-Figures/Zodiac/Aquarius/'+ cmap_name +
>                         'Aquarius_test.png', dpi=DPI)
> # End Aquarius_test.png
>
> DG
>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>    
Hmm,
Appears that you have mixed your indices when creating part2plot. If you 
this line instead it works:
part2plot = argW[j*nx/4:(j+1)*nx/4, i*ny/4:(i+1)*ny/4]


I found that by looking the shape of the part2plot array that is 
component of the argW array.

The shape of argW is (320, 200). So in your loops to find part2plot you 
eventually exceed 200 and eventually the index to the second axis is 
greater than 200 causing everything to crash:
When i=0 or 1 then the shape of part2plot is (50, 80)
when i=2 then the shape of part2plot is (50, 40)
when i=3 then the shape of part2plot is (50,  0) # crash



Bruce

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/numpy-discussion/attachments/20100308/cd5e025c/attachment.html 


More information about the NumPy-Discussion mailing list