[Numpy-discussion] numpy ndarray questions

[Jochen]

On Tue, 2009-01-27 at 12:37 +0100, Sturla Molden wrote:
> On 1/27/2009 1:26 AM, Jochen wrote:
> 
> > a = fftw3.AlignedArray(1024,complex)
> > 
> > a = a+1
> 
> = used this way is not assignment, it is name binding.
> 
> It is easy to use function's like fftw_malloc with NumPy:
> 
> 
> import ctypes
> import numpy
> 
> fftw_malloc = ctypes.cdll.fftw.fftw_malloc
> fftw_malloc.argtypes = [ctypes.c_ulong,]
> fftw_malloc.restype = ctypes.c_ulong
> 
> def aligned_array(N, dtype):
>      d = dtype()
>      address = fftw_malloc(N * d.nbytes) # restype = ctypes.c_ulong
>      if (address = 0): raise MemoryError, 'fftw_malloc returned NULL'
>      class Dummy(object): pass
>      d = Dummy()
>      d.__array_interface__ = {
>          'data' = (address, False),
>          'typestr' : dtype.str,
>          'descr' : dtype.descr,
>          'shape' : shape,
>          'strides' : strides,
>          'version' : 3
>      }
>      return numpy.asarray(d)
> 

I actually do it slightly different, because I want to free the memory
using fftw_free. So I'm subclassing ndarray and in the __new__ function
of the subclass I create a buffer object from fftw_malloc using
PyBuffer_FromReadWriteMemory and pass that to ndarray.__new__. In
__del__ I check if the .base is a buffer object and do an fftw_free.

> 
> If you have to check for a particular alignment before calling fftw, 
> that is trivial as well:
> 
> 
> def is_aligned(array, boundary):
>      address = array.__array_interface__['data'][0]
>      return not(address % boundary)
> 
> 

Yes I knew that, btw do you know of any systems which need alignment to
boundaries other than 16byte for simd operations?

> 
> 
> > there a way that I get a different object type?  Or even better is there
> > a way to prevent operations like a=a+1 or make them automatically
> > in-place operations?
> 
> a = a + 1 # rebinds the name 'a' to another array
> 
> a[:] = a + 1 # fills a with the result of a + 1
>
I knew about this one, this is what I have been doing.
> 
> This has to do with Python syntax, not NumPy per se. You cannot overload 
> the behaviour of Python's name binding operator (=).
> 
Yes I actually thought about it some more yesterday when going home and
realised that this would really not be possible. I guess I just have to
document it clearly so that people don't run into it.

Cheers
Jochen
> 
> Sturla Molden
> 
> 
> 
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion at scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion