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

[David Goldsmith]

On Mon, Mar 8, 2010 at 6:52 AM, Bruce Southey <bsouthey at 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 at 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 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20100308/898f941c/attachment.html>