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

Re: [Numpy-discussion] questions regarding assignement and copy

[Robert K. <rob...@gm...>]

David Cournapeau wrote:

>     bar   += 1
>     print bar is foo
> 
>     prints True
> 
>     But if I do bar = bar + 1, then bar is not a copy of foo anymore. Is 
> this intended ?

Yes, very much so! If it were otherwise Guido would never have added the 
augmented assignment operators. From the reference manual:

http://docs.python.org/ref/augassign.html

"""An augmented assignment expression like x += 1 can be rewritten as x = x + 1 
to achieve a similar, but not exactly equal effect. In the augmented version, x 
is only evaluated once. Also, when possible, the actual operation is performed 
in-place, meaning that rather than creating a new object and assigning that to 
the target, the old object is modified instead."""

> This looks really confusing to me, and I would like to 
> know what the precise rules about copy vs alias are ?

I'm not going to go through all of the functions in numpy, but here are some 
(unorganized) rules of thumb:

* Any operations with regular operators + - / * // ** ~ ^ & | << >> < > == <=, 
etc. will put their results into new arrays. If the expression is on the RHS of 
a regular assignment, just ignore the LHS of the assignment; even if the name on 
the LHS is the same as a name on the RHS, the expression is evaluated, the new 
array is created and *then* the new array is assigned to the LHS name.

* Augmented assignment operators will operate in-place.

* Functions and methods should state that they operate in-place in their docstrings.

* The asarray(), asanyarray(), etc. functions will return the input array if it 
is already suitable.

* Slicing, reshaping, transposing and similar structural operations create new 
array objects, but they are generally views, not copies of the data.

-- 
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

Re: [Numpy-discussion] sort(axis=-1)

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

Christian Kristukat wrote:
> Hi,
>
> it seems that -1 as axis parameter is interpreted like in array indexing, i.e.
> -1 means the last axis rather than meaning the only axis of the flattened
> representation, like for example with take().
> Is that intendend?
>
>   

I don't understand the question.   Yes, sort(axis=-1) sorts along the 
last dimension.

AFAIK, take has the same meaning for axis.

-Travis

Re: [Numpy-discussion] sort(axis=-1)

[Robert K. <rob...@gm...>]

Christian Kristukat wrote:
> Hi,
> 
> it seems that -1 as axis parameter is interpreted like in array indexing, i.e.
> -1 means the last axis rather than meaning the only axis of the flattened
> representation, like for example with take().
> Is that intendend?

take(axis=-1) also goes over the last axis. axis=None operates over the 
flattened array. axis=None is the default for take(), but axis=-1 is default for 
sort() since axis=None doesn't make much sense for an inplace operation and is 
not permitted.

-- 
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

Re: [Numpy-discussion] questions regarding assignement and copy

[Charles R H. <cha...@gm...>]

On 10/18/06, David Cournapeau <da...@ar...> wrote:
>
> Sven Schreiber wrote:
> >
> > Yes it's intended; as far as I understand the python/numpy syntax, <+>
> > is an operator, and that triggers assignment by copy (even if you do
> > something trivial as bar = +foo, you get a copy, if I'm not mistaken),
> >
> So basically, whenever you have
>
> foo = expr
>
> with expr is a numpy expression containing foo, you trigger a copy ?


Yes. Numpy generates a temporary where the expression results are kept, the
temporary is then assigned to foo. This can be inefficient, but knowing when
you may overwrite data can be complicated. For instance:

foo = dot(foo,foo)

Can't be done in place because values needed later on in the dot evaluation
would be overwritten.

Chuck

Re: [Numpy-discussion] The NumPy Fortran-ordering quiz

[Charles R H. <cha...@gm...>]

On 10/17/06, Travis Oliphant <oli...@ie...> wrote:
>
> Stefan van der Walt wrote:
> > One last case, which confuses me still (probably because it is
> > 04:16am):
>
> Please continue to question.  All the code needs as much review as it
> can get.


I am really starting to wonder why we need an order keyword at all except
when order matters for interfacing, i.e., in contiguos arrays where one
wants fortran contiguous or c contiguous. For output type stuff:

tofile
tostring
ravel
flatten

These all need the keyword so that they can be used to produce files and
arrays to interface to fortran. Right now the fortran interface can be
achieved by:

a.T.tofile
a.T.tostring
a.T.ravel
a.T.flatten

The order keyword is just syntatic sugar that makes the intent clearer. For
the input type stuff

fromfile
fromstring
reshape-1d-array (unravel?)

Currently, the key operation is reshape, which only needs to return a view
in fortran order and doesn't even need to mark the resulting array as
fortran order because, well, because it works just fine in numpy as is, it
just isn't contiguous. If the other functions took shape and order, reshape
wouldn't even need the order keyword.

I don't see why the array constructor needs the order keyword, it doesn't
*do* anything. For instance

a = array([[1,2,3],[4,5,6]], order='F')

doesn't produce a fortran contiguous array, it produces the same array as
the 'C' form, just sets the fortran flag and marks contiguous as False. What
is the use of that? It is just a generic non-contiguous numpy array. And

In [131]: ascontiguousarray(array([[1,2,3],[4,5,6]], dtype=int8,
order='F')).flags
Out[131]:
  CONTIGUOUS : True
  FORTRAN : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False

Doesn't produce a fortran contiguous array, so what use was the flag? And

In [141]: array([1,2,3,4,5,6], dtype=int8).reshape((2,3),
order='F').astype(int16).flags
Out[141]:
  CONTIGUOUS : True
  FORTRAN : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False

reorders stuff in memory, so is a bug looking to happen in a fortran
interface.

mmapped files are the only thing I can think of where one might want vary an
operation depending on Fortran ordering because seeking out of order is very
expensive. But that means adapting algorithms depending on order type,
better I think to just stick to using the small strided dimensions when
appropriate.

It would be helpful in debugging all this order stuff if it was clear what
was supposed to happen in every case. Ravel, for instance, ignores the
FORTRAN flag, again begging the question as to why we *have* the flag.

Chuck

[Numpy-discussion] No forums found for Numerical Python

[Ray S. <ra...@bl...>]

FYI:
at http://numpy.org/, the link "Forums" (http://sourceforge.net/forum/?group_id=1369) gives

No forums found for Numerical Python

!

Re: [Numpy-discussion] questions regarding assignement and copy

[A. M. A. <per...@gm...>]

On 18/10/06, David Cournapeau <da...@ar...> wrote:
> Sven Schreiber wrote:
> >
> > Yes it's intended; as far as I understand the python/numpy syntax, <+>
> > is an operator, and that triggers assignment by copy (even if you do
> > something trivial as bar = +foo, you get a copy, if I'm not mistaken),
> >
> So basically, whenever you have
>
> foo = expr
>
> with expr is a numpy expression containing foo, you trigger a copy ?

No. "=" never copies, but "+" does: when you do "A+B" you produce a
new, freshly-allocated array, which you can then store (a reference
to) in any variable you like, including A or B. So

M = M+1

makes a copy because "M+1" is a new matrix, which you are assigning to
M. Unfortunately, this means if you do:

M = 2*(M+1)

python makes the new matrix M+1, then makes the new matrix 2*(M+1),
discards M+1, and sets M to point to 2*(M+1). If you want to avoid all
this copying, you can use python's in-place operators:

M += 1
M *= 2

This is actually a standard difficulty people have with python, made
more obvious because you're working with mutable arrays rather than
immutable scalars.

A. M. Archibald

Re: [Numpy-discussion] The NumPy Fortran-ordering quiz

[A. M. A. <per...@gm...>]

On 18/10/06, Travis Oliphant <oli...@ie...> wrote:
>
> I just committed a change to the check-for-copy code to SVN.   A copy
> will occur if an actual reshaping is taking place that involves more
> than just adding or deleting ones from the old shape and if
>
> 1) self is not "single-segment".  Single-segment means either
> Fortran-order or C-order.  If we don't have a single-segment array, then
> we have to get a (Fortran or C-) contiguous buffer to do any reshaping
> (there may be special cases when it's possible, but I haven't figure
> them out...)

Indeed it is often possible to reshape an array that is not
contiguous; for example zeros(100)[::10] is not contiguous but can be
reshaped to any shape without any copying.

More generally, there are all sorts of peculiar arrangements that can
be reshaped witout a copy, for example:

input has lengths (8,2,3) and strides (39,9,3); output should have
lengths (4,4,3,2), so take strides (156,39,6,3)

Whether this sort of thing actually *occurs* very often, I don't know...

> 2) self is single-segment but self.squeeze().ndim > 1 and it is in the
> "wrong" order from what was requested in the order argument (i.e. self
> is in C-contiguous order but a Fortran-order reshape was requested or
> self is in F-contiguous order but a C-order reshape was requested).
>
> If there are any cases satisfying these rules where a copy does not have
> to occur then let me know.

There are. I came up with a general condition for describing exactly
when you need to copy an array, but it's difficult to put into words.
Here goes:

(0) Dimensions of length 1 are an irrelevant annoyance. Conceptually
they should be squeezed out before beginning and newaxised back in at
the end. (What strides should they have?) Assume that I've done so
from now on. (Actual code can just skip them appropriately.)

(1) If the array looks like an evenly-strided 1D array that has been
reshaped (so strides[i+-1]=lengths[i]*strides[i] for all i) it can
always be reshaped without a copy. (For example, zeros(1000)[::10] can
be reshaped arbitrarily and as many times as desired without copying.)

(2) If the array does *not* look like an evenly-strided 1D array that
has been reshaped, you can still reshape it without a copy if it looks
like an array of such arrays, each of which can be reshaped
separately. (For example, you can view an array of length (8,2,3) that
you have to resize to an array of size (4,4,3,2) as two separate
pieces: an array of size (2,3) that you have to resize to (3,2) and an
array of size (8,) that you need to resize to size (4,4). Each of
these pieces separately needs to look like a reshaped 1d array, but
there need be no relation between them.) (The condition for when you
can pull a smaller resize operation off the end is precisely when you
can find a tail segment of each tuple that has the same product; in
the case above, 2*3=3*2, so we could stop and do that reshape
separately.)

(3) If the array cannot be reshaped without a copy using the rules
above, it cannot be reshaped without a copy at all.


The python code in my previous message actually implements this rule,
if you want a less vague description.

I should also point out that views are not always more efficient than
copies (because of cache concerns, and because a small view can
prevent a giant block of memory from being deallocated); nevertheless
it should be up to the user to copy when it helps.

A. M. Archibald

Re: [Numpy-discussion] The NumPy Fortran-ordering quiz

[A. M. A. <per...@gm...>]

On 18/10/06, Charles R Harris <cha...@gm...> wrote:
>
>
> On 10/17/06, A. M. Archibald <per...@gm...> wrote:
> > On 17/10/06, Charles R Harris <cha...@gm...> wrote:
> > >
> > >
> > > On 10/17/06, Travis Oliphant <oli...@ie... > wrote:
> >
> > > > Thus, reshape does the equivalent of a Fortran ravel to [1,4,2,5,3,6]
> > > > and then a Fortran-order based fill of an empty (3,2) array:  giving
> you
> > > > the result.
> > >
> > > Why a Fortran ravel? I am thinking of it as preserving the memory
> layout,
> > > just messing with how it is addressed.
> >
> > Because, to first order, the memory layout is totally irrelevant to a
> > python user.
> >
> > array([[1,2,3],[4,5,6]],order='C')
> > array([[1,2,3],[4,5,6]],order='F')
> >
> > should be nearly identical to the python user.
>
> So what is the point of having a fortran layout if things are not actually
> layed out in fortran order? As far as I can see, memory layout is the only
> reason for fortran ordering. Now I can sort of see where the fortran ravel
> comes from, because the result can be passed to a fortran routine. And it
> does look like a.ravel('F') is the same as a.T.ravel(). Hmmm. Now I wonder
> about this:

Thins are laid out in Fortran order if you request Fortran order upon
array creation. You just can't see it, normally. Numpy tries hard to
make arrays of all different orders look identical.

Of course, the reason I said "to first order" is that when determining
when to copy and when not to, or when passing arrays to C or Fortran
functions, sometimes it does matter. For example, a typical Fortran
linear algebra routine needs its arrays to be in Fortran order in
memory.

> In [62]: array([[1,2,3],[4,5,6]], dtype=int8, order='F').flags
>  Out[62]:
>   CONTIGUOUS : False
>    FORTRAN : True
>   OWNDATA : True
>    WRITEABLE : True
>   ALIGNED : True
>    UPDATEIFCOPY : False

Unless you're working with C extensions, it doesn't make much sense to
ever look at flags. The memory layout does not affect normal array
behaviour (including reshape).

A. M. Archibald

Re: [Numpy-discussion] questions regarding assignement and copy

[Sven S. <sve...@gm...>]

David Cournapeau schrieb:
> Sven Schreiber wrote:
>> Yes it's intended; as far as I understand the python/numpy syntax, <+>
>> is an operator, and that triggers assignment by copy (even if you do
>> something trivial as bar = +foo, you get a copy, if I'm not mistaken),
>>   
> So basically, whenever you have
> 
> foo = expr
> 
> with expr is a numpy expression containing foo, you trigger a copy ?
> 

David, I would tend to confirm your conjecture, but given its generality
  I should not pretend to really know all about that stuff.

It could be, for example, that the user may specify inside the
expression that she wants a view, not a copy, But I don't know if that's
really possible.

Anybody else available to help out?
-sven

Re: [Numpy-discussion] questions regarding assignement and copy

[David C. <da...@ar...>]

Sven Schreiber wrote:
>
> Yes it's intended; as far as I understand the python/numpy syntax, <+>
> is an operator, and that triggers assignment by copy (even if you do
> something trivial as bar = +foo, you get a copy, if I'm not mistaken),
>   
So basically, whenever you have

foo = expr

with expr is a numpy expression containing foo, you trigger a copy ?

David

[Numpy-discussion] Question about boolean types comparisons

[Francesc A. <fa...@ca...>]

Hi,

I've seen that bool type in numpy seems to work more consistently than
its python counterpart:

In [28]:True+True
Out[28]:2
In [29]:numpy.bool_(True)+numpy.bool_(True)
Out[29]:True

Good! but...

In [30]:True > False
Out[30]:True
In [31]:numpy.bool_(True) > numpy.bool_(False)
Out[31]:True

Perhaps raising an error saying something like "boolean types cannot be
compared" would be nice. Not too important, but worth to notice, IMO.

-- 
>0,0<   Francesc Altet    |  Be careful about using the following code --
V   V   Carabos Coop. V.  |  I've only proven that it works, 
 "-"    Enjoy Data        |  I haven't tested it. -- Donald Knuth

Re: [Numpy-discussion] questions regarding assignement and copy

[Sven S. <sve...@gm...>]

David Cournapeau schrieb:
> Hi there,
> 
>     I've just managed to nail down a bug which took me nearly two whole 
> days to find: this is coming from an unexpected (at least from me) 
> behaviour of numpy.

You have all my sympathy, I tripped over something similar not too long
ago, so welcome to the club.

>     Now, if I do:
> 
>     bar   += 1
>     print bar is foo
> 
>     prints True
> 
>     But if I do bar = bar + 1, then bar is not a copy of foo anymore. Is 
> this intended ? This looks really confusing to me, and I would like to 
> know what the precise rules about copy vs alias are ?
>

Yes it's intended; as far as I understand the python/numpy syntax, <+>
is an operator, and that triggers assignment by copy (even if you do
something trivial as bar = +foo, you get a copy, if I'm not mistaken),
while <+=> is syntactically not really an operator. AFAIK it's sole
purpose is really to _not_ make a copy to save memory (well, it also
looks nice ;-).

You'll need somebody else to tell you the _precise_ rules, though...

-sven

[Numpy-discussion] questions regarding assignement and copy

[David C. <da...@ar...>]

Hi there,

    I've just managed to nail down a bug which took me nearly two whole 
days to find: this is coming from an unexpected (at least from me) 
behaviour of numpy.
    I understand that if foo is a numpy array, doing bar = foo makes no 
copy, and whenever you change the value of foo, you change the value of bar:

    import numpy as N
    foo   = N.linspace(0, 4, 5)
    bar   = foo
    print foo[1]
    bar[1] = -1
    print foo[1]

    => prints 1 and -1.

    Now, if I do:

    bar   += 1
    print bar is foo

    prints True

    But if I do bar = bar + 1, then bar is not a copy of foo anymore. Is 
this intended ? This looks really confusing to me, and I would like to 
know what the precise rules about copy vs alias are ?

    Cheers,

    David

[Numpy-discussion] =?iso-2022-jp?b?GyRCPlA0aSROQUdFKCRKPXcbKEI=?= =?iso-2022-jp?b?GyRCJE47UhsoQg==?=

[<GR...@HO...>]

完全無料サイト大特集【秋の旬】
http://hlgliu.jeeran.com/egao

[Numpy-discussion] NumPy 1.0rc3 will be out in a few hours

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

I'm going to release NumPy 1.0rc3 in a few hours.   It requires a 
re-build of extension modules from 1.0rc2.  This gives us a chance to 
start fresh with NumPy 1.0 and also test the non-dying version of the 
version_checking code.  So, it's not all bad.

Official NumPy 1.0 will be released on Oct. 24 to many tears of rejoicing...

Please file all tickets / requests by this weekend if you hope to have 
them included in the 1.0 release.

-Travis

Re: [Numpy-discussion] recent ubuntu package

[Andrew S. <str...@as...>]

Sven Schreiber wrote:
> Hi,
> as suggested on the website I use the kindly provided pre-built
> (unofficial) ubuntu debs. Recently there is a new one available with
> version numbe 1.0+~dev3336-0ads1.
>
> Apart from the slightly strange +~ thing in there, it very much seems to
> be based on trac changeset 3336, which is somewhere between rc2 and now,
> right?
>   
Yes, it is from svn 3336 with one minor patch. The strange +~ is to make
the version numbers sort properly. (I didn't think about this when
releasing "1.0b2".
> And just checking, what are the matplotlib packages that are compatible?
> The most recent from the same package source seems to be
> 0.87.5-r2781-0ads2, is that ok with the above numpy? Will the mpl ubuntu
> package before that break with the latest numpy ubuntu package?
>   
The matplotlib .deb on my website is working fine for me with the latest
numpy .deb there. If there are any recent patches or anything you are
missing, please let me know -- it's not really a big deal to update
them, although it might take a couple of days for me to find the time.

Re: [Numpy-discussion] The NumPy Fortran-ordering quiz

[A. M. A. <per...@gm...>]

On 17/10/06, Travis Oliphant <oli...@ie...> wrote:
> Charles R Harris wrote:
> >
> > In [3]: array([[1,2,3,4,5,6]]).reshape((3,2),order='F')
> > Out[3]:
> > array([[1, 4],
> >        [2, 5],
> >        [3, 6]])
> >
> > I also don't understand why a copy is returned if 'F' just fiddles
> > with the indices and strides; the underlying data should be the same,
> > just the view changes. FWIW, I think both examples should be returning
> > views.
>
> You are right, it doesn't need to.   My check is not general enough.
>
> It can be challenging to come up with a general way to differentiate the
> view-vs-copy situation and I struggled with it.  In this case, it's the
> fact that while self->nd > 1, the other dimensions are only of shape 1
> and so don't really matter.   If you could come up with some kind of
> striding check that would distinguish the two cases, I would appreciate it.

I wrote, just for exercise I guess, a routine that checks to see if a
copy is needed for a reshape (and if not, it computes the resulting
strides) for an arbitrary array. (It currently only does a C-order
reshape, but F-order would be an easy addition.) It does, however,
successfully handle arbitrarily strided arrays with arbitrary
arrangements of "1" dimensions having arbitrary strides. I think it
also correctly handles 0 strides. I don't imagine you'd want to *use*
it, but it does provide a complete solution to copy-less reshaping.

I'd be happy to help with the C implementation built into numpy, but
for the life of me I can't figure out how it works.


A. M. Archibald

P.S. I haven't written down a *proof* that this implementation never
copies unnecessarily, but I'm pretty sure that it doesn't. It has a
reshape_restricted that effectively does a ravel() first; this can't
cope with things like length (5,2,3) -> length (5,3,2) unless the
strides are conveniently equal, so there's a wrapper, reshape(), that
breaks the reshape operation up into pieces that can be done by
reshape_restricted. -AMA

Re: [Numpy-discussion] The NumPy Fortran-ordering quiz

[Charles R H. <cha...@gm...>]

On 10/17/06, Charles R Harris <cha...@gm...> wrote:
>
>
>
> On 10/17/06, A. M. Archibald <per...@gm...> wrote:
> >
> > On 17/10/06, Charles R Harris <cha...@gm...> wrote:
> > >
> > >
> > > On 10/17/06, Travis Oliphant <oli...@ie... > wrote:
>
>
<snip>

Which doesn't seem to be the case here. I am beginning to wonder if we
> really need fortran order, seems that a few well chosen interface routines
> would fill the need and avoid much confusion.
>

For instance, it would be nice if flatten took an order keyword:

In [107]: array([[1,2,3],[4,5,6]], dtype=int8, order='F').flatten()
Out[107]: array([1, 2, 3, 4, 5, 6], dtype=int8)

Chuck

[Numpy-discussion] sort(axis=-1)

[Christian K. <ck...@ho...>]

Hi,

it seems that -1 as axis parameter is interpreted like in array indexing, i.e.
-1 means the last axis rather than meaning the only axis of the flattened
representation, like for example with take().
Is that intendend?

Christian

Re: [Numpy-discussion] The NumPy Fortran-ordering quiz

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

>
>         You are right, it doesn't need to.   My check is not general
>         enough.
>
>         It can be challenging to come up with a general way to
>         differentiate the
>         view-vs-copy situation and I struggled with it.  In this case,
>         it's the
>         fact that while self->nd > 1, the other dimensions are only of
>         shape 1
>         and so don't really matter.   If you could come up with some
>         kind of
>         striding check that would distinguish the two cases, I would
>         appreciate it. 
>
>
>     I suppose the problem is mostly in discontiguous arrays. Hmmm...,
>     this isn't too different that reshaping the transpose.
>
>     a.reshape((m,n),order='F' ~ a.reshape((n,m)).T.reshape(m,n)
>

I just committed a change to the check-for-copy code to SVN.   A copy 
will occur if an actual reshaping is taking place that involves more 
than just adding or deleting ones from the old shape and if

1) self is not "single-segment".  Single-segment means either 
Fortran-order or C-order.  If we don't have a single-segment array, then 
we have to get a (Fortran or C-) contiguous buffer to do any reshaping 
(there may be special cases when it's possible, but I haven't figure 
them out...)

2) self is single-segment but self.squeeze().ndim > 1 and it is in the 
"wrong" order from what was requested in the order argument (i.e. self 
is in C-contiguous order but a Fortran-order reshape was requested or 
self is in F-contiguous order but a C-order reshape was requested).

If there are any cases satisfying these rules where a copy does not have 
to occur then let me know.


-Travis

Re: [Numpy-discussion] The NumPy Fortran-ordering quiz

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

Stefan van der Walt wrote:
> One last case, which confuses me still (probably because it is
> 04:16am):
>
> In [41]: x = N.array([[0,1,2],[3,4,5]],order='F')
>
> In [42]: x
> Out[42]: 
> array([[0, 1, 2],
>        [3, 4, 5]])
>
> I assume the data is now stored in memory as
>
> [0 3 1 4 2 5] (column-wise)
>
> If I now do
>
> x.reshape((3,2),order='C')
>
> i.e. take that block of memory, assume it is in 'C' order, and make
> its shape (3,2), I expect
>
> [[0 3]
>  [1 4]
>  [2 5]]
>
> but get
>
> [[1 2]
>  [3 4]
>  [5 6]]
>
> I'm obviously missing something trivial -- I'll try again tomorrow.
>   

I think I see what is going on and where people are getting tripped up.  
  You have to remember, that to NumPy it doesn't semantically matter 
what the "ordering" of the array is.  There is no guarantee that C- 
order or Fortran-order is *ever* preserved through an operation.    
Because, in fact the general memory model of the array has no defined 
"order".  It's defined by the strides array.  It just so happens that 
two special-cases are tracked so that we can call out to compiled 
routines that expect contiguous arrays more easily.  

So, your mistake is trying to think that the "block" of memory is [0, 3, 
1, 4, 2, 5]  when you pass the Fortran-order array to the reshape 
method.    While this is true, and it means that you will save a copy if 
you passed this off to a Fortran routine, the reshape command does not 
use this information in determining how to "think-about" the input 
array.  In fact,  the reshape method does not allow any way to specify 
the order of the "input" array (self) separately from the order of the 
output array.   The order argument indicates the defined order of both 
input and output.  You might think that the order of self should be used 
as the order of  the input array.  The problem with this is, again, that 
a general array does not have a defined "order".  What should be used as 
the assumed "order" for an un-strided array?  You're left with an 
unresolved question.

To avoid two input order arguments, I just let order indicate the order 
for both the input and the output arrays.  We could provide both, but 
this seems a bit over-done as the same could be accomplished by 
separately raveling the input to 1-d and then specifying the order 
argument on reshape.

Please continue to question.  All the code needs as much review as it 
can get.

Best regards,


-Travis

Re: [Numpy-discussion] The NumPy Fortran-ordering quiz

[Charles R H. <cha...@gm...>]

On 10/17/06, A. M. Archibald <per...@gm...> wrote:
>
> On 17/10/06, Charles R Harris <cha...@gm...> wrote:
> >
> >
> > On 10/17/06, Travis Oliphant <oli...@ie...> wrote:
>
> > > Thus, reshape does the equivalent of a Fortran ravel to [1,4,2,5,3,6]
> > > and then a Fortran-order based fill of an empty (3,2) array:  giving
> you
> > > the result.
> >
> > Why a Fortran ravel? I am thinking of it as preserving the memory
> layout,
> > just messing with how it is addressed.
>
> Because, to first order, the memory layout is totally irrelevant to a
> python user.
>
> array([[1,2,3],[4,5,6]],order='C')
> array([[1,2,3],[4,5,6]],order='F')
>
> should be nearly identical to the python user.


So what is the point of having a fortran layout if things are not actually
layed out in fortran order? As far as I can see, memory layout is the only
reason for fortran ordering. Now I can sort of see where the fortran ravel
comes from, because the result can be passed to a fortran routine. And it
does look like a.ravel('F') is the same as a.T.ravel(). Hmmm. Now I wonder
about this:

In [62]: array([[1,2,3],[4,5,6]], dtype=int8, order='F').flags
Out[62]:
  CONTIGUOUS : False
  FORTRAN : True
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False

Now, either CONTIGUOUS is in error, because it really *is* fortran
contiguous (but not c contiguous), or the array cannot be passed to a
fortran routine that expects fortran ordering, or CONTIGUOUS refers to C
addressing, which we already know is not contiguous, in which case we feel
uncertain. Note that knowing the answer matters if I want to call a fortran
routine with this by pulling out the data pointer. The fortran routine could
be in a library, or maybe in the LaPack wrapper, but whatever, it hasn't
been wrapped by f2py that takes care of such details.

This also looks fishy:

In [93]: asfortranarray(array([[1,2,3],[4,5,6]], dtype=int8))
Out[93]:
array([[1, 2, 3],
       [4, 5, 6]], dtype=int8)

In [96]: asfortranarray(array([[1,2,3],[4,5,6]], dtype=int8)).flags
Out[96]:
  CONTIGUOUS : False
  FORTRAN : True
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False

because the docstring says:

asfortranarray(a, dtype=None)
    Return 'a' as an array laid out in Fortran-order in memory.


Which doesn't seem to be the case here. I am beginning to wonder if we
really need fortran order, seems that a few well chosen interface routines
would fill the need and avoid much confusion.

Chuck

Re: [Numpy-discussion] The NumPy Fortran-ordering quiz

[A. M. A. <per...@gm...>]

On 17/10/06, Charles R Harris <cha...@gm...> wrote:
>
>
> On 10/17/06, Travis Oliphant <oli...@ie...> wrote:

> > Thus, reshape does the equivalent of a Fortran ravel to [1,4,2,5,3,6]
> > and then a Fortran-order based fill of an empty (3,2) array:  giving you
> > the result.
>
> Why a Fortran ravel? I am thinking of it as preserving the memory layout,
> just messing with how it is addressed.

Because, to first order, the memory layout is totally irrelevant to a
python user.

array([[1,2,3],[4,5,6]],order='C')
array([[1,2,3],[4,5,6]],order='F')

should be nearly identical to the python user. If the storage order
mattered, you'd have to know, every time you used reshape, what order
your matrix was stored in (did it come from a transpose operation, for
example?).

As for why a Fortran ravel rather than a C ravel, that's a simplifying
decision. If you like, you can think of reshape as happening in two
steps:

o1='C'
o2='C'
A.reshape((6,),order=o1).reshape((2,3),order=o2)

The design could have allowed, as Travis Oliphant says, o1!=o2, but
explaining that would have been quite difficult (even more than now, I
mean). And in any case you can use the code above to achieve that, if
you really want.

A. M. Archibald

Re: [Numpy-discussion] The NumPy Fortran-ordering quiz

[Charles R H. <cha...@gm...>]

On 10/17/06, Charles R Harris <cha...@gm...> wrote:
>
>
>
> On 10/17/06, Travis Oliphant <oli...@ie...> wrote:
> >
> > Charles R Harris wrote:
> > >
> > >
> > > On 10/17/06, *Lisandro Dalcin* <da...@gm...
> > > <mailto:da...@gm...>> wrote:
> > >
> > >     I was surprised by this
> > >
> > >     In [14]: array([[1,2,3],[4,5,6]]).reshape((3,2),order='F')
> > >     Out[14]:
> > >     array([[1, 5],
> > >            [4, 3],
> > >            [2, 6]])
> > >
> > >
> > > This one still looks wrong.
> > >
> > >     In [15]: array([1,2,3,4,5,6]).reshape((3,2),order='F')
> > >     Out[15]:
> > >     array([[1, 2],
> > >            [3, 4],
> > >            [5, 6]])
> > >
> > >
> > >
> > >  This one is fixed,
> > >
> > > In [3]: array([[1,2,3,4,5,6]]).reshape((3,2),order='F')
> > > Out[3]:
> > > array([[1, 4],
> > >        [2, 5],
> > >        [3, 6]])
> > >
> > > I also don't understand why a copy is returned if 'F' just fiddles
> > > with the indices and strides; the underlying data should be the same,
> > > just the view changes. FWIW, I think both examples should be returning
> >
> > > views.
> >
> > You are right, it doesn't need to.   My check is not general enough.
> >
> > It can be challenging to come up with a general way to differentiate the
> > view-vs-copy situation and I struggled with it.  In this case, it's the
> > fact that while self->nd > 1, the other dimensions are only of shape 1
> > and so don't really matter.   If you could come up with some kind of
> > striding check that would distinguish the two cases, I would appreciate
> > it.
>
>
> I suppose the problem is mostly in discontiguous arrays. Hmmm..., this
> isn't too different that reshaping the transpose.
>
> a.reshape((m,n),order='F' ~ a.reshape((n,m)).T.reshape(m,n)
>

Erm, brainfart. Make that a.reshape((m,n),order='F') ~ a.reshape((n,m)).T

That's how I think of it, anyway. Actual layout in memory is not affected
until something like ascontiguosarray is called.

Chuck