headers
John Salvatier | 1 Dec 2010 04:22
Gravatar

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 <kwgoodman <at> gmail.com> 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
Permalink | Reply | 0 comments |
headers
Felix Schlesinger | 1 Dec 2010 06:46
Favicon

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
Permalink | Reply | 1 comment |
headers
Jean-Luc Menut | 1 Dec 2010 11:23
Picon
Favicon

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.
Permalink | Reply | 0 comments |
headers
John Hornstein | 1 Dec 2010 15:26
Picon
Picon
Favicon

Python versions for NumPy 1.5

Does NumPy 1.5 work with Python 2.7 or Python 3.x?
Permalink | Reply | 0 comments |
headers
Charles R Harris | 1 Dec 2010 15:32
Picon

Re: Python versions for NumPy 1.5



On Wed, Dec 1, 2010 at 7:26 AM, John Hornstein <John.Hornstein <at> nrl.navy.mil> 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
Permalink | Reply | 0 comments |
headers
John Salvatier | 1 Dec 2010 17:07
Gravatar

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 <schlesin <at> cshl.edu> 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
Permalink | Reply | 0 comments |
headers
greg whittier | 1 Dec 2010 20:16
Picon

Re: broadcasting with numpy.interp

On Wed, Nov 24, 2010 at 3:16 PM, Friedrich Romstedt <friedrichromstedt <at> gmail.com> 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
Permalink | Reply | 1 comment |
headers
Ken Basye | 1 Dec 2010 20:18
Favicon

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.
Permalink | Reply | 2 comments |
headers
Keith Goodman | 1 Dec 2010 20:47
Picon
Gravatar

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
Permalink | Reply | 14 comments |
headers
Dag Sverre Seljebotn | 1 Dec 2010 21:00
Picon
Picon

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
>
Permalink | Reply | 0 comments |

Gmane