On 17/10/06, Travis Oliphant <[EMAIL PROTECTED]> 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
from numpy import multiply, array, empty, zeros, add, asarray
from itertools import izip

class MustCopy(Exception):
    pass

def index(indices, strides):
    indices = asarray(indices)
    strides = asarray(strides)
    if len(indices)!=len(strides):
	raise ValueError
    return add.reduce(indices*strides)

def reshape_restricted(from_strides,from_lengths,to_lengths):

    if multiply.reduce(from_lengths)!=multiply.reduce(to_lengths):
	raise ValueError, "Array sizes are not compatible"
    for i in xrange(1,len(from_strides)):
	if not (from_strides[i-1]==from_strides[i]*from_lengths[i] or 
		from_lengths[i]==1): # If a dimension has length one, who cares what its stride is?
	    raise MustCopy

    r = empty(len(to_lengths), dtype=int)
    j = len(from_lengths)-1
    while j>=0 and from_lengths[j]==1:
	j-=1 # Ignore strides on length-0 dimensions
    if j<0: # Both arrays have size 1
	r[-1] = 1
    else:
	r[-1] = from_strides[j]
    i = -1
    while len(to_lengths)+i>0:
	r[i-1]=to_lengths[i]*r[i]
	i -= 1

    return r

def reshape(from_strides,from_lengths,to_lengths):
    if multiply.reduce(from_lengths)!=multiply.reduce(to_lengths):
	raise ValueError
    
    r = empty(len(to_lengths), dtype=int)

    j1 = len(from_lengths)
    j2 = len(to_lengths)

    while j1!=0:
	i1 = j1-1
	i2 = j2-1

	p1 = from_lengths[i1]
	p2 = to_lengths[i2]

	while p1!=p2:
	    if p1<p2:
		i1-=1
		assert i1>=0
		p1*=from_lengths[i1]
	    else:
		i2-=1
		assert i2>=0
		p2*=to_lengths[i2]

	r[i2:j2] = reshape_restricted(from_strides[i1:j1],
		                      from_lengths[i1:j1],
				      to_lengths[i2:j2])
	
	j1 = i1
	j2 = i2

    return r

def test_reshape(from_strides,from_lengths,to_strides,to_lengths):
    if multiply.reduce(from_lengths)!=multiply.reduce(to_lengths):
	raise ValueError
    
    if len(from_strides)!=len(from_lengths) or len(to_strides)!=len(to_lengths):
	raise ValueError

    def walk(lengths):
	i = zeros(len(lengths),int)
	while True:
	    yield i
	    j = len(lengths)-1
	    while i[j]==lengths[j]-1:
		i[j]=0
		j-=1
		if j<0:
		    return
	    i[j] += 1	

    for (i1,i2) in izip(walk(from_lengths), walk(to_lengths)):
	if index(i1,from_strides)!=index(i2,to_strides):
	    raise ValueError
-------------------------------------------------------------------------
Using Tomcat but need to do more? Need to support web services, security?
Get stuff done quickly with pre-integrated technology to make your job easier
Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo
http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion

Reply via email to