Re: [Numpy-discussion] Backout complex changes for next release?

2009-11-07 Thread David Cournapeau
On Sat, Nov 7, 2009 at 7:12 AM, Charles R Harris
charlesr.har...@gmail.com wrote:
 There seems to be a rat's nest of problems showing up in scipy due to the
 recent changes in numpy complex support. The problems are of two basic
 sorts: 1) reserved name conflicts and 2) file conflicts. The first could be
 dealt with, but the second might be trickier, both c/c++ defined complex.h
 but the files are different. C++ has moved on, but not all platforms seem to
 have followed.

The solution for 2 is easy: I should not have included complex.h in
the public headers. I don't think you can include complex.h in C++
without causing trouble.

The scipy C++ code will need to be fixed to avoid using reserved
names, but that should not be difficult,

David
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] dtype 'groups'

2009-11-07 Thread Stéfan van der Walt
2009/11/7 David Warde-Farley d...@cs.toronto.edu:
 Just to confirm my suspicions, what's the preferred way to check if a
 dtype is...

Not sure what the official way is, but this method works well:

 a) an integer type vs. a floating point type vs. a string type? I'm
 assuming dt.kind.

np.issubdtype(dt, float)

Regards
Stéfan
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Vector interpolation on a 2D grid

2009-11-07 Thread Christopher Barker


Pierre GM wrote:
 I have a vector of observations (latitude,longitude,value) that I need  
 to interpolate over a given area.

 Some colleagues complained that the result looked a bit too choppy,  
 meaning that too much weight was given to the actual observations.
 What are my other options ?

Kriging might work well for this. I don't know of any python accessible 
libs, but I'd like to hear about it if anyone else does.

-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

chris.bar...@noaa.gov
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] ctypes and numpy

2009-11-07 Thread Zachary Pincus
Check out this thread:
http://www.mail-archive.com/numpy-discuss...@lists.sourceforge.net/msg01154.html

In shot, it can be done, but it can be tricky to make sure you don't  
leak memory. A better option if possible is to pre-allocate the array  
with numpy and pass that buffer into the C code -- but that's not  
always possible...

Zach


On Nov 6, 2009, at 10:18 PM, Chris Colbert wrote:

 Are you bound to using ctypes as the only option?

 This could be done quite easily in Cython...

 On Fri, Nov 6, 2009 at 1:57 PM, Trevor Clarke tre...@notcows.com  
 wrote:
 I've found some info on ctypes and numpy but it mostly tells me how  
 to pass
 numpy objects to C funcs. I've got a C func which returns a  
 c_void_p to a
 buffer and I'd like to access the data using numpy without creating  
 a copy
 of the data. Could someone point me in the right direction?

 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion


 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Deprecate np.max/np.min ?

2009-11-07 Thread Neil Crighton
Charles R Harris charlesr.harris at gmail.com writes:

 People import these functions -- yes, they shouldn't do that -- and the python
builtin versions are overloaded, causing hard to locate errors.

While I would love less duplication in the numpy namespace, I don't think the 
small gain here is worth the pain of deprecation. 

 OTOH, one can ask, why is
 
   np.min(3, 2)
 
 allowed when
 
   np.min([3], 2)
 
 gives ValueError: axis(=2) out of bounds. It seems to me that
 0-dimensional objects should accept only None as the axis? (Fixing this
 would also make misuse of np.min and np.max more difficult.)

I think it would be better to fix this issue.  np.min(3,2) should also give
ValueError: axis(=2) out of bounds. Fixing this also removes any possibility
of generating hard-to-find errors by overwriting the builtin min/max. (Unless
there's some corner case I'm missing).

Neil

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] How to get rid of the loop?

2009-11-07 Thread Stas K
Can I get rid of the loop in this example? And what is the fastest way  
to get v in the example?

ar = array([1,2,3])
for a in ar:
for b in ar:
v = a**2+b**2
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] How to get rid of the loop?

2009-11-07 Thread josef . pktd
On Sat, Nov 7, 2009 at 1:51 PM, Stas K stanc...@gmail.com wrote:
 Can I get rid of the loop in this example? And what is the fastest way
 to get v in the example?

 ar = array([1,2,3])
 for a in ar:
    for b in ar:
        v = a**2+b**2

 ar[:,None]**2 + ar**2
array([[ 2,  5, 10],
   [ 5,  8, 13],
   [10, 13, 18]])

I think, for this case there is also directly a function in numpy
hypot which should also work.

Josef
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] How to get rid of the loop?

2009-11-07 Thread Alan G Isaac
On 11/7/2009 1:51 PM, Stas K wrote:
 Can I get rid of the loop in this example? And what is the fastest way
 to get v in the example?

 ar = array([1,2,3])
 for a in ar:
  for b in ar:
  v = a**2+b**2


 a2 = a*a
 np.add.outer(a2,a2)
array([[ 2,  5, 10],
[ 5,  8, 13],
[10, 13, 18]])

hth,
Alan Isaac

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] How to get rid of the loop?

2009-11-07 Thread Stas K

Thank you, Josef

It is exactly what I want:

ar[:,None]**2 + ar**2


Do you know something about performance of this? In my real program  
ar  have ~ 10k elements, and expression for v more complicated (it has  
some trigonometric functions)



On 07.11.2009, at 21:57, josef.p...@gmail.com wrote:


On Sat, Nov 7, 2009 at 1:51 PM, Stas K stanc...@gmail.com wrote:
Can I get rid of the loop in this example? And what is the fastest  
way

to get v in the example?

ar = array([1,2,3])
for a in ar:
   for b in ar:
   v = a**2+b**2



ar[:,None]**2 + ar**2

array([[ 2,  5, 10],
  [ 5,  8, 13],
  [10, 13, 18]])

I think, for this case there is also directly a function in numpy
hypot which should also work.

Josef

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] PGI vs. Intel Compiler

2009-11-07 Thread Sebastian Żurek
Hi!

I'm not sure, if it's the right group to ask, but it's a kind of a numeric 
question involving Python, so...

My faculty is going to spend some money on commercial compiler. The
two choices are:
- PGI C++ compiler 
- Intel C++ compiler

I wonder, if anyone as some experience with Python numerical packages  (and 
Python itself) compiled with the compiler listed above? Can anyone provide 
some pro and cons?

How about the basic tools like BLAS/LAPACK behavior compiled with the 
PGI/Intel C++ tools?

Any suggestion is welcome!

Thanks in advance!
Sebastian


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] masked record arrays

2009-11-07 Thread Thomas Robitaille


Pierre GM-2 wrote:
 
 Mmh. With a recent (1.3) version of numpy, you should already be able  
 to mask individual fields of a structured array without problems. If  
 you need fields to be accessed as attributes the np.recarray way, you  
 can give numpy.ma.mrecords.MaskedRecords a try. It's been a while I  
 haven't touched it, so you may run into the occasional bug. FYI, I'm  
 not a big fan of record arrays and tend to prefer structured ones...
 What two implementations were you talking about ?
 In any case, feel free to try and please, report any issue you run  
 into with MaskedRecords.
 Cheers
 

Thanks for the advice! I'm somewhat confused by the difference between
structured and record arrays. My understanding is that record arrays allow
you to access fields by attribute (e.g. r.field_name), but I imagine that
there are much more fundamental differences for the two to be treated
separately in numpy. I find the numpy documentation somewhat confusing in
that respect - if you have a look at this page

http://docs.scipy.org/doc/numpy/user/basics.rec.html

I think the 'aka record arrays' is especially confusing as this would
suggest the two are the same. So is there good information anywhere about
what exactly are the differences between the two? This page is also
confusing:

http://docs.scipy.org/doc/numpy/reference/generated/numpy.recarray.html

as to me Construct an ndarray that allows field access using attributes
suggests that all a recarray is is an ndarray/structured array with
overloaded __getattr__/__setattr__ methods. Is that all recarrays are? If
so, why was a completely separate package developed for masked record arrays
- can one not just use masked structured arrays and overload
getattr/setattr?

Cheers,

Thomas

-- 
View this message in context: 
http://old.nabble.com/masked-record-arrays-tp26237612p26247808.html
Sent from the Numpy-discussion mailing list archive at Nabble.com.

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] How to get rid of the loop?

2009-11-07 Thread Anne Archibald
2009/11/7 Stas K stanc...@gmail.com:
 Thank you, Josef
 It is exactly what I want:

 ar[:,None]**2 + ar**2

 Do you know something about performance of this? In my real program ar  have
 ~ 10k elements, and expression for v more complicated (it has some
 trigonometric functions)

The construction of ar[:,None] (which I prefer to spell
ar[:,np.newaxis]) is cheap, since it just constructs a view into ar.
Computing ar**2 and ar[:,None] require squaring each element of ar,
but this is done in a C loop. (It's done twice in this expression, so
this isn't as efficient as it might be). Only when it comes time to do
the addition do you do n**2 operations. For very large n, this can be
a crushing memory burden, and if you only need these values for an
intermediate calculation it sometimes turns out that it's better to
loop over one of the dimensions. But this expression is probably about
as efficient as you can hope for, given what it does. If you need to
do some more complicated calculation, the main thing to be careful of
is that you do as much calculation as possible on the separate arrays,
while they're still only n elements and not n**2.


Anne

 On 07.11.2009, at 21:57, josef.p...@gmail.com wrote:

 On Sat, Nov 7, 2009 at 1:51 PM, Stas K stanc...@gmail.com wrote:

 Can I get rid of the loop in this example? And what is the fastest way

 to get v in the example?

 ar = array([1,2,3])

 for a in ar:

    for b in ar:

        v = a**2+b**2

 ar[:,None]**2 + ar**2

 array([[ 2,  5, 10],
   [ 5,  8, 13],
   [10, 13, 18]])

 I think, for this case there is also directly a function in numpy
 hypot which should also work.

 Josef

 ___

 NumPy-Discussion mailing list

 NumPy-Discussion@scipy.org

 http://mail.scipy.org/mailman/listinfo/numpy-discussion

 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion


 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Random int64 and float64 numbers

2009-11-07 Thread Sturla Molden
David Cournapeau wrote:
 On Fri, Nov 6, 2009 at 6:54 AM, David Goldsmith d.l.goldsm...@gmail.com 
 wrote:
   
 Interesting thread, which leaves me wondering two things: is it documented
 somewhere (e.g., at the IEEE site) precisely how many *decimal* mantissae
 are representable using the 64-bit IEEE standard for float representation
 (if that makes sense); and are such decimal mantissae uniformly distributed
 

 They are definitely not uniformly distributed: that's why two numbers
 are close around 1 when they have only a few EPS difference, but
 around 1e100, you have to add quite a few EPS to even get a different
 number at all.

 That may be my audio processing background, but I like to think about
 float as numbers which have the same relative precision at any level
 - a kind of dB scale. If you want numbers with a fixed number of
 decimals, you need a fixed point representation.
   

David Godsmith was asking about the mantissae. For a double, that is a 
53 bit signed integer. I.e. you have 52 bit fractional part (bit 0-51), 
11 bit exponent (bit 52-62), and one sign bit (bit 63). The mantissae is 
uniformly distributed like any signed integer. The mantissae of a double 
have 2**53 different integer values: -2**52 to 2**52-1.

But the value of a floating point number is

   value = (-1)**signbit * 2**(exponent - bias) * (1 - fraction)

with bias = 1023 for a double. Thus, floating point numbers are not 
uniformly distributed, but the mantissae is.

For numerical illiterates this might come as a surprise. But in 
numerical mathematics, the resolution is in the number of significant 
digits, not in the number of decimals. 101 and .00201 have the same 
numerical precision.

A decimal, on the other hand, can be thought of as a floating point 
number using base-10 instead of base-2 for the exponent:

   value = (-1)**signbit * 10**(exponent - bias) * (1 - fraction)

Decimals and floats are not fundamentally different. There are number 
exactly representable with a decimal that cannot be exactly represented 
with a float. But numerical computation do not become more precise with 
a decimal than a float.

















___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] Accessing LAPACK and BLAS from the numpy C API

2009-11-07 Thread Jake VanderPlas
Hello,
I'm working on wrapping a set of C++ routines for manifold learning
(LLE, Isomap, LTSA, etc) in python.  In the LLE routine, it is
necessary to loop through the input points and perform an svd of each
local covariance matrix.  Presently I have my own C-LAPACK wrapper
that I call within a C loop, but this seems non-ideal because numpy
already wraps LAPACK/ATLAS for linear algebra.  Does anybody know a
way to directly access the numpy.linalg routines from a C extension,
without the overhead of a python callback?  Thanks for the help.
   -Jake
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Accessing LAPACK and BLAS from the numpy C API

2009-11-07 Thread Sturla Molden
Jake VanderPlas wrote:
 Does anybody know a
 way to directly access the numpy.linalg routines from a C extension,
 without the overhead of a python callback?  Thanks for the help.
   

You find a C function pointer wrapped in a CObject in the ._cpointer 
attribute.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] Use-case for np.choose

2009-11-07 Thread David Goldsmith
Hi, all!  I'm working to clarify the docstring for np.choose (Stefan pointed
out to me that it is pretty unclear, and I agreed), and, now that (I'm
pretty sure that) I've figured out what it does in its full generality
(e.g., when the 'choices' array is greater than 2-D), I'm curious as to
use-cases.  (Not that I'm suggesting that we deprecate it, but I am curious
as to whether or not anyone is presently using it and if so, how they're
using it, esp. if anyone *is* using it with 'choices' arrays 3-D or
greater.)

Also, my experimenting suggests that the index array ('a', the first
argument in the func. sig.) *must* have shape (choices.shape[-1],) - someone
please let me know ASAP if this is not the case, and please furnish me w/ a
counterexample because I was unable to generate one myself.

Thanks,

DG
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Use-case for np.choose

2009-11-07 Thread Anne Archibald
2009/11/7 David Goldsmith d.l.goldsm...@gmail.com:
 Hi, all!  I'm working to clarify the docstring for np.choose (Stefan pointed
 out to me that it is pretty unclear, and I agreed), and, now that (I'm
 pretty sure that) I've figured out what it does in its full generality
 (e.g., when the 'choices' array is greater than 2-D), I'm curious as to
 use-cases.  (Not that I'm suggesting that we deprecate it, but I am curious
 as to whether or not anyone is presently using it and if so, how they're
 using it, esp. if anyone *is* using it with 'choices' arrays 3-D or
 greater.)

It's a generalized np.where(), allowing more than two options:

In [10]: a = np.arange(2*3*5).reshape(2,3,5) % 3

In [11]: o = np.ones((2,3,5))

In [12]: np.choose(a,(o,0.5*o,0.1*o))
Out[12]:
array([[[ 1. ,  0.5,  0.1,  1. ,  0.5],
[ 0.1,  1. ,  0.5,  0.1,  1. ],
[ 0.5,  0.1,  1. ,  0.5,  0.1]],

   [[ 1. ,  0.5,  0.1,  1. ,  0.5],
[ 0.1,  1. ,  0.5,  0.1,  1. ],
[ 0.5,  0.1,  1. ,  0.5,  0.1]]])

 Also, my experimenting suggests that the index array ('a', the first
 argument in the func. sig.) *must* have shape (choices.shape[-1],) - someone
 please let me know ASAP if this is not the case, and please furnish me w/ a
 counterexample because I was unable to generate one myself.

It seems like a and each of the choices must have the same shape (with
the exception that choices acn be scalars), but I would consider this
a bug. Really, a and all the choices should be broadcast to the same
shape. Or maybe it doesn't make sense to broadcast a - it could be
valuable to know that the result is always exactly the same shape as a
- but broadcasting all the choice arrays presents an important
improvement of choose over fancy indexing. There's a reason choose
accepts a sequence of arrays as its second argument, rather than a
higher-dimensional array.

Anne

 Thanks,

 DG

 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Deprecate np.max/np.min ?

2009-11-07 Thread Pauli Virtanen
la, 2009-11-07 kello 18:27 +, Neil Crighton kirjoitti:
[clip]
 I think it would be better to fix this issue.  np.min(3,2) should also give
 ValueError: axis(=2) out of bounds. Fixing this also removes any possibility
 of generating hard-to-find errors by overwriting the builtin min/max. (Unless
 there's some corner case I'm missing).

Fixed in r7697, r7698, which change the behavior of scalars in many
functions using axis=xxx.

I don't think this breaks anyone's code -- using out-of-bounds values of
axis is almost certainly an error.

I left axis=-1 and axis=0 allowed, in addition to axis=None. These
seemed to be required by at least the masked arrays unit tests...

-- 
Pauli Virtanen



___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] ctypes and numpy

2009-11-07 Thread Trevor Clarke
this looks like what I need...I'm not concerned with leaking memory as it's
a borrowed pointer which will be cleaned up in C code later. Thanks for the
pointer.

On Sat, Nov 7, 2009 at 12:36 PM, Zachary Pincus zachary.pin...@yale.eduwrote:

 Check out this thread:

 http://www.mail-archive.com/numpy-discuss...@lists.sourceforge.net/msg01154.html

 In shot, it can be done, but it can be tricky to make sure you don't
 leak memory. A better option if possible is to pre-allocate the array
 with numpy and pass that buffer into the C code -- but that's not
 always possible...

 Zach


 On Nov 6, 2009, at 10:18 PM, Chris Colbert wrote:

  Are you bound to using ctypes as the only option?
 
  This could be done quite easily in Cython...
 
  On Fri, Nov 6, 2009 at 1:57 PM, Trevor Clarke tre...@notcows.com
  wrote:
  I've found some info on ctypes and numpy but it mostly tells me how
  to pass
  numpy objects to C funcs. I've got a C func which returns a
  c_void_p to a
  buffer and I'd like to access the data using numpy without creating
  a copy
  of the data. Could someone point me in the right direction?
 
  ___
  NumPy-Discussion mailing list
  NumPy-Discussion@scipy.org
  http://mail.scipy.org/mailman/listinfo/numpy-discussion
 
 
  ___
  NumPy-Discussion mailing list
  NumPy-Discussion@scipy.org
  http://mail.scipy.org/mailman/listinfo/numpy-discussion

 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] masked record arrays

2009-11-07 Thread Pierre GM

On Nov 7, 2009, at 2:26 PM, Thomas Robitaille wrote:

 Thanks for the advice! I'm somewhat confused by the difference between
 structured and record arrays. My understanding is that record arrays  
 allow
 you to access fields by attribute (e.g. r.field_name), but I imagine  
 that
 there are much more fundamental differences for the two to be treated
 separately in numpy.

Actually, no. recarray is just ndarray w/ a special __getattribute__/ 
__setattr__ . They bring the convenience of exposing fields as  
properties, but they come to the cost of overloading __getattribute__

 I find the numpy documentation somewhat confusing in
 that respect - if you have a look at this page

 http://docs.scipy.org/doc/numpy/user/basics.rec.html

 I think the 'aka record arrays' is especially confusing as this would
 suggest the two are the same.

Not the most fortunate formulation, true...


 So is there good information anywhere about
 what exactly are the differences between the two? This page is also
 confusing:

 http://docs.scipy.org/doc/numpy/reference/generated/ 
 numpy.recarray.html

 as to me Construct an ndarray that allows field access using  
 attributes
 suggests that all a recarray is is an ndarray/structured array with
 overloaded __getattr__/__setattr__ methods. Is that all recarrays are?

Yep.

 If
 so, why was a completely separate package developed for masked  
 record arrays
 - can one not just use masked structured arrays and overload
 getattr/setattr?

Mostly historical reasons. Initially, there was only limited support  
for structured masked arrays and masked records filled the gap (albeit  
experimentally). With the 1.3 release, MaskedArray fully supports  
structured type, giving the possibility to mask individual fields and  
masked records became less useful. Still, it's cleaner to have a  
specific module where to  store functions like fromrecords, fromarrays  
and so forth. Note that the doc of numpy.ma.mrecords is a tad  
outdated, any help to this regard would be welcome.


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Vector interpolation on a 2D grid (with realistic results)

2009-11-07 Thread Pierre GM

On Nov 6, 2009, at 5:45 PM, Pauli Virtanen wrote:

 pe, 2009-11-06 kello 17:20 -0500, Pierre GM kirjoitti:
 All,
 I have a vector of observations (latitude,longitude,value) that I
 need
 to interpolate over a given area.


 You could try to use linear interpolation from the delaynay package
 (only supports slicing, so is a bit tricky to use). I may be slightly
 smoother than the natural neighbour one.

 For scattered data, when I wasn't happy with delaunay, I've also used
 the radial basis functions (Rbf) from scipy.interpolate -- something
 like

interpolator = Rbf(x, y, z, function='thin-plate', smooth=1e-9)

 should give an interpolator that is reasonably smooth between the data
 points. (I'd guess you can give it leeway in trading interpolation to
 smoothness by increasing the smooth parameter.) Other Rbf functions  
 than
 thin-plate splines might also work.

Linear interpolation with the delaunay package doesn't work great for  
my data. I played with the radial basis functions, but I'm afraid  
they're leading me down the dark, dark path of parameter fiddling. In  
particular, I'm not sure how to prevent my interpolated values to be  
bounded by the min and max of my actual observations.
Ralf' suggestion of smoothing the values afterwards is tempting, but  
sounds a bit too ad-hoc (btw, Ralf, those were relative differences of  
monthly average precipitation between El Niño and Neutral phases for  
November).

Chris, I gonna poke around and try to find some kriging algorithms.  
I'll report in a few. In the meantime, if anybody has anythng already  
implemented, please just let us know.


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Vector interpolation on a 2D grid (with realistic results)

2009-11-07 Thread Ryan May
On Sat, Nov 7, 2009 at 5:38 PM, Pierre GM pgmdevl...@gmail.com wrote:
 Linear interpolation with the delaunay package doesn't work great for
 my data. I played with the radial basis functions, but I'm afraid
 they're leading me down the dark, dark path of parameter fiddling. In
 particular, I'm not sure how to prevent my interpolated values to be
 bounded by the min and max of my actual observations.
 Ralf' suggestion of smoothing the values afterwards is tempting, but
 sounds a bit too ad-hoc (btw, Ralf, those were relative differences of
 monthly average precipitation between El Niño and Neutral phases for
 November).

That was me, not Ralf. :)  And I agree, I the interpolated field does
look a bit noisy for such data.  I've been doing the smoothing on top
of natural neighbor for doing some of my own meteorological analysis.
Using the Gaussian kernel isn't really *that* ad hoc considering the
prevalence of Barnes/Cressman weighting for spatial averaging
typically used in meteorology.  And if you have no idea what I'm
talking about, Google them, and you'll see. :)

Ryan

-- 
Ryan May
Graduate Research Assistant
School of Meteorology
University of Oklahoma
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Vector interpolation on a 2D grid (with realistic results)

2009-11-07 Thread Ryan May
On Sat, Nov 7, 2009 at 5:38 PM, Pierre GM pgmdevl...@gmail.com wrote:
 Linear interpolation with the delaunay package doesn't work great for
 my data. I played with the radial basis functions, but I'm afraid
 they're leading me down the dark, dark path of parameter fiddling. In
 particular, I'm not sure how to prevent my interpolated values to be
 bounded by the min and max of my actual observations.
 Ralf' suggestion of smoothing the values afterwards is tempting, but
 sounds a bit too ad-hoc (btw, Ralf, those were relative differences of
 monthly average precipitation between El Niño and Neutral phases for
 November).

That was me, not Ralf. :)  And I agree, I the interpolated field does
look a bit noisy for such data.  I've been doing the smoothing on top
of natural neighbor for doing some of my own meteorological analysis.
Using the Gaussian kernel isn't really *that* ad hoc considering the
prevalence of Barnes/Cressman weighting for spatial averaging
typically used in meteorology.  And if you have no idea what I'm
talking about, Google them, and you'll see. :)

Ryan

-- 
Ryan May
Graduate Research Assistant
School of Meteorology
University of Oklahoma
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] initializing an array of lists

2009-11-07 Thread alan
I want to build a 2D array of lists, and so I need to initialize the
array with empty lists :

myarray = array([[[],[],[]] ,[[],[],[]]])

Is there a clever way to do this? I could define the array

myarray = zeros( (xdim,ydim), dtype=object)
and then iterate through the elements initializing then to empty lists, but 
surely there is a better way.

-- 
---
| Alan K. Jackson| To see a World in a Grain of Sand  |
| a...@ajackson.org  | And a Heaven in a Wild Flower, |
| www.ajackson.org   | Hold Infinity in the palm of your hand |
| Houston, Texas | And Eternity in an hour. - Blake   |
---
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] initializing an array of lists

2009-11-07 Thread Alan G Isaac
On 11/7/2009 10:56 PM, a...@ajackson.org wrote:
 I want to build a 2D array of lists, and so I need to initialize the
 array with empty lists :

 myarray = array([[[],[],[]] ,[[],[],[]]])


[[[] for i in range(3)] for j in range(2) ]

fwiw,
Alan Isaac

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Use-case for np.choose

2009-11-07 Thread josef . pktd
On Sat, Nov 7, 2009 at 7:53 PM, David Goldsmith d.l.goldsm...@gmail.com wrote:
 Thanks, Anne.

 On Sat, Nov 7, 2009 at 1:32 PM, Anne Archibald peridot.face...@gmail.com
 wrote:

 2009/11/7 David Goldsmith d.l.goldsm...@gmail.com:

 snip


  Also, my experimenting suggests that the index array ('a', the first
  argument in the func. sig.) *must* have shape (choices.shape[-1],) -
  someone
  please let me know ASAP if this is not the case, and please furnish me
  w/ a
  counterexample because I was unable to generate one myself.

 It seems like a and each of the choices must have the same shape

 So in essence, at least as it presently functions, the shape of 'a'
 *defines* what the individual choices are within 'choices`, and if 'choices'
 can't be parsed into an integer number of such individual choices, that's
 when an exception is raised?


 (with

 the exception that choices acn be scalars), but I would consider this
 a bug.

 OK, then we definitely need more people to opine on this, because, if the
 the two don't match, our established policy is to document *desired*
 behavior, not extant behavior (and file a bug ticket).


 Really, a and all the choices should be broadcast to the same
 shape. Or maybe it doesn't make sense to broadcast a - it could be

 Thus begging the question: does anyone actually have an extant, specific
 use-case?


 valuable to know that the result is always exactly the same shape as a
 - but broadcasting all the choice arrays presents an important
 improvement of choose over fancy indexing.

 Then perhaps we need either another function, or a flag specifying which
 behavior this one should exhibit.


 There's a reason choose
 accepts a sequence of arrays as its second argument, rather than a
 higher-dimensional array.

 And that reason is probably supposed to be transparent above, but I've
 confused it by this point, so can you please reiterate it here, in so many
 words. :-)

From looking at a few special cases, I think that full broadcasting rules 
apply.
First a and all choice array are broadcast to the same shape,
then the selection is done according to the elements of (the broadcasted) a.

For broadcasting it doesn't matter whether they are scalars or 1d or 2d or
a 2d single column array. (I haven't tried more than 2 dimensions)

The examples look a bit messy, but broadcasting is relatively straightforward.

(I think, np.where is a bit easier to use because `a` is just a
condition and doesn't
require an index array)

Josef

 np.choose(1, (3,4))
4
 np.choose(0, (3,4))
3
 np.choose(0, (np.arange(3)[:,None],np.arange(4),0))
array([[0, 0, 0, 0],
   [1, 1, 1, 1],
   [2, 2, 2, 2]])
 np.choose(2, (np.arange(3)[:,None],np.arange(4),0))
array([[0, 0, 0, 0],
   [0, 0, 0, 0],
   [0, 0, 0, 0]])
 np.choose(1, (np.arange(3)[:,None],np.arange(4),0))
array([[0, 1, 2, 3],
   [0, 1, 2, 3],
   [0, 1, 2, 3]])
 np.choose([1,2,0,0], (np.arange(3)[:,None],np.arange(4),0))
array([[0, 0, 0, 0],
   [0, 0, 1, 1],
   [0, 0, 2, 2]])
 np.choose(np.array([[1,2,0,0]]), (np.arange(3)[:,None],np.arange(4),0))
array([[0, 0, 0, 0],
   [0, 0, 1, 1],
   [0, 0, 2, 2]])
 np.choose(np.array([[1,2,0]]).T, (np.arange(3)[:,None],np.arange(4),0))
array([[0, 1, 2, 3],
   [0, 0, 0, 0],
   [2, 2, 2, 2]])





 Thanks again,

 DG

 Anne

  Thanks,
 
  DG
 
  ___
  NumPy-Discussion mailing list
  NumPy-Discussion@scipy.org
  http://mail.scipy.org/mailman/listinfo/numpy-discussion
 
 
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion


 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Use-case for np.choose

2009-11-07 Thread Anne Archibald
2009/11/7 David Goldsmith d.l.goldsm...@gmail.com:
 Thanks, Anne.

 On Sat, Nov 7, 2009 at 1:32 PM, Anne Archibald peridot.face...@gmail.com
 wrote:

 2009/11/7 David Goldsmith d.l.goldsm...@gmail.com:

 snip


  Also, my experimenting suggests that the index array ('a', the first
  argument in the func. sig.) *must* have shape (choices.shape[-1],) -
  someone
  please let me know ASAP if this is not the case, and please furnish me
  w/ a
  counterexample because I was unable to generate one myself.

 It seems like a and each of the choices must have the same shape

 So in essence, at least as it presently functions, the shape of 'a'
 *defines* what the individual choices are within 'choices`, and if 'choices'
 can't be parsed into an integer number of such individual choices, that's
 when an exception is raised?

Um, I don't think so.

Think of it this way: you provide np.choose with a selector array, a,
and a list (not array!) [c0, c1, ..., cM] of choices. You construct an
output array, say r, the same shape as a (no matter how many
dimensions it has). The (i0, i1, ..., iN) element of the output array
is obtained by looking at the (i0, i1, ..., iN) element of a, which
should be an integer no larger than M; say j. Then r[i0, i1, ..., iN]
= cj[i0, i1, ..., iN]. That is, each element of the selector array
determines which of the choice arrays to pull the corresponding
element from.

For example, suppose that you are processing an array C, and have
constructed a selector array A the same shape as C in which a value is
0, 1, or 2 depending on whether the C value is too small, okay, or too
big respectively. Then you might do something like:

C = np.choose(A, [-inf, C, inf])

This is something you might want to do no matter what shape A and C
have. It's important not to require that the choices be an array of
choices, because they often have quite different shapes (here, two are
scalars) and it would be wasteful to broadcast them up to the same
shape as C, just to stack them.

Anne
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion