1 Dec 2010 04:22

### Re: A faster median (Wirth's method)

I think last time I looked into how to apply a function along an axis I thought that the PyArray_IterAllButAxis would not work for that task (http://docs.scipy.org/doc/numpy/reference/c-api.array.html#PyArray_IterAllButAxis), but I think perhaps I misunderstood it. I'm looking into how to use it.

On Tue, Nov 30, 2010 at 12:06 PM, Keith Goodman wrote:
On Tue, Nov 30, 2010 at 11:58 AM, Matthew Brett <matthew.brett <at> gmail.com> wrote:
> Hi,
>
> On Tue, Nov 30, 2010 at 11:35 AM, Keith Goodman <kwgoodman <at> gmail.com> wrote:
>> On Tue, Nov 30, 2010 at 11:25 AM, John Salvatier
>> <jsalvati <at> u.washington.edu> wrote:
>>> I am very interested in this result. I have wanted to know how to do an
>>
>> My first thought was to write the reducing function like this
>>
>> cdef np.float64_t namean(np.ndarray[np.float64_t, ndim=1] a):
>>
>> but cython doesn't allow np.ndarray in a cdef.
>
> Sorry for the ill-considered hasty reply, but do you mean that this:
>
> import numpy as np
> cimport numpy as cnp
>
> cdef cnp.float64_t namean(cnp.ndarray[cnp.float64_t, ndim=1] a):
>    return np.nanmean(a)  # just a placeholder
>
> is not allowed?  It works for me.  Is it a cython version thing?
> (I've got 0.13),

Oh, that's nice! I'm using 0.11.2. OK, time to upgrade.
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion <at> scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

```_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion <at> scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
```
1 Dec 2010 06:46

### Re: A faster median (Wirth's method)

```> > import numpy as np
> > cimport numpy as cnp
>
> > cdef cnp.float64_t namean(cnp.ndarray[cnp.float64_t, ndim=1] a):
> >    return np.nanmean(a)  # just a placeholder
>
> > is not allowed?  It works for me.  Is it a cython version thing?
> > (I've got 0.13),
>
> Oh, that's nice! I'm using 0.11.2. OK, time to upgrade.

Oh wow, does that mean that http://trac.cython.org/cython_trac/ticket/177
is fixed? I couldn't find anything in the release notes about that,
but it would be great news. Does the cdef function acquire and hold
the buffer?

Felix
```
1 Dec 2010 11:23

### Re: numpy speed question

```Le 26/11/2010 17:48, Bruce Sherwood a écrit :
> Although this was mentioned earlier, it's worth emphasizing that if
> you need to use functions such as cosine with scalar arguments, you
> should use math.cos(), not numpy.cos(). The numpy versions of these
> functions are optimized for handling array arguments and are much
> slower than the math versions for scalar arguments.

Yes I understand that. I just want to stress that it was not a benchmark
(nor a critic) but a test to know if it was interesting to translate
directly an IDL code into python/numpy before trying to optimize it (I
know more python than IDL). I expected  to have approximatively the same
speed for both, was surprised by the result, and wanted to know if there
was an obvious reason besides the unoptimization for scalars.
```
1 Dec 2010 15:26

### Python versions for NumPy 1.5

```Does NumPy 1.5 work with Python 2.7 or Python 3.x?
```
1 Dec 2010 15:32

### Re: Python versions for NumPy 1.5

On Wed, Dec 1, 2010 at 7:26 AM, John Hornstein wrote:
Does NumPy 1.5 work with Python 2.7 or Python 3.x?

Yes, both. NumPy 1.5.1 fixes some small bugs and that is what you should use.

Chuck

```_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion <at> scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
```
1 Dec 2010 17:07

### Re: A faster median (Wirth's method)

<at> Keith Goodman

I think I figured it out. I believe something like the following will do what you want, iterating across one axis specially, so it can apply a median function along an axis. This code in particular is for calculating a moving average and seems to work (though I haven't checked my math). Let me know if you find any problems.

def ewma(a, d,  int axis = -1):

out = np.empty(a.shape, dtype)

cdef np.flatiter ita, ito

ita = np.PyArray_IterAllButAxis(a,   &axis)
ito = np.PyArray_IterAllButAxis(out, &axis)

cdef int i
cdef int axis_length = a.shape[axis]
cdef int a_axis_stride = a.strides[axis]/a.itemsize
cdef int o_axis_stride = out.strides[axis]/out.itemsize

cdef double avg = 0.0
cdef double weight = 1.0 - np.exp(-d)

while np.PyArray_ITER_NOTDONE(ita):

avg = 0.0
for i in range(axis_length):

avg += (<dtype_t*>np.PyArray_ITER_DATA (ita))[i * a_axis_stride ] * weight + avg * (1 - weight)
(<dtype_t*>np.PyArray_ITER_DATA (ito))[i * o_axis_stride ] = avg

np.PyArray_ITER_NEXT(ita)
np.PyArray_ITER_NEXT(ito)

return out

On Tue, Nov 30, 2010 at 9:46 PM, Felix Schlesinger wrote:
> > import numpy as np
> > cimport numpy as cnp
>
> > cdef cnp.float64_t namean(cnp.ndarray[cnp.float64_t, ndim=1] a):
> >    return np.nanmean(a)  # just a placeholder
>
> > is not allowed?  It works for me.  Is it a cython version thing?
> > (I've got 0.13),
>
> Oh, that's nice! I'm using 0.11.2. OK, time to upgrade.

Oh wow, does that mean that http://trac.cython.org/cython_trac/ticket/177
is fixed? I couldn't find anything in the release notes about that,
but it would be great news. Does the cdef function acquire and hold
the buffer?

Felix

_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion <at> scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

```_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion <at> scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
```
1 Dec 2010 20:16

### Re: broadcasting with numpy.interp

On Wed, Nov 24, 2010 at 3:16 PM, Friedrich Romstedt wrote:

2010/11/16 greg whittier <gregwh <at> gmail.com>:
> I'd like to be able to speed up the following code.
>
> def replace_dead(cube, dead):
>   # cube.shape == (320, 640, 1200)
>   # dead.shape == (320, 640)
>   # cube[i,j,:] are bad points to be replaced via interpolation if
> dead[i,j] == True
>
>    bands = np.arange(0, cube.shape[0])
>    for line in range(cube.shape[1]):
>        dead_bands = bands[dead[:, line] == True]
>        good_bands = bands[dead[:, line] == False]
>        for sample in range(cube.shape[2]):
>            # interp returns fp[0] for x < xp[0] and fp[-1] for x > xp[-1]
>            cube[dead_bands, line, sample] = \
>                np.interp(dead_bands,
>                          good_bands,
>                          cube[good_bands, line, sample])

I assume you just need *some* interpolation, not that specific one?
In that case, I'd suggest the following:

1)  Use a 2d interpolation, taking into account all nearest neighbours.
2)  For this, use a looped interpolation in this nearest-neighbour sense:
a)  Generate sums of all unmasked nearest-neighbour values
b)  Generate counts for the nearest neighbours present
c)  Replace the bad values by the sums divided by the count.
d)  Continue at (a) if there are bad values left

Bad values which are neighbouring each other (>= 3) need multiple
passes through the loop.  It should be pretty fast.

If this is what you have in mind, maybe we (or I) can make up some code.

Friedrich

Thanks so much for the response!  Sorry I didn't respond earlier.  I put it aside until I found time to try and understand part 2 of your response and forgot about it.  I'm not really looking for 2d interpolation at the moment, but I can see needing it in the future.  Right now, I just want to interpolate along one of the three axes.  I think what you're suggesting might work for 1d or 2d depending on how you find the nearest neighbors.  What routine would you use?  Also, when you say "unmasked" do you mean literally using masked arrays?

Thanks,
Greg

```_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion <at> scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
```
1 Dec 2010 20:18

### printoption to allow hexified floats?

```Hi Numpy folks,
When working with floats, I prefer to have exact string
representations in doctests and other reference-based testing; I find it
helps a lot to avoid chasing cross-platform differences that are really
about the string conversion rather than about numerical differences.
Since Python 2.6, the, the hex() method on floats has been available and
it gives an exact representation.  Is there any way to have Numpy arrays
of floats printed using this representation?  If not, would there be
interest in adding that?
On a somewhat related note, is there a table someplace which shows
which versions of Python are supported in each release of Numpy?  I
found an FAQ that mentioned 2.4 and 2.5, but since it didn't mention 2.6
or 2.7 (much less 3.1), I assume it's out of date.  This relates to the
above since it would be harder to support a new hex printoption for
Pythons before 2.6.
Thanks,
Ken B.
```
1 Dec 2010 20:47

### A Cython apply_along_axis function

```It's hard to write Cython code that can handle all dtypes and
arbitrary number of dimensions. The former is typically dealt with
using templates, but what do people do about the latter?

I'm trying to take baby steps towards writing an apply_along_axis
function that takes as input a cython function, numpy array, and axis.
I'm using the following numpy ticket as a guide but I'm really just
copying and pasting without understanding:

http://projects.scipy.org/numpy/attachment/ticket/1213/_selectmodule.pyx

Can anyone spot why I get a segfault on the call to nanmean_1d in
apply_along_axis?

import numpy as np
cimport numpy as np
import cython

cdef double NAN = <double> np.nan
ctypedef np.float64_t (*func_t)(void *buf, np.npy_intp size, np.npy_intp s)

def apply_along_axis(np.ndarray[np.float64_t, ndim=1] a, int axis):

cdef func_t nanmean_1d
cdef np.npy_intp stride, itemsize
cdef int ndim = a.ndim
cdef np.float64_t out

itemsize = <np.npy_intp> a.itemsize

if ndim == 1:
stride = a.strides[0] // itemsize # convert stride bytes --> items
out = nanmean_1d(a.data, a.shape[0], stride)
else:
raise ValueError("Not yet coded")

return out

cdef np.float64_t nanmean_1d(void *buf, np.npy_intp n, np.npy_intp s):
"nanmean of buffer."
cdef np.float64_t *a = <np.float64_t *> buf #
cdef np.npy_intp i, count = 0
cdef np.float64_t asum, ai
if s == 1:
for i in range(n):
ai = a[i]
if ai == ai:
asum += ai
count += 1
else:
for i in range(n):
ai = a[i*s]
if ai == ai:
asum += ai
count += 1
if count > 0:
return asum / count
else:
return NAN
```
1 Dec 2010 21:00

### Re: A Cython apply_along_axis function

```On 12/01/2010 08:47 PM, Keith Goodman wrote:
> It's hard to write Cython code that can handle all dtypes and
> arbitrary number of dimensions. The former is typically dealt with
> using templates, but what do people do about the latter?
>

What you typically do is to use the C-level iterator API. In fact
there's a recent thread on cython-users that does exactly this ("How can
I use PyArray_IterAllButAxis..."). Of course, make sure you take the
comments of that thread into account (!).

I feel that is easier to work with than what you do below. Not saying it
couldn't be easier, but it's not too bad once you get used to it.

Dag Sverre

> I'm trying to take baby steps towards writing an apply_along_axis
> function that takes as input a cython function, numpy array, and axis.
> I'm using the following numpy ticket as a guide but I'm really just
> copying and pasting without understanding:
>
> http://projects.scipy.org/numpy/attachment/ticket/1213/_selectmodule.pyx
>
> Can anyone spot why I get a segfault on the call to nanmean_1d in
> apply_along_axis?
>
> import numpy as np
> cimport numpy as np
> import cython
>
> cdef double NAN =<double>  np.nan
> ctypedef np.float64_t (*func_t)(void *buf, np.npy_intp size, np.npy_intp s)
>
> def apply_along_axis(np.ndarray[np.float64_t, ndim=1] a, int axis):
>
>      cdef func_t nanmean_1d
>      cdef np.npy_intp stride, itemsize
>      cdef int ndim = a.ndim
>      cdef np.float64_t out
>
>      itemsize =<np.npy_intp>  a.itemsize
>
>      if ndim == 1:
>          stride = a.strides[0] // itemsize # convert stride bytes -->  items
>          out = nanmean_1d(a.data, a.shape[0], stride)
>      else:
>          raise ValueError("Not yet coded")
>
>      return out
>
> cdef np.float64_t nanmean_1d(void *buf, np.npy_intp n, np.npy_intp s):
>      "nanmean of buffer."
>      cdef np.float64_t *a =<np.float64_t *>  buf #
>      cdef np.npy_intp i, count = 0
>      cdef np.float64_t asum, ai
>      if s == 1:
>          for i in range(n):
>              ai = a[i]
>              if ai == ai:
>                  asum += ai
>                  count += 1
>      else:
>          for i in range(n):
>              ai = a[i*s]
>              if ai == ai:
>                  asum += ai
>                  count += 1
>      if count>  0:
>          return asum / count
>      else:
>          return NAN
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion <at> scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
```

Gmane