Re: [Numpy-discussion] any interest in including asecond-ordergradient?

2008-10-30 Thread Stéfan van der Walt
Hi Fernando,

Thanks for your input.

2008/10/29 Fernando Perez [EMAIL PROTECTED]:
 I think it's fine to ask for functions that compute higher order
 derivatives of n-d arrays: we already have diff(), which operates on a
 single direction, and a hessian could make sense (with the caveats
 David points out).   But with higher order derivatives there are many
 more combinations to worry about, and I really think it's a bad idea
 to lump those issues into the definition of gradient, which was a
 perfectly unambiguous object up until this point.

Maybe we should focus on writing a decent 'deriv' function then.  I
know Konrad Hinsen's Scientific had a derivatives package
(Scientific.Functions.Derivatives) that implemented automatic
differentiation:

http://en.wikipedia.org/wiki/Automatic_differentiation

Cheers
Stéfan
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Adding the ability to clone a few fields from a data-type

2008-10-30 Thread Stéfan van der Walt
2008/10/30 Robert Kern [EMAIL PROTECTED]:
 Provide a kind of inaccessible and invisible dtype for implementing
 dummy fields. This is useful in other places like file parsing. At the
 same time, implement a function that uses this capability to make
 views with a subset of the fields of a structured array. I'm not sure
 that people need an API for replacing the fields of a dtype like this.

I have thought about adding such an invisible dtype before
(incidentally, it was while parsing files!).  I would be interested in
exploring this route further.

Regards
Stéfan
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] passing a C array to embedded Python from C code

2008-10-30 Thread Matthieu Brucher
 Does this approach make sense?  Is there a better way to go about it?
 Maybe calling a custom module from the Python code that does the C
 array to NumPy translation using Cython/pyrex/swig/etc.  Would it be
 possible to use the same C arrays from here without copying them?

Hi,

Your case seems to fit the array interface. The goal is to create a C
structure with some additional information that Numpy can understand,
and then your array will be treated as a Numpy array. If you can
follow a French tutorial, you can go on
http://matthieu-brucher.developpez.com/tutoriels/python/swig-numpy/#LV
to have a skeletton for your issue.

Matthieu
-- 
Information System Engineer, Ph.D.
Website: http://matthieu-brucher.developpez.com/
Blogs: http://matt.eifelle.com and http://blog.developpez.com/?blog=92
LinkedIn: http://www.linkedin.com/in/matthieubrucher
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] Numpy matrix multiplication slow even though ATLAS linked

2008-10-30 Thread Jan-Willem van de Meent
Dear all,

This is my first post to this list. I am having perfomance issues with with 
numpy/atlas. Doing dot(a,a) for a 2000x2000 matrix takes  about 1m40s, even 
though numpy is appears to link to my atlas libraries:

I did a quick benchmark by running the following script: 

#! /usr/bin/env python
import numpy
import time

try:
   import numpy.core._dotblas
   print 'Using ATLAS:'
except ImportError:
   print 'No ATLAS:'

t = time.time()
x = numpy.random.random((1000,1000))
y = numpy.random.random((1000,1000))
z = numpy.dot(x, y)

print time.time()-t

My laptop is a Dell D620 Core Duo T2300 1.66 Ghz, running Archlinux with GCC 
4.3.2, atlas 3.8.2, python 2.5.2 and numpy 1.2.1. Output of the script above 
is:

Using ATLAS:
7.99549412727

A department desktop PC, Pentium D 3.00 Ghz, running Scientific Linux, with 
GCC 4.1.2, atlas 3.7.30, python 2.5.1 and numpy 1.1.0, runs this test 24 
times faster:

Using ATLAS:
0.337520122528

So even though _dotblas.so exists, matrix multiplication appears to run at 
pretty much the same speed as if atlas were not  available. Running ldd on 
_dotblas.so suggests  that numpy is indeed linking to the atlas libs:

ldd /usr/lib/python2.5/site-packages/numpy/core/_dotblas.so
linux-gate.so.1 =  (0xb7fcf000)
libatlas.so = /usr/lib/libatlas.so (0xb7cb5000)
liblapack.so = /usr/lib/liblapack.so (0xb77ab000)
libcblas.so = /usr/lib/libcblas.so (0xb778b000)
libf77blas.so = /usr/lib/libf77blas.so (0xb776f000)
libpython2.5.so.1.0 = /usr/lib/libpython2.5.so.1.0 (0xb763)
libpthread.so.0 = /lib/libpthread.so.0 (0xb7618000)
libc.so.6 = /lib/libc.so.6 (0xb74d6000)
libm.so.6 = /lib/libm.so.6 (0xb74b)
libgfortran.so.3 = /usr/lib/libgfortran.so.3 (0xb73ff000)
libdl.so.2 = /lib/libdl.so.2 (0xb73fb000)
libutil.so.1 = /lib/libutil.so.1 (0xb73f6000)
/lib/ld-linux.so.2 (0xb7fd)

Something appears to be going wrong in the compilation process. In particular, 
the compile process appears unable to determine the version of the atlas libs 
installed (I get NO_ATLAS_INFO). I've pasted some snippets from the build 
output which may be relevant below.

Can anyone help me figure this out? 

Many thanks in advance,
Jan-Willem

-- 
Jan-Willem van de Meent
Research Student
Goldstein Lab, DAMTP
University of Cambridge

--
Snippets from setup.py output
...

atlas_threads_info:
Setting PTATLAS=ATLAS
  libraries lapack_atlas not found in /usr/lib
numpy.distutils.system_info.atlas_threads_info
Setting PTATLAS=ATLAS
Setting PTATLAS=ATLAS
  FOUND:
libraries = 
['lapack', 'atlas', 'lapack', 'cblas', 'f77blas', 'ptcblas', 'ptf77blas']
library_dirs = ['/usr/lib']
language = f77
include_dirs = ['/usr/include']

...

customize Gnu95FCompiler using config
compiling '_configtest.c':

/* This file is generated from numpy/distutils/system_info.py */
void ATL_buildinfo(void);
int main(void) {
  ATL_buildinfo();
  return 0;
}
C compiler: 
gcc -pthread -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -march=i686 
-mtune=generic -O2 -pipe -fPIC

compile options: '-c'
gcc: _configtest.c
gcc -pthread 
_configtest.o -L/usr/lib -llapack -latlas -llapack -lcblas -lf77blas -lptcblas 
-lptf77blas -o 
_configtest
/usr/bin/ld: _configtest: hidden symbol `__powidf2' 
in /usr/lib/gcc/i686-pc-linux-gnu/4.3.2/libgcc.a(_powidf2.o) is referenced by 
DSO
/usr/bin/ld: final link failed: Nonrepresentable section on output
collect2: ld returned 1 exit status
/usr/bin/ld: _configtest: hidden symbol `__powidf2' 
in /usr/lib/gcc/i686-pc-linux-gnu/4.3.2/libgcc.a(_powidf2.o) is referenced by 
DSO
/usr/bin/ld: final link failed: Nonrepresentable section on output
collect2: ld returned 1 exit status
failure.
removing: _configtest.c _configtest.o
Status: 255
Output: 
  FOUND:
libraries = 
['lapack', 'atlas', 'lapack', 'cblas', 'f77blas', 'ptcblas', 'ptf77blas']
library_dirs = ['/usr/lib']
language = f77
define_macros = [('NO_ATLAS_INFO', 2)]
include_dirs = ['/usr/include']
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Adding the ability to clone a few fields from a data-type

2008-10-30 Thread Francesc Alted
A Thursday 30 October 2008, Robert Kern escrigué:
 On Wed, Oct 29, 2008 at 19:05, Travis E. Oliphant

 [EMAIL PROTECTED] wrote:
  Hi all,
 
  I'd like to add to NumPy the ability to clone a data-type object so
  that only a view fields are copied over but that it retains the
  same total size.
 
  This would allow, for example, the ability to select out a few
  records from a structured array using
 
  subarr = arr.view(cloned_dtype)
 
  Right now, it is hard to do this because you have to at least add a
  dummy field at the end.  A simple method on the dtype class
  (fromfields or something) would be easy to add.

 I'm not sure what this accomplishes. Would the dummy fields that fill
 in the space be inaccessible? E.g. tuple(subarr[i,j,k]) gives a tuple
 with no numpy.void scalars? That would be a novel feature, but I'm
 not

 sure it fits the problem. On the contrary:
  It was thought in the past to do this with indexing
 
  arr['field1', 'field2']
 
  And that would still be possible (and mostly implemented) if this
  feature is added.

 This appears more like the interface that people want. Except that I
 think people were thinking that it would follow fancy indexing
 syntax:

   arr[['field1', 'field2']]

I've thought about that too.  That would be a great thing to have, IMO.


 I guess there are two ways to implement this. One is to make a new
 array that just contains the desired fields. Another is to make a
 view that just points to the desired fields in the original array
 provided that we have a new feature for inaccessible dummy fields.
 One point for the former approach is that it is closer to fancy
 indexing which must always make a copy. The latter approach breaks
 that connection.

Yeah.  I'd vote for avoid the copy.

 OTOH, now that I think about it, I don't think there is really any
 coherent way to mix field selection with any other indexing
 operations. At least, not within the same brackets. Hmm. So maybe the
 link to fancy indexing can be ignored as, ahem, fanciful.

Well, one can always check that fields in the fancy list are either 
strings (map to name fields) or integers (map to positional fields).  
However, I'm not sure if this check would be too expensive.

 Overall, I guess, I would present the feature slightly differently.
 Provide a kind of inaccessible and invisible dtype for implementing
 dummy fields. This is useful in other places like file parsing. At
 the same time, implement a function that uses this capability to make
 views with a subset of the fields of a structured array. I'm not sure
 that people need an API for replacing the fields of a dtype like
 this.

Mmh, not sure on what you are proposing there.  You mean something like:

In [21]: t = numpy.dtype([('f0','i4'),('f1', 'f8'), ('f2', 'S20')])

In [22]: nt = t.astype(['f2', 'f0'])

In [23]: ra = numpy.zeros(10, dtype=t)

In [24]: nra = ra.view(nt)

In [25]: ra
Out[25]:
array([(0, 0.0, ''), (0, 0.0, ''), (0, 0.0, ''), (0, 0.0, ''),
   (0, 0.0, ''), (0, 0.0, ''), (0, 0.0, ''), (0, 0.0, ''),
   (0, 0.0, ''), (0, 0.0, '')],
  dtype=[('f0', 'i4'), ('f1', 'f8'), ('f2', '|S20')])

In [26]: nra
Out[26]:
array([('', 0), ('', 0), ('', 0), ('', 0), ('', 0), ('', 0), ('', 0),
   ('', 0), ('', 0), ('', 0)],
  dtype=[('f2', '|S20'), ('f0', 'i4')])

?

In that case, that would be a great feature to add.

-- 
Francesc Alted
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] numpy array change notifier?

2008-10-30 Thread Christan K .
Erik Tollerud erik.tollerud at gmail.com writes:

 
 Is there any straightforward way of notifying on change of a numpy
 array that leaves the numpy arrays still efficient?
 

You can instantiate the following ndarray derived class with a callable 
argument. Any change to the data will call the callback function.

In [5]: def notify(arg):
   ...: print 'array changed'
   ...:
   ...:

In [6]: x = cbarray([1,2,3,4],cb=notify)

In [7]: x[2] = -1
array changed

alternatively I have used the pubsub module available as part of wxPython (see 
the commented line in _notify())

class cbarray(N.ndarray):
def __new__(subtype, data, cb=None, dtype=None, copy=False):
subtype.__defaultcb = cb

if copy:
data = N.array(data,dtype=dtype)
else:
data = N.asarray(data,dtype=dtype)

data = data.view(subtype)
return data

def _notify(self):
if self.cb is not None:
self.cb()
#Publisher().sendMessage(('changed'))

def _get_shape(self):
return super(cbarray, self).shape
shape = property(_get_shape)

def __setitem__(self, item, val):
N.ndarray.__setitem__(self, item, val)
self._notify()

def __array_finalize__(self,obj):
if not hasattr(self, cb):
# The object does not already have a `.cb` attribute
self.cb = getattr(obj,'cb',self.__defaultcb)

def __reduce__(self):
object_state = list(N.ndarray.__reduce__(self))
subclass_state = (self.cb,)
object_state[2] = (object_state[2],subclass_state)
return tuple(object_state)

def __setstate__(self,state):
nd_state, own_state = state
N.ndarray.__setstate__(self,nd_state)

cb, = own_state
self.cb = cb

Hope thqat helps, Christian


___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Adding the ability to clone a few fields from a data-type

2008-10-30 Thread Travis E. Oliphant

 I'm not sure what this accomplishes. Would the dummy fields that fill
 in the space be inaccessible? E.g. tuple(subarr[i,j,k]) gives a tuple
 with no numpy.void scalars? That would be a novel feature, but I'm not
 sure it fits the problem. On the contrary:
   

Yes, that was the idea.   You can do it now, but only in C.   The real 
problem right now from my point of view is that there is no way to tell 
the dtype constructor to pad the itemsize to x bytes.If that were 
changed, then many things would be possible. 

 OTOH, now that I think about it, I don't think there is really any
 coherent way to mix field selection with any other indexing
 operations. At least, not within the same brackets. Hmm. So maybe the
 link to fancy indexing can be ignored as, ahem, fanciful.
   
Yeah,  I was wondering how to do it well, myself, and couldn't come up 
with anything which is why I went the .view route with another dtype.   

By inaccessible and invisible dtype do you mean something like the 
basic built-in void data type, but which doesn't try to report itself 
when the dtype prints?

That sounds interesting but I'm not sure it's necessary because the 
field specification can already skip bytes (just not bytes at the end 
--- which is what I would like to fix).Perhaps what is needed is a 
pseudo-dtype (something like 'c' compared to 'S1') which doesn't 
actually create a new dtype but which is handled differently when the 
dtype is created with the [('field1', type), ('field2', type2)] 
approach.   Specifically, it doesn't add an entry to the fields 
dictionary nor an entry to the names but does affect the itemsize of the 
element (and the offset of follow-on fields).

So, let's assume the character is 'v':

If we have an array with underlying dtype:

od = [('date', 'S10'), ('high', 'f4'), ('low', 'f4'), ('close', 'f4'), 
('volume', 'i4')]

Then, we could define a new dtype

nd = [('date', 'S10'), ('', 'v8'), ('close', 'f4'), ('', 'v4')]

and  arr.view(nd)   would provide a view of the array where element 
selection would be a tuple with just the date and close elements but the 
itemsize would be exactly the same but nd.names would be ['date', 'close']

I like this approach.  It impacts the API the very least but provides 
the desired functionality.

-Travis



-Travis



 Overall, I guess, I would present the feature slightly differently.
 Provide a kind of inaccessible and invisible dtype for implementing
 dummy fields. This is useful in other places like file parsing. At the
 same time, implement a function that uses this capability to make
 views with a subset of the fields of a structured array. I'm not sure
 that people need an API for replacing the fields of a dtype like this.

   

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] numpy array change notifier?

2008-10-30 Thread Steve Schmerler
On Oct 27 16:43 -0400, Pierre GM  wrote:
 
 Erik, may be you could try the trick presented here : 
 http://www.scipy.org/Subclasses
 in the __array_wrap__ section.

Stupid question: How do I find pages like
http://www.scipy.org/Subclasses on scipy.org?

I can find them with the search function, but only when I know that it's
there :)

best,
steve
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] numpy array change notifier?

2008-10-30 Thread Pierre GM
Steve,
Right, there's not a lot of visibility for this one. You could try 
google 'numpy subclass'. there's also an entry on subclassing in the numpy 
user guide (through docs.scipy.org).
Cheers
P.

On Thursday 30 October 2008 09:57:03 Steve Schmerler wrote:
 On Oct 27 16:43 -0400, Pierre GM  wrote:
  Erik, may be you could try the trick presented here :
  http://www.scipy.org/Subclasses
  in the __array_wrap__ section.

 Stupid question: How do I find pages like
 http://www.scipy.org/Subclasses on scipy.org?

 I can find them with the search function, but only when I know that it's
 there :)

 best,
 steve
 ___
 Numpy-discussion mailing list
 Numpy-discussion@scipy.org
 http://projects.scipy.org/mailman/listinfo/numpy-discussion


___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] numpy array change notifier?

2008-10-30 Thread Steve Schmerler
On Oct 30 10:06 -0400, Pierre GM  wrote:
 Steve,
 Right, there's not a lot of visibility for this one. 

Yes, wouldn't it be better placed in the Cookbook section (or any other
suitable place that's visible to anyone entering scipy.org)?  Other
pages are unfortunately also not visible, like
http://www.scipy.org/EricsBroadcastingDoc .

best,
steve
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] Automatic differentiation (was Re: second-order gradient)

2008-10-30 Thread Rob Clewley
 Maybe we should focus on writing a decent 'deriv' function then.  I
 know Konrad Hinsen's Scientific had a derivatives package
 (Scientific.Functions.Derivatives) that implemented automatic
 differentiation:

 http://en.wikipedia.org/wiki/Automatic_differentiation

That would be great, but wouldn't that be best suited as a utility
requiring Sympy? You'll want to take advantage of all sorts of
symbolic classes, especially for any source code transformation
approach. IMO Hinsen's implementation isn't a very efficient or
attractive solution to AD given the great existing C/C++ codes out
there. Maybe we should be looking to provide a python interface to an
existing open source package such as ADOL-C, but I'm all in favour of
a new pure python approach too. What would be perfect is to have a
single interface to a python AD package that would support a faster
implementation if the user wished to install a C/C++ package,
otherwise would default to a pure python equivalent.

-Rob
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] any interest in including asecond-ordergradient?

2008-10-30 Thread Rob Clewley
 2008/10/29 Fernando Perez [EMAIL PROTECTED]:
 I think it's fine to ask for functions that compute higher order
 derivatives of n-d arrays: we already have diff(), which operates on a
 single direction, and a hessian could make sense (with the caveats
 David points out).   But with higher order derivatives there are many
 more combinations to worry about, and I really think it's a bad idea
 to lump those issues into the definition of gradient, which was a
 perfectly unambiguous object up until this point.

I'm basically in favour of Fernando's suggestion to keep gradient
simple and add a hessian function. Higher numerical derivatives from
data aren't very reliable anyway. You're much better off interpolating
with a polynomial and then differentiating that.

 Maybe we should focus on writing a decent 'deriv' function then.  I
 know Konrad Hinsen's Scientific had a derivatives package
 (Scientific.Functions.Derivatives) that implemented automatic
 differentiation:

Improving the support for a gradient of array data is an entirely
independent project in my mind - but I like this idea and I replied in
a new thread.

-Rob
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Adding the ability to clone a few fields from a data-type

2008-10-30 Thread Robert Kern
On Thu, Oct 30, 2008 at 08:27, Travis E. Oliphant
[EMAIL PROTECTED] wrote:

 I'm not sure what this accomplishes. Would the dummy fields that fill
 in the space be inaccessible? E.g. tuple(subarr[i,j,k]) gives a tuple
 with no numpy.void scalars? That would be a novel feature, but I'm not
 sure it fits the problem. On the contrary:


 Yes, that was the idea.   You can do it now, but only in C.   The real
 problem right now from my point of view is that there is no way to tell
 the dtype constructor to pad the itemsize to x bytes.If that were
 changed, then many things would be possible.

 OTOH, now that I think about it, I don't think there is really any
 coherent way to mix field selection with any other indexing
 operations. At least, not within the same brackets. Hmm. So maybe the
 link to fancy indexing can be ignored as, ahem, fanciful.

 Yeah,  I was wondering how to do it well, myself, and couldn't come up
 with anything which is why I went the .view route with another dtype.

 By inaccessible and invisible dtype do you mean something like the
 basic built-in void data type, but which doesn't try to report itself
 when the dtype prints?

The field doesn't report itself when the *values* print, is what I'm
concerned with. The dtype should display the dummy fields such that
repr() can accurately reconstruct the dtype.

 That sounds interesting but I'm not sure it's necessary because the
 field specification can already skip bytes (just not bytes at the end
 --- which is what I would like to fix).Perhaps what is needed is a
 pseudo-dtype (something like 'c' compared to 'S1') which doesn't
 actually create a new dtype but which is handled differently when the
 dtype is created with the [('field1', type), ('field2', type2)]
 approach.   Specifically, it doesn't add an entry to the fields
 dictionary nor an entry to the names but does affect the itemsize of the
 element (and the offset of follow-on fields).

 So, let's assume the character is 'v':

 If we have an array with underlying dtype:

 od = [('date', 'S10'), ('high', 'f4'), ('low', 'f4'), ('close', 'f4'),
 ('volume', 'i4')]

 Then, we could define a new dtype

 nd = [('date', 'S10'), ('', 'v8'), ('close', 'f4'), ('', 'v4')]

To do this, we would also have to fix the current behavior of
converting ''s to 'f0', 'f1', etc., when these are passed to the
dtype() constructor.

 and  arr.view(nd)   would provide a view of the array where element
 selection would be a tuple with just the date and close elements but the
 itemsize would be exactly the same but nd.names would be ['date', 'close']

 I like this approach.  It impacts the API the very least but provides
 the desired functionality.

-- 
Robert Kern

I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth.
  -- Umberto Eco
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Adding the ability to clone a few fields from a data-type

2008-10-30 Thread Robert Kern
On Thu, Oct 30, 2008 at 04:33, Francesc Alted [EMAIL PROTECTED] wrote:
 A Thursday 30 October 2008, Robert Kern escrigué:
 On Wed, Oct 29, 2008 at 19:05, Travis E. Oliphant

 [EMAIL PROTECTED] wrote:
  Hi all,
 
  I'd like to add to NumPy the ability to clone a data-type object so
  that only a view fields are copied over but that it retains the
  same total size.
 
  This would allow, for example, the ability to select out a few
  records from a structured array using
 
  subarr = arr.view(cloned_dtype)
 
  Right now, it is hard to do this because you have to at least add a
  dummy field at the end.  A simple method on the dtype class
  (fromfields or something) would be easy to add.

 I'm not sure what this accomplishes. Would the dummy fields that fill
 in the space be inaccessible? E.g. tuple(subarr[i,j,k]) gives a tuple
 with no numpy.void scalars? That would be a novel feature, but I'm
 not

 sure it fits the problem. On the contrary:
  It was thought in the past to do this with indexing
 
  arr['field1', 'field2']
 
  And that would still be possible (and mostly implemented) if this
  feature is added.

 This appears more like the interface that people want. Except that I
 think people were thinking that it would follow fancy indexing
 syntax:

   arr[['field1', 'field2']]

 I've thought about that too.  That would be a great thing to have, IMO.

 I guess there are two ways to implement this. One is to make a new
 array that just contains the desired fields. Another is to make a
 view that just points to the desired fields in the original array
 provided that we have a new feature for inaccessible dummy fields.
 One point for the former approach is that it is closer to fancy
 indexing which must always make a copy. The latter approach breaks
 that connection.

 Yeah.  I'd vote for avoid the copy.

 OTOH, now that I think about it, I don't think there is really any
 coherent way to mix field selection with any other indexing
 operations. At least, not within the same brackets. Hmm. So maybe the
 link to fancy indexing can be ignored as, ahem, fanciful.

 Well, one can always check that fields in the fancy list are either
 strings (map to name fields) or integers (map to positional fields).
 However, I'm not sure if this check would be too expensive.

That's not my concern. The problem is that the field-indexing applies
to the entire array, not just an axis. So what would the following
mean?

  a[['foo', 'bar'], [1,2,3]]

Compared to

  a[[5,8,10], [1,2,3]]

 Overall, I guess, I would present the feature slightly differently.
 Provide a kind of inaccessible and invisible dtype for implementing
 dummy fields. This is useful in other places like file parsing. At
 the same time, implement a function that uses this capability to make
 views with a subset of the fields of a structured array. I'm not sure
 that people need an API for replacing the fields of a dtype like
 this.

 Mmh, not sure on what you are proposing there.  You mean something like:

 In [21]: t = numpy.dtype([('f0','i4'),('f1', 'f8'), ('f2', 'S20')])

 In [22]: nt = t.astype(['f2', 'f0'])

 In [23]: ra = numpy.zeros(10, dtype=t)

 In [24]: nra = ra.view(nt)

 In [25]: ra
 Out[25]:
 array([(0, 0.0, ''), (0, 0.0, ''), (0, 0.0, ''), (0, 0.0, ''),
   (0, 0.0, ''), (0, 0.0, ''), (0, 0.0, ''), (0, 0.0, ''),
   (0, 0.0, ''), (0, 0.0, '')],
  dtype=[('f0', 'i4'), ('f1', 'f8'), ('f2', '|S20')])

 In [26]: nra
 Out[26]:
 array([('', 0), ('', 0), ('', 0), ('', 0), ('', 0), ('', 0), ('', 0),
   ('', 0), ('', 0), ('', 0)],
  dtype=[('f2', '|S20'), ('f0', 'i4')])

 ?

 In that case, that would be a great feature to add.

That's what Travis is proposing. I would like to see a function that
does this (however it is implemented under the covers):

  nra = subset_fields(ra, ['f0', 'f2'])

With the view, I don't think you can reorder the fields as in your example.

-- 
Robert Kern

I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth.
  -- Umberto Eco
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Automatic differentiation (was Re: second-order gradient)

2008-10-30 Thread Stéfan van der Walt
2008/10/30 Rob Clewley [EMAIL PROTECTED]:
 http://en.wikipedia.org/wiki/Automatic_differentiation

 That would be great, but wouldn't that be best suited as a utility
 requiring Sympy? You'll want to take advantage of all sorts of
 symbolic classes, especially for any source code transformation
 approach. IMO Hinsen's implementation isn't a very efficient or
 attractive solution to AD given the great existing C/C++ codes out
 there. Maybe we should be looking to provide a python interface to an
 existing open source package such as ADOL-C, but I'm all in favour of
 a new pure python approach too. What would be perfect is to have a
 single interface to a python AD package that would support a faster
 implementation if the user wished to install a C/C++ package,
 otherwise would default to a pure python equivalent.

In your experience, is this functionality enough to start a separate
package, or should we try to include it somewhere else?  Otherwise we
could think of a new SciKit.

Regards
Stéfan
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Automatic differentiation (was Re: second-order gradient)

2008-10-30 Thread Rob Clewley
 In your experience, is this functionality enough to start a separate
 package, or should we try to include it somewhere else?  Otherwise we
 could think of a new SciKit.

I confess to knowing no details about scikits so I don't know what the
difference really is between a new package and a scikit. To do this
properly you'd end up with a sizable body of code, and the potential
dependency on Sympy would also suggest making it somewhat separate.
I'd defer to others on that point, although I don't really see what
other package it would naturally fit with because AD has multiple
applications.

-Rob
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Adding the ability to clone a few fields from a data-type

2008-10-30 Thread Francesc Alted
A Thursday 30 October 2008, Robert Kern escrigué:
[clip]
  OTOH, now that I think about it, I don't think there is really any
  coherent way to mix field selection with any other indexing
  operations. At least, not within the same brackets. Hmm. So maybe
  the link to fancy indexing can be ignored as, ahem, fanciful.
 
  Well, one can always check that fields in the fancy list are either
  strings (map to name fields) or integers (map to positional
  fields). However, I'm not sure if this check would be too
  expensive.

 That's not my concern. The problem is that the field-indexing applies
 to the entire array, not just an axis. So what would the following
 mean?

   a[['foo', 'bar'], [1,2,3]]

 Compared to

   a[[5,8,10], [1,2,3]]

Well, as I see them, fields are like another axis, just that it is 
always the leading one.  In order to cope with them we could use a 
generalization of what it works already:

In [15]: ra = numpy.zeros((3,4), i4,f4)

In [16]: ra['f1'][[1,2],[0,3]]  # this already works
Out[16]: array([ 0.,  0.], dtype=float32)

In [17]: ra[['f1','f2']][[1,2],[0,3]]   # this could be make to work
Out[17]:
array([(0, 0.0), (0, 0.0)],
  dtype=[('f0', 'i4'), ('f1', 'f4')])

  Overall, I guess, I would present the feature slightly
  differently. Provide a kind of inaccessible and invisible dtype
  for implementing dummy fields. This is useful in other places like
  file parsing. At the same time, implement a function that uses
  this capability to make views with a subset of the fields of a
  structured array. I'm not sure that people need an API for
  replacing the fields of a dtype like this.
 
  Mmh, not sure on what you are proposing there.  You mean something
  like:
 
  In [21]: t = numpy.dtype([('f0','i4'),('f1', 'f8'), ('f2', 'S20')])
 
  In [22]: nt = t.astype(['f2', 'f0'])
 
  In [23]: ra = numpy.zeros(10, dtype=t)
 
  In [24]: nra = ra.view(nt)
 
  In [25]: ra
  Out[25]:
  array([(0, 0.0, ''), (0, 0.0, ''), (0, 0.0, ''), (0, 0.0, ''),
(0, 0.0, ''), (0, 0.0, ''), (0, 0.0, ''), (0, 0.0, ''),
(0, 0.0, ''), (0, 0.0, '')],
   dtype=[('f0', 'i4'), ('f1', 'f8'), ('f2', '|S20')])
 
  In [26]: nra
  Out[26]:
  array([('', 0), ('', 0), ('', 0), ('', 0), ('', 0), ('', 0), ('',
  0), ('', 0), ('', 0), ('', 0)],
   dtype=[('f2', '|S20'), ('f0', 'i4')])
 
  ?
 
  In that case, that would be a great feature to add.

 That's what Travis is proposing. I would like to see a function that
 does this (however it is implemented under the covers):

   nra = subset_fields(ra, ['f0', 'f2'])

Interesting.

 With the view, I don't think you can reorder the fields as in your
 example.

That's a pity.  Providing a dtype with the notion of an internal reorder 
can be very powerful in some situations.  But I guess that implementing 
this would be complicated.

-- 
Francesc Alted
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Adding the ability to clone a few fields from a data-type

2008-10-30 Thread Bruce Southey
Francesc Alted wrote:
 A Thursday 30 October 2008, Robert Kern escrigué:
 [clip]
   
 OTOH, now that I think about it, I don't think there is really any
 coherent way to mix field selection with any other indexing
 operations. At least, not within the same brackets. Hmm. So maybe
 the link to fancy indexing can be ignored as, ahem, fanciful.
 
 Well, one can always check that fields in the fancy list are either
 strings (map to name fields) or integers (map to positional
 fields). However, I'm not sure if this check would be too
 expensive.
   
 That's not my concern. The problem is that the field-indexing applies
 to the entire array, not just an axis. So what would the following
 mean?

   a[['foo', 'bar'], [1,2,3]]

 Compared to

   a[[5,8,10], [1,2,3]]
 

 Well, as I see them, fields are like another axis, just that it is 
 always the leading one.  In order to cope with them we could use a 
 generalization of what it works already:

 In [15]: ra = numpy.zeros((3,4), i4,f4)

 In [16]: ra['f1'][[1,2],[0,3]]  # this already works
 Out[16]: array([ 0.,  0.], dtype=float32)

 In [17]: ra[['f1','f2']][[1,2],[0,3]]   # this could be make to work
 Out[17]:
 array([(0, 0.0), (0, 0.0)],
   dtype=[('f0', 'i4'), ('f1', 'f4')])

   
 Overall, I guess, I would present the feature slightly
 differently. Provide a kind of inaccessible and invisible dtype
 for implementing dummy fields. This is useful in other places like
 file parsing. At the same time, implement a function that uses
 this capability to make views with a subset of the fields of a
 structured array. I'm not sure that people need an API for
 replacing the fields of a dtype like this.
 
 Mmh, not sure on what you are proposing there.  You mean something
 like:

 In [21]: t = numpy.dtype([('f0','i4'),('f1', 'f8'), ('f2', 'S20')])

 In [22]: nt = t.astype(['f2', 'f0'])

 In [23]: ra = numpy.zeros(10, dtype=t)

 In [24]: nra = ra.view(nt)

 In [25]: ra
 Out[25]:
 array([(0, 0.0, ''), (0, 0.0, ''), (0, 0.0, ''), (0, 0.0, ''),
   (0, 0.0, ''), (0, 0.0, ''), (0, 0.0, ''), (0, 0.0, ''),
   (0, 0.0, ''), (0, 0.0, '')],
  dtype=[('f0', 'i4'), ('f1', 'f8'), ('f2', '|S20')])

 In [26]: nra
 Out[26]:
 array([('', 0), ('', 0), ('', 0), ('', 0), ('', 0), ('', 0), ('',
 0), ('', 0), ('', 0), ('', 0)],
  dtype=[('f2', '|S20'), ('f0', 'i4')])

 ?

 In that case, that would be a great feature to add.
   
 That's what Travis is proposing. I would like to see a function that
 does this (however it is implemented under the covers):

   nra = subset_fields(ra, ['f0', 'f2'])
 

 Interesting.

   
 With the view, I don't think you can reorder the fields as in your
 example.
 

 That's a pity.  Providing a dtype with the notion of an internal reorder 
 can be very powerful in some situations.  But I guess that implementing 
 this would be complicated.

   
In general I agree with the idea but this starts sounding like R's data 
frames. So, is part of the goal to replicate some of the function of R's 
data frames?
For example the extract function 
(http://rweb.stat.umn.edu/R/library/base/html/Extract.data.frame.html)
(there is also the cookbook example of setting a name to null to remove 
it, see http://www.r-cookbook.com/node/50).



Bruce

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] passing a C array to embedded Python from C code

2008-10-30 Thread Christopher Barker

A few comments:

Chris LeBlanc wrote:
  I'm thinking I could create a
 NumPy array that uses the same C array (1d, 2d, or 3d) that our
 program is using instead of copying the memory.

yes, you can do that.

 At this point, I'm thinking the python code wouldn't need to return
 any objects because it would be modifying the seismic data (and trace
 headers) in memory.

There are a number of operations in Python that would require copies of 
the data, and you may want to be able to create data in python, and 
process it with your C code, so you probably do want to support passing 
data arrays from Python to the C code as well.


 Maybe calling a custom module from the Python code that does the C
 array to NumPy translation using Cython/pyrex/swig/etc.  Would it be
 possible to use the same C arrays from here without copying them?

yes, and you probably do want to use one of Cython, SWIG, or Ctypes -- 
it's really easier not to have to handle the reference counting, etc, 
yourself.

There are a set of SWIG typemaps and docs in the numpy distribution (in 
the docs dir, I believe), and there are various Wiki pages describing 
how  to do it with Cython and ctypes as well.

How to choose?

My thoughts on SWIG: SWIG is a pretty complex system, so you need to 
learn what is essentially yet another language. However, the numpy.i 
typemaps work well for the simple cases, so it may be easy to get 
started. The big advantage of SWIG is that it auto-generates the 
wrapper, once you have defined the typemaps you need. This gives a real 
advantage if you need to :

- provide wrappers for more than one language (PERL, Java, etc..)

- are wrapping a substantial library, particularly one that is under 
active development.

-Chris



-- 
Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/ORR(206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115   (206) 526-6317   main reception

[EMAIL PROTECTED]
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] passing a C array to embedded Python from C code

2008-10-30 Thread Christopher Barker
Matthieu Brucher wrote:
  If you can
 follow a French tutorial, you can go on
 http://matthieu-brucher.developpez.com/tutoriels/python/swig-numpy/#LV
 to have a skeletton for your issue.

That looks very useful -- any chance of an English translation? My one 
year of high school French is proving useless. Otherwise, the code 
itself is still quite helpful.

-Chris




-- 
Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/ORR(206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115   (206) 526-6317   main reception

[EMAIL PROTECTED]
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] passing a C array to embedded Python from C code

2008-10-30 Thread Matthieu Brucher
2008/10/30 Christopher Barker [EMAIL PROTECTED]:
 Matthieu Brucher wrote:
  If you can
 follow a French tutorial, you can go on
 http://matthieu-brucher.developpez.com/tutoriels/python/swig-numpy/#LV
 to have a skeletton for your issue.

 That looks very useful -- any chance of an English translation? My one
 year of high school French is proving useless. Otherwise, the code
 itself is still quite helpful.

 -Chris

I thought I put it on my blog, but no :|
I may find one day the time to translate it and to improve/enhance it.

Matthieu
-- 
Information System Engineer, Ph.D.
Website: http://matthieu-brucher.developpez.com/
Blogs: http://matt.eifelle.com and http://blog.developpez.com/?blog=92
LinkedIn: http://www.linkedin.com/in/matthieubrucher
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] passing a C array to embedded Python fromC code

2008-10-30 Thread Anthony Floyd
Hi Chris,

 Matthieu Brucher wrote:
   If you can
  follow a French tutorial, you can go on
  
 http://matthieu-brucher.developpez.com/tutoriels/python/swig-numpy/#LV
  to have a skeletton for your issue.
 
 That looks very useful -- any chance of an English 
 translation? My one 
 year of high school French is proving useless. Otherwise, the code 
 itself is still quite helpful.
 

Google does a pretty good job on this one:

http://translate.google.com/translate?u=http%3A%2F%2Fmatthieu-brucher.developpez.com%2Ftutoriels%2Fpython%2Fswig-numpy%2F%23LVsl=frtl=enhl=enie=UTF-8

Anthony.

--
Anthony Floyd, PhD
Convergent Manufacturing Technologies Inc.
6190 Agronomy Rd, Suite 403
Vancouver BC  V6T 1Z3
CANADA

Email: [EMAIL PROTECTED] | Tel:   604-822-9682 x102
WWW:   http://www.convergent.ca| Fax:   604-822-9659  

CMT is hiring: See http://www.convergent.ca for details
 
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Numpy matrix multiplication slow even though ATLAS linked

2008-10-30 Thread Charles R Harris
On Thu, Oct 30, 2008 at 5:19 AM, Jan-Willem van de Meent 
[EMAIL PROTECTED] wrote:

 Dear all,

 This is my first post to this list. I am having perfomance issues with with
 numpy/atlas. Doing dot(a,a) for a 2000x2000 matrix takes  about 1m40s, even
 though numpy is appears to link to my atlas libraries:

 I did a quick benchmark by running the following script:

#! /usr/bin/env python
import numpy
import time

try:
   import numpy.core._dotblas
   print 'Using ATLAS:'
except ImportError:
   print 'No ATLAS:'

t = time.time()
x = numpy.random.random((1000,1000))
y = numpy.random.random((1000,1000))
z = numpy.dot(x, y)

print time.time()-t

 My laptop is a Dell D620 Core Duo T2300 1.66 Ghz, running Archlinux with
 GCC
 4.3.2, atlas 3.8.2, python 2.5.2 and numpy 1.2.1. Output of the script
 above
 is:

Using ATLAS:
7.99549412727

 A department desktop PC, Pentium D 3.00 Ghz, running Scientific Linux, with
 GCC 4.1.2, atlas 3.7.30, python 2.5.1 and numpy 1.1.0, runs this test 24
 times faster:

Using ATLAS:
0.337520122528


About .40 here with numpy from svn.



 So even though _dotblas.so exists, matrix multiplication appears to run at
 pretty much the same speed as if atlas were not  available. Running ldd on
 _dotblas.so suggests  that numpy is indeed linking to the atlas libs:

 ldd /usr/lib/python2.5/site-packages/numpy/core/_dotblas.so
linux-gate.so.1 =  (0xb7fcf000)
libatlas.so = /usr/lib/libatlas.so (0xb7cb5000)
liblapack.so = /usr/lib/liblapack.so (0xb77ab000)
libcblas.so = /usr/lib/libcblas.so (0xb778b000)
libf77blas.so = /usr/lib/libf77blas.so (0xb776f000)
libpython2.5.so.1.0 = /usr/lib/libpython2.5.so.1.0 (0xb763)
libpthread.so.0 = /lib/libpthread.so.0 (0xb7618000)
libc.so.6 = /lib/libc.so.6 (0xb74d6000)
libm.so.6 = /lib/libm.so.6 (0xb74b)
libgfortran.so.3 = /usr/lib/libgfortran.so.3 (0xb73ff000)
libdl.so.2 = /lib/libdl.so.2 (0xb73fb000)
libutil.so.1 = /lib/libutil.so.1 (0xb73f6000)
/lib/ld-linux.so.2 (0xb7fd)


What's in /usr/local/lib? Do you have a 64 bit system? What does locate
libatlas return?

Chuck
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] Complete LAPACK needed

2008-10-30 Thread Frank Lagor
Dear all,

I need to use functions in scipy which depend on having a complete lapack
library.  However, I am having a bit of trouble installing numpy and
referencing a complete lapack library that I built.  I have a few questions
that I am hoping someone can help me answer:

1) This machine is a cluster, so do I need to be sure that the lapack
library is a .so file?  (I think only if I want the processors to access the
file at the same time, which is my need...) If so, how do I make a shared
library file?
2) If I can tell numpy to use some default version of lapack, will it be a
complete version?
3) I was unable to getting numpy installed by referencing the lapack.a
library from the site.cfg file. I also tried to set the LAPACK environment
variable in my .bashrc file, but this did not seem to work either.  The
problem is that when I configure numpy, it just does not recognize this file
as a lapack library.  It ignores it and tries to use an lapack lite library
located somewhere else (which doesn't even work anyway).

Any ideas? suggestions?
Thanks in advance!
Frank
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Numpy matrix multiplication slow even though ATLAS linked

2008-10-30 Thread Jan-Willem van de Meent
On Thursday 30 October 2008 18:41:51 Charles R Harris wrote:
 On Thu, Oct 30, 2008 at 5:19 AM, Jan-Willem van de Meent 

 [EMAIL PROTECTED] wrote:
  Dear all,
 
  This is my first post to this list. I am having perfomance issues with
  with numpy/atlas. Doing dot(a,a) for a 2000x2000 matrix takes  about
  1m40s, even though numpy is appears to link to my atlas libraries:
 
  I did a quick benchmark by running the following script:
 
 #! /usr/bin/env python
 import numpy
 import time
 
 try:
import numpy.core._dotblas
print 'Using ATLAS:'
 except ImportError:
print 'No ATLAS:'
 
 t = time.time()
 x = numpy.random.random((1000,1000))
 y = numpy.random.random((1000,1000))
 z = numpy.dot(x, y)
 
 print time.time()-t
 
  My laptop is a Dell D620 Core Duo T2300 1.66 Ghz, running Archlinux with
  GCC
  4.3.2, atlas 3.8.2, python 2.5.2 and numpy 1.2.1. Output of the script
  above
  is:
 
 Using ATLAS:
 7.99549412727
 
  A department desktop PC, Pentium D 3.00 Ghz, running Scientific Linux,
  with GCC 4.1.2, atlas 3.7.30, python 2.5.1 and numpy 1.1.0, runs this
  test 24 times faster:
 
 Using ATLAS:
 0.337520122528

 About .40 here with numpy from svn.

  So even though _dotblas.so exists, matrix multiplication appears to run
  at pretty much the same speed as if atlas were not  available. Running
  ldd on _dotblas.so suggests  that numpy is indeed linking to the atlas
  libs:
 
  ldd /usr/lib/python2.5/site-packages/numpy/core/_dotblas.so
 linux-gate.so.1 =  (0xb7fcf000)
 libatlas.so = /usr/lib/libatlas.so (0xb7cb5000)
 liblapack.so = /usr/lib/liblapack.so (0xb77ab000)
 libcblas.so = /usr/lib/libcblas.so (0xb778b000)
 libf77blas.so = /usr/lib/libf77blas.so (0xb776f000)
 libpython2.5.so.1.0 = /usr/lib/libpython2.5.so.1.0 (0xb763)
 libpthread.so.0 = /lib/libpthread.so.0 (0xb7618000)
 libc.so.6 = /lib/libc.so.6 (0xb74d6000)
 libm.so.6 = /lib/libm.so.6 (0xb74b)
 libgfortran.so.3 = /usr/lib/libgfortran.so.3 (0xb73ff000)
 libdl.so.2 = /lib/libdl.so.2 (0xb73fb000)
 libutil.so.1 = /lib/libutil.so.1 (0xb73f6000)
 /lib/ld-linux.so.2 (0xb7fd)

 What's in /usr/local/lib? Do you have a 64 bit system? What does locate
 libatlas return?

 Chuck

Thanks for the response. I'm on 32 bit and all ATLAS files are installed 
in /usr/lib (headers are in /usr/include/atlas)

/usr/lib/libatlas.a
/usr/lib/libatlas.so
/usr/lib/libblas.so
/usr/lib/libblas.so.3
/usr/lib/libblas.so.3.0.3
/usr/lib/libcblas.a
/usr/lib/libcblas.so
/usr/lib/libf77blas.a
/usr/lib/libf77blas.so
/usr/lib/liblapack.a
/usr/lib/liblapack.so
/usr/lib/liblapack.so.3
/usr/lib/libptcblas.a
/usr/lib/libptf77blas.a

I posted a full build log earlier today, but I think it is still awaiting 
approval of moderators because of file size limitations.

JW
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Numpy matrix multiplication slow even though ATLAS linked

2008-10-30 Thread Robert Kern
On Thu, Oct 30, 2008 at 17:19, Jan-Willem van de Meent
[EMAIL PROTECTED] wrote:
 I posted a full build log earlier today, but I think it is still awaiting
 approval of moderators because of file size limitations.

There are no moderators to approve it. The file size limitation is
just a limitation. If you have web space to place the log, you can
give us a URL.

-- 
Robert Kern

I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth.
  -- Umberto Eco
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] Large array using 4 times as much memory as it should

2008-10-30 Thread Rick Giuly
Hello All,

I find that python is using about four times as much memory as it should 
need for arrays. This is problematic as I need to use all available 
memory for large 3D imaging datasets. Is there a way to get around this 
problem? Am I making a mistake? Is it a bug?

(I'm running windowsXP 32bit with 760M of memory Available according 
to the Performance pane of the task manager.)

versions: numpy 1.2.0 with python 2.5.2

Any help is appreciated


-Rick




**
Details of my testing:

Each test was run from the command line and for each test python was 
restarted.

Testing a 50M array:
a = numpy.ones((1024,1024,50), dtype=numpy.uint32)
The available memory dropped by 200M


Testing a 100M array:
a = numpy.ones((1024,1024,100), dtype=numpy.uint32)
The available memory dropped by 400M


Testing a 200M array:
a = numpy.ones((1024,1024,200), dtype=numpy.uint32)
The available memory dropped by 750M


Testing a 300M array:
a = numpy.ones((1024,1024,300), dtype=numpy.uint32)
an error occurs:
Traceback (most recent call last):
   File stdin, line 1, in module
   File 
o:\software\pythonxy\python\lib\site-packages\numpy\core\numeric.py, li
ne 1445, in ones
 a = empty(shape, dtype, order)
MemoryError







___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Large array using 4 times as much memory as it should

2008-10-30 Thread Robert Kern
On Thu, Oct 30, 2008 at 20:41, Rick Giuly [EMAIL PROTECTED] wrote:
 Hello All,

 I find that python is using about four times as much memory as it should
 need for arrays. This is problematic as I need to use all available
 memory for large 3D imaging datasets. Is there a way to get around this
 problem? Am I making a mistake? Is it a bug?

You are making a mistake. uint32s take up 4 bytes each. So
(1024,1024,50) takes up 4*50*1024*1024 bytes == 200 Mb.

-- 
Robert Kern

I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth.
  -- Umberto Eco
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] [Fwd: Large array using 4 times as much memory as it should]

2008-10-30 Thread Rick Giuly
Please disregard the last message.

I just realised 8*4=32, very sorry about this. Apparently I need some sleep.

-Rick

 Original Message 
Subject: Large array using 4 times as much memory as it should
Date: Thu, 30 Oct 2008 18:41:44 -0700
From: Rick Giuly [EMAIL PROTECTED]
To: numpy-discussion@scipy.org

Hello All,

I find that python is using about four times as much memory as it should
need for arrays. This is problematic as I need to use all available
memory for large 3D imaging datasets. Is there a way to get around this
problem? Am I making a mistake? Is it a bug?

(I'm running windowsXP 32bit with 760M of memory Available according
to the Performance pane of the task manager.)

versions: numpy 1.2.0 with python 2.5.2

Any help is appreciated


-Rick




**
Details of my testing:

Each test was run from the command line and for each test python was
restarted.

Testing a 50M array:
a = numpy.ones((1024,1024,50), dtype=numpy.uint32)
The available memory dropped by 200M


Testing a 100M array:
a = numpy.ones((1024,1024,100), dtype=numpy.uint32)
The available memory dropped by 400M


Testing a 200M array:
a = numpy.ones((1024,1024,200), dtype=numpy.uint32)
The available memory dropped by 750M


Testing a 300M array:
a = numpy.ones((1024,1024,300), dtype=numpy.uint32)
an error occurs:
Traceback (most recent call last):
   File stdin, line 1, in module
   File
o:\software\pythonxy\python\lib\site-packages\numpy\core\numeric.py, li
ne 1445, in ones
 a = empty(shape, dtype, order)
MemoryError







___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] [Fwd: Large array using 4 times as much memory as it should]

2008-10-30 Thread Fernando Perez
On Thu, Oct 30, 2008 at 6:48 PM, Rick Giuly [EMAIL PROTECTED] wrote:
 Please disregard the last message.

 I just realised 8*4=32, very sorry about this. Apparently I need some sleep.

 -Rick

  Original Message 
 Subject: Large array using 4 times as much memory as it should
 Date: Thu, 30 Oct 2008 18:41:44 -0700
 From: Rick Giuly [EMAIL PROTECTED]
 To: numpy-discussion@scipy.org

 Hello All,

 I find that python is using about four times as much memory as it should

You may find this handy for looking at your variable's memory use:

In [12]: a = numpy.ones((1024,1024,50), dtype=numpy.uint32)

In [13]: print 'approx size in Mb:',(a.size*a.itemsize)/1e6
approx size in Mb: 209.7152

numpy.who() gives you a nice summary:

In [14]: numpy.who()
NameShape  BytesType
=

a   1024 x 1024 x 50   209715200uint32

Upper bound on total bytes  =   209715200


Cheers,

f
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] unable to build atlas and lapack correctly

2008-10-30 Thread Alexandre Lacoste
Hi I'm trying to install numpy on a x86_64 ubuntu hardy and I'm having a
hard time to get atlas to be linked correctly.

Installing it through apt-get install with the astraw repository worked on
my intel dual core. But now I need to install it in on a x86_64 and the
astraw binaries will link against libatlas3gf-base instead of
libatlas3gf-3dnow because the 3dnow version simply doesn't exist for x86_64
arch. So I still get slow matrix multiplications ...

Then I went through the instructions for building atlas by hand. Then I
build numpy. Everything went ok... Now when comes the time to build scipy, I
get the following error :

from numpy.linalg import lapack_lite
ImportError: /usr/local/atlas/lib/liblapack.so: undefined symbol: ztbsv_

$ uname -a
Linux graal1 2.6.24-21-generic #1 SMP Mon Aug 25 16:57:51 UTC 2008 x86_64
GNU/Linux

$ ldd liblapack.so
linux-vdso.so.1 =  (0x7fffe29fe000)
libgfortran.so.2 = /usr/lib/libgfortran.so.2 (0x7febd9ee3000)
libm.so.6 = /lib/libm.so.6 (0x7febd9c62000)
libc.so.6 = /lib/libc.so.6 (0x7febd98ff000)
/lib64/ld-linux-x86-64.so.2 (0x7febda93b000)


I'm getting clueless and googling for the error message return noting useful
any thoughts ?
thanks...

ªŁ€×
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] unable to build atlas and lapack correctly

2008-10-30 Thread David Cournapeau
Alexandre Lacoste wrote:
 Hi I'm trying to install numpy on a x86_64 ubuntu hardy and I'm having
 a hard time to get atlas to be linked correctly.

 Installing it through apt-get install with the astraw repository
 worked on my intel dual core. But now I need to install it in on a
 x86_64 and the astraw binaries will link against libatlas3gf-base
 instead of libatlas3gf-3dnow because the 3dnow version simply doesn't
 exist for x86_64 arch. So I still get slow matrix multiplications ...

You can simply use the sse2 version of atlas: SSE2 is available on any
x86_64 CPU.

Note that on Ubuntu, you can link to libatals3fg-base (the gfortran,
basic version of ATLAS), and after building numpy, installing atlas for
SSE2 and still get the speed up, thanks to the hwcap capability of the
GNU loader. The hwcap means that when the loaded code (numpy) request a
library (atlas), it will first look for a specially optimized one in
some directory (/usr/lib/sse2, for example), and will use the default
one if the optimized one is not there.



 Then I went through the instructions for building atlas by hand. Then
 I build numpy. Everything went ok... Now when comes the time to build
 scipy, I get the following error :

 from numpy.linalg import lapack_lite
 ImportError: /usr/local/atlas/lib/
 liblapack.so: undefined symbol: ztbsv_

You made a mistake when building lapack and/or atlas. Those are
difficult to build correctly, you should avoid build them by yourself as
much as possible,

cheers,

David
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] (no subject)

2008-10-30 Thread frank wang

Hi,
 
In my work, I want to implement a fir filter with an input array. Since 
performing the filter on each input sample is slow, are there fast way to 
perform the fir filter operation?  Are there ways to convert input into an 
array and perform the array multipication?
 
Thanks
 
Frank
_
When your life is on the go—take your life with you.
http://clk.atdmt.com/MRT/go/115298558/direct/01/___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] (no subject)

2008-10-30 Thread David Cournapeau
frank wang wrote:
 Hi,
  
 In my work, I want to implement a fir filter with an input array.
 Since performing the filter on each input sample is slow, are there
 fast way to perform the fir filter operation?  Are there ways to
 convert input into an array and perform the array multipication?

look at scipy.signal.lfilter,

David
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] unable to build atlas and lapack correctly

2008-10-30 Thread Alexandre Lacoste
Thanks for your answer


 You can simply use the sse2 version of atlas: SSE2 is available on any
 x86_64 CPU.


Are you sure ? I can't find it...
http://packages.ubuntu.com/search?arch=amd64searchon=nameskeywords=atlas%20sse2



 Note that on Ubuntu, you can link to libatals3fg-base (the gfortran,
 basic version of ATLAS), and after building numpy, installing atlas for
 SSE2 and still get the speed up, thanks to the hwcap capability of the
 GNU loader. The hwcap means that when the loaded code (numpy) request a
 library (atlas), it will first look for a specially optimized one in
 some directory (/usr/lib/sse2, for example), and will use the default
 one if the optimized one is not there.


That's what I tried in the first place but I must have done something wrong
on that part too ;)



 You made a mistake when building lapack and/or atlas. Those are
 difficult to build correctly, you should avoid build them by yourself as
 much as possible,


I did it a couple of time in the past ... even on 64 bits machines... I
guess I was lucky back then.
I really need to get it work. I'll give it a couple of other tries. It is
only annoying to have to wait half an hour each time. If you can give me a
couple of other hints, maybe I'll target the problem faster.

Thanks again


ªŁ€×
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] help to speed up the python code

2008-10-30 Thread frank wang

Hi, Bob,
 
The problem is that I want to resample my data with another sampling rate. the 
two rates is very close. I use the formula:
 
s(t)=sum(a_k*sinc(t-kTs)).
 
the new sampling rate is Ts', so I have 
s(nTs')=sum(a_k*sinc(nTs'-kTs)). The sum index k is over the (-P, P), Centered 
at n. The n is start from zero. THe code is using two for loops and it is slow. 
The length of s(nTs) is very long, so it takes quite long time to do it. 
 
Thanks
 
Frank Date: Sun, 26 Oct 2008 22:45:42 -0500 From: [EMAIL PROTECTED] To: 
numpy-discussion@scipy.org Subject: Re: [Numpy-discussion] help to speed up 
the python code  On Fri, Oct 24, 2008 at 11:30, frank wang [EMAIL 
PROTECTED] wrote:  Hi,   I have to send this request second time since 
my first message contains the  attached data file which is too big and was 
blocked by the system. So this  time I will not attach the data file.   I 
have converted a matlab function to python using numpy. both matlab and  
python run slow. I know that numpy has a lot of features, so I hope some  
experts can help me to speed up the code.  Can you describe in higher level 
terms what you are trying to do? I'm having trouble following the code.  -- 
 Robert Kern  I have come to believe that the whole world is an enigma, a 
harmless enigma that is made terrible by our own mad attempt to interpret it 
as though it had an underlying truth. -- Umberto Eco 
___ Numpy-discussion mailing list 
Numpy-discussion@scipy.org 
http://projects.scipy.org/mailman/listinfo/numpy-discussion
_
Stay organized with simple drag and drop from Windows Live Hotmail.
http://windowslive.com/Explore/hotmail?ocid=TXT_TAGLM_WL_hotmail_102008___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] help to speed up the python code

2008-10-30 Thread David Cournapeau
frank wang wrote:
 Hi, Bob,
  
 The problem is that I want to resample my data with another sampling
 rate. the two rates is very close. I use the formula:
  
 s(t)=sum(a_k*sinc(t-kTs)).
  
 the new sampling rate is Ts', so I have
 s(nTs')=sum(a_k*sinc(nTs'-kTs)). The sum index k is over the (-P, P),
 Centered at n. The n is start from zero. THe code is using two for
 loops and it is slow. The length of s(nTs) is very long, so it takes
 quite long time to do it.

If you want to do high quality resampling, you may want to look at the
samplerate scikits. It uses sinc interpolation for resampling (the
scikits is just a wrapper around SRC: http://www.mega-nerd.com/SRC/, if
you want more details on the implementation)

The documentation is really light unfortunately, but there is an example:

http://www.ar.media.kyoto-u.ac.jp/members/david/softwares/samplerate/index.html

The implementation of SRC is in C, and you can choose between 3 modes
for quality (plus other interpolations like linear, but those are ill
suited for anything serious if you can afford the computational cost),

cheers,

David
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion