numpy-discussion — Discussion list for all users of Numerical Python

Re: [Numpy-discussion] storage for records

[Stefan v. d. W. <st...@su...>]

On Sat, Feb 18, 2006 at 06:20:50PM -0700, Travis Oliphant wrote:
> Stefan van der Walt wrote:
>=20
> >I am probably trying to do something silly, but still:
> >
> >In [1]: import numpy as N
> >
> >In [2]: N.__version__
> >Out[2]: '0.9.6.2127'
> >
> >In [3]: P =3D N.array(N.zeros((2,2)), N.dtype((('f4',3), {'names':=20
> >['x','y','z'], 'formats': ['f4','f4','f4']})))
> >*** glibc detected *** malloc(): memory corruption: 0x0830bb48 ***
> >Aborted
> >
> >Regards
> >St=E9fan
> >=20
> >
> This code found a bug that's been there for a while in the=20
> PyArray_CastTo code (only seen on multiple copies) which is being done=20
> here as the 2x2 array of zeros is being cast to a 2x2x3 array of=20
> floating-point zeros.
>=20
> The bug should be fixed in SVN, now.

Thank you very much for fixing this! (It works now).

> Despite the use of fields, the base-type is ('f4',3) which is equivalen=
t=20
> to (tack on a 3 to the shape of the array of 'f4').  So, on array=20
> creation the fields will be lost and you will get a 2x2x3 array of=20
> float32.   Types like ('f4', 3) are really only meant to be used in=20
> records.  If they are used "by themselves" they simply create an array=20
> of larger dimension. =20

<snip insightful explanation>

Exactly what I needed for my application!  I'll write this up and put
it on the wiki.

Cheers
St=E9fan

[Numpy-discussion] Re: number ranges (was Re: Matlab page on scipy wiki)

[Colin J. W. <cj...@sy...>]

Bryan Cole wrote:

>>        
>>
>>First, I think the range() function in python is ugly to begin with.
>>Why can't python just support range notation directly like 'for a in
>>0:10'.  Or with 0..10 or 0...10 syntax.  That seems to make a lot more
>>sense to me than having to call a named function.   Anyway, that's a
>>python pet peeve, and python's probably not going to change something
>>so fundamental... 
>>    
>>
>
>There was a python PEP on this. It got rejected as having too many
>'issues'. Pity, in my view.
>
>see http://www.python.org/peps/pep-0204.html
>
>BC
>  
>
This decision appears to have been made nearly six years ago.  It would
be a good idea to revisit the decision, particularly since the reasons
for rejection are not clearly spelled out.

The conditional expression (PEP 308) was rejected but is currently being
implemented in Python version 2.5.  It will have the syntax A if C else B.

I have felt, as Gary Ruben says above, that the range structure is ugly.

Two alternatives have been suggested:
a) a:b:c
b) a..b..c

Do either of these create parsing problems?
for i in a:b:
Should this be treated as an error, with a missing c (the increment)
  print i

for i in 2..5:
  print i                                                           We
don't know whether the 2 is an integer until the second period is scanned.

It would be good if the range and slice could be merged in some way,
although the extended slice is rather complicated - I don't understand it.

    The semantics for an extended slicing are as follows. The primary
    must evaluate to a mapping object, and it is indexed with a key that
    is constructed from the slice list, as follows. If the slice list
    contains at least one comma, the key is a tuple containing the
    conversion of the slice items; otherwise, the conversion of the lone
    slice item is the key. The conversion of a slice item that is an
    expression is that expression. The conversion of an ellipsis slice
    item is the built-in |Ellipsis| object. The conversion of a proper
    slice is a slice object (see section section 4.2 The standard type
    hierarchy <http://www.network-theory.co.uk/docs/pylang/ref_30.html>)
    whose |start|, |stop| and |step| attributes are the values of the
    expressions given as lower bound, upper bound and stride,
    respectively, substituting |None| for missing expressions.

[source: http://www.network-theory.co.uk/docs/pylang/ref_60.html]

The seems to be a bit of a problem with slicing that needs sorting out.

The syntax for a slice list appears to allow multiple slices in a list:

    extended_slicing 	::= 	primary <primaries.html#tok-primary> "["
    slice_list <slicings.html#tok-slice_list> "]"
    slice_list 	::= 	slice_item <slicings.html#tok-slice_item> (","
    slice_item <slicings.html#tok-slice_item>)* [","]

but the current interpreter reports an error:

     >>> a= range(20)
     >>> a[slice(3, 9, 2)]
    [3, 5, 7]
     >>> a[slice(3, 9, 2), slice(5, 10)]
    Traceback (most recent call last):
      File "<interactive input>", line 1, in ?
    TypeError: list indices must be integers
     >>>

I have taken the liberty of cross posting this to c.l.p.

Colin W.

Re: [Numpy-discussion] storage for records

[Travis O. <oli...@ie...>]

Stefan van der Walt wrote:

>I am probably trying to do something silly, but still:
>
>In [1]: import numpy as N
>
>In [2]: N.__version__
>Out[2]: '0.9.6.2127'
>
>In [3]: P = N.array(N.zeros((2,2)), N.dtype((('f4',3), {'names': ['x','y','z'], 'formats': ['f4','f4','f4']})))
>*** glibc detected *** malloc(): memory corruption: 0x0830bb48 ***
>Aborted
>
>Regards
>Stéfan
>  
>
This code found a bug that's been there for a while in the 
PyArray_CastTo code (only seen on multiple copies) which is being done 
here as the 2x2 array of zeros is being cast to a 2x2x3 array of 
floating-point zeros.

The bug should be fixed in SVN, now.

Despite the use of fields, the base-type is ('f4',3) which is equivalent 
to (tack on a 3 to the shape of the array of 'f4').  So, on array 
creation the fields will be lost and you will get a 2x2x3 array of 
float32.   Types like ('f4', 3) are really only meant to be used in 
records.  If they are used "by themselves" they simply create an array 
of larger dimension.  

By the way,  the N.dtype in the array constructor is unnecessary as that 
is essentially what is done to the second argument anyway

You can get two different views of the same data (which it seems you are 
after) like this:

P = N.zeros((2,2), {'names': ['x','y','z'], 'formats': ['f4','f4','f4']})
Q = P.view(('f4',3))

Then

Q[...,0] = 10
print P['x']

If you want the field to vary in the first dimension, then you really 
want a FORTRAN array.  So,

P = N.zeros((2,2), {'names': ['x','y','z'], 'formats': 
['f4','f4','f4']},fortran=1)
Q = P.view(('f4',3))

Then

Q[0] = 20
print P['x']


Best,

-Travis

Re: [Numpy-discussion] optimizing power() for complex and real cases

[Tim H. <tim...@co...>]

OK, I now have a faily clean implementation in C of:

def __pow__(self, p):
    if p is not a scalar:
       return power(self, p)
    elif p == 1:
       return p
    elif p == 2:
       return square(self)
#    elif p == 3:
#       return cube(self)
#    elif p == 4:
#       return power_4(self)
#    elif p == 0:
#       return ones(self.shape, dtype=self.dtype)
#    elif p == -1:
#       return 1.0/self
    elif p == 0.5:
       return sqrt(self)


First a couple of technical questions, then on to the philosophical portion of this note.

1. Is there a nice fast way to get a matrix filled with ones from C. I've been tempted to write a ufunc 'ones_like', but I'm afraid that might be considered inappropriate.

2. Are people aware that array_power is sometimes passed non arrays as its first argument? Despite having the signature:

array_power(PyArrayObject *a1, PyObject *o2) 

This caused me almost no end of headaches, not to mention crashes during numpy.test().

I'll check this into the power_optimization branch RSN, hopefully with a fix for the zero power case. Possibly also after extending it to inplace power as well.


OK, now on to more important stuff. As I've been playing with this my opinion has gone in circles a couple of times. I now think the issue of optimizing integer powers of complex numbers and integer powers of floats are almost completely different. Because complex powers are quite slow and relatively inaccurate, it is appropriate to optimize them for integer powers at the level of nc_pow. This should be just a matter of liberal borrowing from complexobject.c, but I haven't tried it yet.

On the other hand, real powers are fast enough that doing anything at the single element level is unlikely to help. So in that case we're left with either optimizing the cases where the dimension is zero as David has done, or optimizing at the __pow__ (AKA array_power) level as I've done now based on David's original suggestion. This second approach is faster because it avoids the mysterious python scalar -> zero-D array conversion overhead. However, it suffers if we want to optimize lots of different powers since one needs a ufunc for each one. So the question becomes, which powers should we optimize?

My latest thinking on this is that we should optimize only those cases where the optimized result is no less accurate than that produced by pow. I'm going to assume that all C operations are equivalently accurate, so pow(x,2) has roughly the same amount of error as x*x. (Something on the order of 0.5 ULP I'd guess). In that case:
   pow(x, -1) -> 1 / x
   pow(x, 0) -> 1
   pow(x, 0.5) -> sqrt(x)
   pow(x, 1) -> x
   pow(x, 2) -> x*x
can all be implemented in terms of multiply or divide with the same accuracy as the original power methods. Once we get beyond these, the error will go up progressively.

The minimal set described above seems like it should be relatively uncontroversial and it's what I favor. Once we get beyond this basic set, we would need to reach some sort of consensus on how much additional error we are willing to tolerate for optimizing these extra cases. You'll notice that I've changed my mind, yet again, over whether to optimize A**0.5. Since the set of additional ufuncs needed in this case is relatively small, just square and inverse (==1/x), this minimal set works well if optimizing in pow as I've done.

That's the state of my thinking on this at this exact moment. I'd appreciate any comments and suggestions you might have.

[Numpy-discussion] A case for rank-0 arrays

[Sasha <nd...@ma...>]

I have reviewed mailing list discussions of rank-0 arrays vs. scalars
and I concluded that the current implementation that contains both is
(almost) correct.  I will address the "almost" part with a concrete
proposal at the end of this post (search for PROPOSALS if you are only
interested in the practical part).

The main criticism of supporting both scalars and rank-0 arrays is
that it is "unpythonic" in the sense that it provides two almost
equivalent ways to achieve the same result.  However, I am now
convinced that this is the case where practicality beats purity.

If you take the "one way" rule to it's logical conclusion, you will
find that once your language has functions, it does not need numbers
or any other data type because they all can be represented by
functions (see http://en.wikipedia.org/wiki/Church_numeral).   Another
example of core python violating the "one way rule" is the presence of
scalars and length-1 tuples.  In S+, for example, scalars are
represented by single element lists.

The situation with ndarrays is somewhat similar. A rank-N array is
very similar to a function with N arguments, where each argument has a
finite domain (i-th domain of a is range(a.shape[i])).  A rank-0 array
is just a function with no arguments and as such it is quite different
from a scalar.  Just as a function with no arguments cannot be
replaced by a constant in the case when a value returned may change
during the run of the program, rank-0 array cannot be replaced by an
array scalar because it is mutable.  (See
http://projects.scipy.org/scipy/numpy/wiki/ZeroRankArray for use
cases).

Rather than trying to hide rank-0 arrays from the end-user and treat
it as an implementation artifact, I believe numpy should emphasize the
difference between rank-0 arrays and scalars and have clear rules on
when to use what.

PROPOSALS
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D

Here are three suggestions:

1. Probably the most controversial question is what getitem should
return. I believe that most of the confusion comes from the fact that
the same syntax implements two different operations: indexing and
projection (for the lack of better name).  Using the analogy between
ndarrays and functions, indexing is just the application of the
function to its arguments and projection is the function projection
((f, x) -> lambda (*args): f(x, *args)).

The problem is that the same syntax results in different operations
depending on the rank of the array.

Let

>>> x =3D ones((2,2))
>>> y =3D ones(2)

then x[1] is projection and type(x[1]) is ndarray, but y[1] is
indexing and type(y[1]) is int32.  Similarly, y[1,...] is indexing,
while x[1,...] is projection.

I propose to change numpy rules so that if ellipsis is present inside
[], the operation is always projection and both y[1,...] and
x[1,1,...] return zero-rank arrays.  Note that I have previously
rejected Francesc's idea that x[...] and x[()] should have different
meaning for zero-rank arrays.  I was wrong.

2. Another source of ambiguity is the various "reduce" operations such
as sum or max.  Using the previous example, type(x.sum(axis=3D0)) is
ndarray, but type(y.sum(axis=3D0)) is int32.  I propose two changes:

   a. Make x.sum(axis)  return ndarray unless axis is None, making
type(y.sum(axis=3D0)) is ndarray true in the example.

   b. Allow axis to be a sequence of ints and make
x.sum(axis=3Drange(rank(x))) return rank-0 array according to the rule
2.a above.

   c. Make x.sum() raise an error for rank-0 arrays and scalars, but
allow x.sum(axis=3D()) to return x.  This will make numpy sum consistent
with the built-in sum that does not work on scalars.

3. This is a really small change currently

>>> empty(())
array(0)

but

>>> ndarray(())
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
ValueError: need to give a valid shape as the first argument


I propose to make shape=3D() valid in ndarray constructor.

Re: [Numpy-discussion] maxima of masked arrays

[Sasha <nd...@ma...>]

On 2/17/06, Sasha <nd...@ma...> wrote:
> On 2/17/06, George Nurser <ag...@no...> wrote:
> > But with a masked array there are problems.
>
> What you see is a bug in ma.

Fixed in SVN.

Re: [Numpy-discussion] Behavior of array scalars

[Sasha <nd...@ma...>]

On 2/17/06, Travis Oliphant <oli...@ie...> wrote:
> Sorry we worked on the same code.
Not a problem.  My code was just a proof of concept anyway.

>  I already comitted a changed that solves the problem.
I've seen your change on the timeline.  Thanks for extensive comments.

> It doesn't change the C-API function but instead
> changes the _find_wrap code which needed changing anyway ...

You are right, my patch could be fooled by an output arg with an
__array_wrap__. However, I am not sure calling output argument's =20
__array_wrap__ is a good idea: it looks like it may lead to "(a is
add(a, 2, a)) =3D=3D False)" in some circumstances.  Another concern is
that it looks like what output arguments are supplied is determined
twice: in _find_wrap and in construct_matrices.

Re: [Numpy-discussion] maxima of masked arrays

[Sasha <nd...@ma...>]

On 2/17/06, George Nurser <ag...@no...> wrote:
> But with a masked array there are problems.

What you see is a bug in ma.

> Are there any workarounds for this?

For now you can use

>>> ma.maximum.reduce(amask, 0)
array(data =3D
 [ 9 10 11 12],
      mask =3D
 [False False False False],
      fill_value=3D999999)

Re: [Numpy-discussion] Behavior of array scalars

[Travis O. <oli...@ie...>]

Sasha wrote:

>On 2/17/06, Travis Oliphant <oli...@ee...> wrote:
>  
>
>>>>from numpy 
>>>>
>
>The patch passes all the tests, but I would like to hear from others
>before commit.  Personally, I am unhappy that I had to change C-API
>function.
>  
>
Sorry we worked on the same code.  I already comitted a changed that 
solves the problem.  It doesn't change the C-API function but instead 
changes the _find_wrap code which needed changing anyway so that other 
objects passed in as output arrays work cause the returned object to be

obj.__array_wrap__(<ndarray_ufunc_result>)

no matter what the array priority was.

Now. 

a += b

doesn't change type no matter what a is.

-Travis

Re: [Numpy-discussion] Behavior of array scalars

[Sasha <nd...@ma...>]

On 2/17/06, Travis Oliphant <oli...@ee...> wrote:
> Sasha wrote:
>
> >Maybe we can change ufunc logic so that when the output argument is
> >supplied it is returned without scalar conversion.
> >
> >
> That seems sensible.

Attached patch implements this idea.  With the patch applied:

>>> from numpy import *
>>> x =3D array(5)
>>> add(x,5,x)
array(10)
>>> x+=3D5
>>> x
array(15)

The patch passes all the tests, but I would like to hear from others
before commit.  Personally, I am unhappy that I had to change C-API
function.

Re: [Numpy-discussion] Behavior of array scalars

[Sasha <nd...@ma...>]

On 2/17/06, Christopher Barker <Chr...@no...> wrote:
> >>>>x[()] =3D 10
>
> I can't do that. HAs that been added since version: '0.9.2' ?

Yes, you need 0.9.4 or later.

Re: [Numpy-discussion] Behavior of array scalars

[Travis O. <oli...@ee...>]

Christopher Barker wrote:

> Travis Oliphant wrote:
>
>>  They just don't *stay* 0-d arrays and are converted to array scalars 
>> at almost every opportunity....
>
>
> I'm still confused as to what the difference is. This (recent) 
> convesation started with my desire for a mutable scalar. CAn array 
> scalars fill this role?

No.  array scalars are immutable (well except for the void array scalar...)

>
> What I mean by that role is some way to do:
>
> x += 5 # (and friends)


This now works in SVN.  In-place operations on 0-d arrays don't change 
on you.

-Travis

Re: [Numpy-discussion] Behavior of array scalars

[Christopher B. <Chr...@no...>]

Travis Oliphant wrote:

>  They just don't 
> *stay* 0-d arrays and are converted to array scalars at almost every 
> opportunity....

I'm still confused as to what the difference is. This (recent) 
convesation started with my desire for a mutable scalar. CAn array 
scalars fill this role?

What I mean by that role is some way to do:

x += 5 # (and friends)
x[()] = 45 # or some other notation


And have x be the same object throughout.

Heck even something like:

x.set(45)

would work for me.

Alexander Belopolsky wrote:
>>>>x[...] = 10
> 
> or
> 
>>>>x[()] = 10

I can't do that. HAs that been added since version: '0.9.2' ?

-Chris



-- 
Christopher Barker, Ph.D.
Oceanographer
                                     		
NOAA/OR&R/HAZMAT         (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chr...@no...

Re: [Numpy-discussion] Behavior of array scalars

[Travis O. <oli...@ee...>]

Sasha wrote:

>It would be nice to collect the motivations behind the current state
>of affairs with rank-0 arrays in one place.  Due to the "hard-hat"
>nature of the issue, I would suggest to do it at
>http://projects.scipy.org/scipy/numpy/wiki/ZeroRankArray .
>
>Travis' Numeric3 design document actually leaves the issue open
>
>  
>
This document is old.  Please don't refer to it too stringently.   It 
reflected my thinking at the start of the project.  There are mailing 
list discussions that have more relevance.   The source reflects what 
was actually done.  What was done is introduce scalar array types for 
every data-type and return those.

I had originally thought that the pure Python user would *never* see 
rank-0 arrays.  That's why PyArray_Return is called all over the place 
in the code.  The concept that practicality beats purity won out and 
there are a few limited wasy you can get zero-dimensional arrays (i.e. 
using array(5) which used to return an array scalar).  They just don't 
*stay* 0-d arrays and are converted to array scalars at almost every 
opportunity....  I have been relaxing this over time, however.  

I can't say I have some grand understanding that is guiding the 
relaxation of this rule, however, except that I still think array 
scalars are *better* to deal with (I think this will be especially 
obvious when we get scalar math implemented).   So, I relunctantly give 
visibility to 0-d arrays when particular use-cases emerge.

>In any case I will collect all these thoughts on the ZeroRankArray
>page unless I hear that this belongs to the main wiki.
>  
>
It's a good start.   This particular use case of course is actually 
showing us a deeper flaw in our use of output arguments in the ufunc 
which needs changing.

-Travis

Re: [Numpy-discussion] Behavior of array scalars

[Sasha <nd...@ma...>]

Sorry for a truncated post.  Here is what I intended.

On 2/17/06, Travis Oliphant <oli...@ee...> wrote:
> NumPy (starting with Numeric) has always had this love-hate relationship
> with zero-dimensional arrays.  We use them internally to simplify the
> code, but try not to expose them to the user.  Ultimately, we couldn't
> figure out how to do that cleanly and so we have the current compromise
> situation where 0-d arrays are available but treated as second-class
> citizens.   Thus, we still get funny behavior in certain circumstances.

It would be nice to collect the motivations behind the current state
of affairs with rank-0 arrays in one place.  Due to the "hard-hat"
nature of the issue, I would suggest to do it at
http://projects.scipy.org/scipy/numpy/wiki/ZeroRankArray .

Travis' Numeric3 design document actually leaves the issue open

"""
What does single element indexing return? Scalars or rank-0 arrays?
Right now, a scalar is returned if there is a direct map to a Python
type, otherwise a rank-0 array (Numeric scalar) is returned. But, in
problems which reduce to an array of arbitrary size, this can lead to
a lot of code that basically just checks to see if the object is a
scalar. There are two ways I can see to solve this: 1) always return
rank-0 arrays (never convert to Python scalars) and 2) always use
special functions (like alen) that handle Python scalars correctly.
I'm open to both ideas, but probably prefer #1 (never convert to
Python scalars) unless requested.
""" <http://web.archive.org/web/20050207124620/www.scipy.org/wikis/numdesig=
n/>

I can think of two compelling reasons in favor of scalar array types:

1. Rank-0 arrays cannot be used as indices to tuples.
2. Rank-0 arrays cannot be used as keys in dicts.

Neither of these resons is future proof.  It looks like python 2.5
will introduce __index__ slot that will fix #1 and #2 is probably
better solved by introduction of "frozen" ndarray.

In any case I will collect all these thoughts on the ZeroRankArray
page unless I hear that this belongs to the main wiki.

Re: [Numpy-discussion] Behavior of array scalars

[<co...@ph...>]

Travis Oliphant <oli...@ee...> writes:

> Consider...
>
> a = rand(5,5)
> b = mat(a)
>
> a += b
>
> What do you think the type of a now is?  What should it be?
>
> Currently,  this code would change a from an array to a matrix because
>
> add(a,b,a) returns a matrix.
>
> I'm thinking that we should establish the rule that if output arrays
> are given, then what is returned should just be those output arrays...

+1.

-- 
|>|\/|<
/--------------------------------------------------------------------------\
|David M. Cooke                      http://arbutus.physics.mcmaster.ca/dmc/
|co...@ph...

Re: [Numpy-discussion] Behavior of array scalars

[Travis O. <oli...@ee...>]

Sasha wrote:

>It is actually being translated to "a = add(a,10,a)" by virtue of
>array_inplace_add supplied in the inplace_add slot.  Here is the
>proof:
>  
>

I think we do need to fix something.  Because the problem is even more 
apparent when you in-place add
an array to a matrix.

Consider...

a = rand(5,5)
b = mat(a)

a += b

What do you think the type of a now is?  What should it be?

Currently,  this code would change a from an array to a matrix because

add(a,b,a) returns a matrix.

I'm thinking that we should establish the rule that if output arrays are 
given, then what is returned should just be those output arrays...

This seems to make consistent sense and it will make the inplace 
operators work as expected (not changing the type).  We are currently 
not letting in-place operators change the data-type.  Our argument 
against that behavior is weakened if we do let them change the Python 
type....

-Travis

Re: [Numpy-discussion] Behavior of array scalars

[Travis O. <oli...@ee...>]

Sasha wrote:

>Maybe we can change ufunc logic so that when the output argument is
>supplied it is returned without scalar conversion.
>  
>
That seems sensible.  Any objections?   It is PyArray_Return that 
changes things from 0-d array's to scalars.  It's all that function has 
every really done....

Notice that this behavior was always in Numeric...

a = Numeric.array(5)
a += 10
type(a)

<type 'int'>

 >>> a = Numeric.array(5)
 >>> type(a) == type(a+a)
False

But...

 >>> a = Numeric.array(5,'f')
 >>> type(a) == type(a+a)
True

So, we've been dealing with these issues (poorly) for a long time....

-Travis

Re: [Numpy-discussion] Behavior of array scalars

[Sasha <nd...@ma...>]

On 2/17/06, Travis Oliphant <oli...@ee...> wrote:
> ...
> Perhaps what is going on is that
> a +=3D 10
> is begin translated to
> a =3D a + 10
> rather than
> add(a,10,a)
> I'll have to look deeper to see why.

It is actually being translated to "a =3D add(a,10,a)" by virtue of
array_inplace_add supplied in the inplace_add slot.  Here is the
proof:

>>> a =3D array(0)
>>> a =3D b =3D array(0)
>>> a +=3D 10
>>> b
array(10)
>>> a
10

Another way to explain it is to note that a +=3D 10 is equivalent to "a
=3D a.__iadd__(10)" and a.__iadd__(10) is equivalent to "add(a, 10, a)".

This is not easy to fix because the real culprit is

>>> a =3D array(0)
>>> type(a) is type(a+a)
False

Maybe we can change ufunc logic so that when the output argument is
supplied it is returned without scalar conversion.

Re: [Numpy-discussion] Behavior of array scalars

[Travis O. <oli...@ee...>]

Christopher Barker wrote:

> Hi all,
>
> It just dawned on my that the numpy array scalars might give something 
> I have wanted once in a while: mutable scalars. However, it seems that 
> we almost, but no quite, have them. A few questions:

NumPy (starting with Numeric) has always had this love-hate relationship 
with zero-dimensional arrays.  We use them internally to simplify the 
code, but try not to expose them to the user.  Ultimately, we couldn't 
figure out how to do that cleanly and so we have the current compromise 
situation where 0-d arrays are available but treated as second-class 
citizens.   Thus, we still get funny behavior in certain circumstances.  
I think you found another such quirky area. 

I'm open to suggestions.


To analyze this particular case...

The a+= 10

operation should be equivalent to add(a,10,a). 

Note that explicitly writing

add(a,10,a)

returns a scalar (all ufuncs return scalars if 0-d arrays are the 
result).   But, a is modified in-place as you wanted.

Perhaps what is going on is that

a += 10

is begin translated to

a = a + 10

rather than

add(a,10,a)

I'll have to look deeper to see why.

-Travis

Re: [Numpy-discussion] optimizing power() for complex and real cases

[Tim H. <tim...@co...>]

David M. Cooke wrote:

>Tim Hochberg <tim...@co...> writes:
>
>  
>
>>Here's a little progress report: I now have A**2 running as fast as
>>square(A). This is done by special casing stuff in array_power so that
>>A**2 acutally calls square(A) instead of going through power(A,2).
>>Things still need a bunch of cleaning up (in fact right now A**1
>>segfaults, but I know why and it should be an easy fix). However, I
>>think I've discovered why you couldn't get your special cased power to
>>run as fast for A**2 as square(A) or A*A. It appears that the overhead
>>of creating a new array object from the integer 2 is the bottleneck. I
>>was running into the same mysterious overhead, even when dispatching
>>from array_power, until I special cased on PyInt to avoid the array
>>creation in that case.
>>    
>>
>
>Hmm, if that's true about the overhead, that'll hit all computations
>of the type op(x, <some scalar>). Something to look at. That ufunc code
>for casting the arguments is pretty big and hairy, so I'm not going to
>look at right now :-)
>  
>
Well, it's just a guess based on the fact that the extra time went away 
when I stopped calling PyArray_EnsureArray(o2) for python ints. For what 
it's worth, numpy scalars seem to have much less overhead. As indicated 
below (note that numpy scalars are not currently special cased like 
PyInts are). The overhead from PyInts was closer to 75% versus about 15% 
for numpy scalars. Of course, the percentage of overhead is going to go 
up for smaller arrays.

 >>> Timer('a**2', 'from numpy import arange;a = arange(10000.); b = 
a[2]').timeit(10000)
0.28345055027943999
 >>> Timer('a**b', 'from numpy import arange;a = arange(10000.); b = 
a[2]').timeit(10000)
0.32190487897204889
 >>> Timer('a*a', 'from numpy import arange;a = arange(10000.); b = 
a[2]').timeit(10000)
0.27305732991204223
 >>> Timer('square(a)', 'from numpy import arange, square;a = 
arange(10000.); b = a[2]').timeit(10000)
0.27989618792332749

-tim

Re: [Numpy-discussion] Behavior of array scalars

[Sasha <nd...@ma...>]

On 2/17/06, Christopher Barker <Chr...@no...> wrote:
>  >>> x +=3D 5
>
> I expect this to change the object in place.
>
>  >>> x
> 10
>
> but what is this? is it no longer an array?

I would say it is a bug, but here is an easy work-around

>>> x =3D array(5)
>>> id(x)
6425088
>>> x[()]+=3D5
>>> id(x)
6425088
>>> x
array(10)

You can also use

>>> x[...]+=3D5
>>> x
array(15)

With an additional benefit that the same syntax works for any shape.

> Is there a way to set the value in place, without resorting to:


>>> x[...] =3D 10
or
>>> x[()] =3D 10


You can see more on this feature at
http://projects.scipy.org/scipy/numpy/wiki/ZeroRankArray

Re: [Numpy-discussion] Behavior of array scalars

[Alexander B. <ale...@gm...>]

On 2/17/06, Christopher Barker <Chr...@no...> wrote:
>  >>> x +=3D 5
>
> I expect this to change the object in place.
>
>  >>> x
> 10
>
> but what is this? is it no longer an array?

I would say it is a bug, but here is an easy work-around

>>> x =3D array(5)
>>> id(x)
6425088
>>> x[()]+=3D5
>>> id(x)
6425088
>>> x
array(10)

You can also use

>>> x[...]+=3D5
>>> x
array(15)

With an additional benefit that the same syntax works for any shape.

> Is there a way to set the value in place, without resorting to:


>>> x[...] =3D 10
or
>>> x[()] =3D 10


You can see more on this feature at
http://projects.scipy.org/scipy/numpy/wiki/ZeroRankArray

[Numpy-discussion] Behavior of array scalars

[Christopher B. <Chr...@no...>]

Hi all,

It just dawned on my that the numpy array scalars might give something I 
have wanted once in a while: mutable scalars. However, it seems that we 
almost, but no quite, have them. A few questions:

 >>> import numpy as N
 >>> N.__version__
'0.9.2'
 >>> N.array(5)
array(5)
 >>>
 >>> x = N.array(5)
 >>> x.shape
()

So it looks like a scalar.
 >>> y = x

Now I have two names bound to the same object.

 >>> x += 5

I expect this to change the object in place.

 >>> x
10

but what is this? is it no longer an array?

 >>> y
array(10)

y changed, so it looks like the object has changed in place.

 >>> type(x)
<type 'int32_arrtype'>
 >>> type(y)
<type 'numpy.ndarray'>

So why did x += 5 result in a different type of object?

Also:

I can see that we could use += and friends to mutate an array scalar, 
but what if I want to set it's value, as a mutation, like:

 >>> x = N.array((5,))
 >>> x
array([5])
 >>> x[0] = 10
 >>> x
array([10])

but I can't so that with an array scalar:

 >>> x = N.array(5)
 >>> x
array(5)
 >>> x[0] = 10
Traceback (most recent call last):
   File "<stdin>", line 1, in ?
IndexError: 0-d arrays can't be indexed.
 >>> x[] = 10
   File "<stdin>", line 1
     x[] = 10
       ^
SyntaxError: invalid syntax
 >>> x[:] = 10
Traceback (most recent call last):
   File "<stdin>", line 1, in ?
ValueError: cannot slice a scalar

Is there a way to set the value in place, without resorting to:

 >>> x *= 0
 >>> x += 34

I think it would be really handy to have a full featured, mutable scalar.

-Chris









-- 
Christopher Barker, Ph.D.
Oceanographer
                                     		
NOAA/OR&R/HAZMAT         (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chr...@no...

[Numpy-discussion] Re: number ranges (was Re: Matlab page on scipy wiki)

[Bryan C. <br...@co...>]

>         
> 
> First, I think the range() function in python is ugly to begin with.
> Why can't python just support range notation directly like 'for a in
> 0:10'.  Or with 0..10 or 0...10 syntax.  That seems to make a lot more
> sense to me than having to call a named function.   Anyway, that's a
> python pet peeve, and python's probably not going to change something
> so fundamental... 

There was a python PEP on this. It got rejected as having too many
'issues'. Pity, in my view.

see http://www.python.org/peps/pep-0204.html

BC