Here's a new version of 'reparray'.
I tested the previous version against 'repmat' for the 2d case and
found it to be as much as 3x slower in some instances. Ick. So I
redid it using the approach repmat uses -- reshape() and repeat()
rather than concatenate(). Now it's very nearly as fast as repmat for
the 2-d case.
def reparray(A, tup):
if numpy.isscalar(tup):
tup = (tup,)
d = len(tup)
c = numpy.array(A,copy=False,subok=True,ndmin=d)
shape = list(c.shape)
n = c.size
for i, nrep in enumerate(tup):
if nrep!=1:
c = c.reshape(-1,n).repeat(nrep,0)
dim_in = shape[i]
dim_out = dim_in*nrep
shape[i] = dim_out
n /= dim_in
return c.reshape(shape)
The full file with tests and timing code is attached.
One thing I noticed while doing this was that repmat doesn't preserve
subclasses. Which is funny considering 'mat' is part of the name and
it only works on 2-d arrays. Using asanyarray() instead of asarray()
in the first line should fix that.
--Bill
import numpy
def reparray(A, tup):
"""Repeat an array the number of times given in the integer tuple, tup.
Similar to repmat, but works for arrays of any dimension.
reparray(A,(m,n)) is equivalent to repmat(A,m,n)
If tup has length d, the result will have dimension of max(d, A.ndim).
If tup is scalar it is treated as a 1-tuple.
If A.ndim < d, A is promoted to be d-dimensional by prepending new axes.
So a shape (3,) array is promoted to (1,3) for 2-D replication,
or shape (1,1,3) for 3-D replication.
If this is not the desired behavior, promote A to d-dimensions manually
before calling this function.
If d < A.ndim, tup is promoted to A.ndim by appending 1's to it. Thus
for an A.shape of (2,3,4,5), a tup of (2,2) is treated as (2,2,1,1)
Examples:
>>> a = array([0,1,2])
>>> reparray(a,2)
array([0, 1, 2, 0, 1, 2])
>>> reparray(a,(1,2))
array([[0, 1, 2, 0, 1, 2]])
>>> reparray(a,(2,2))
array([[0, 1, 2, 0, 1, 2],
[0, 1, 2, 0, 1, 2]])
>>> reparray(a,(2,1,2))
array([[[0, 1, 2, 0, 1, 2]],
[[0, 1, 2, 0, 1, 2]]])
See Also:
repmat, repeat
"""
if numpy.isscalar(tup):
tup = (tup,)
d = len(tup)
c = numpy.array(A,copy=False,subok=True,ndmin=d)
shape = list(c.shape)
n = c.size
for i, nrep in enumerate(tup):
if nrep!=1:
c = c.reshape(-1,n).repeat(nrep,0)
dim_in = shape[i]
dim_out = dim_in*nrep
shape[i] = dim_out
n /= dim_in
return c.reshape(shape)
def _reparray_alt(A, tup):
"""Alternate implementation of reparray, that turns out to be slower.
"""
if numpy.isscalar(tup):
tup = (tup,)
d = len(tup)
A = numpy.array(A,copy=False,subok=True,ndmin=d)
for i,k in enumerate(tup):
A = numpy.concatenate([A]*k, axis=i)
return A
def _test():
a = numpy.array([0,1,2])
b = numpy.array([[0,1,2],[3,4,5]])
c = numpy.array([[[0,1,2],[3,4,5]]])
d = numpy.array([[[0,1,2]],[[3,4,5]]])
shapes = [(1,1),(2,1),(1,2),(2,2),(2,3),(3,2)]
# double check results against repmat
for s in shapes:
r = reparray(b,s)
r0 = numpy.repmat(b,s[0],s[1])
assert( numpy.all(r == r0) )
shapes += [(1,),(2,),(1,1,1),(2,1,1),(1,2,1),(1,1,2),(1,2,2),(2,1,2),(2,2,1),(2,2,2)]
# cross check results between implementations
for arr in (a,b,c,d):
for s in shapes:
r = reparray(arr,s)
r0 = _reparray_alt(arr,s)
assert( numpy.all(r == r0) )
def speed_compare():
from timeit import Timer
print '10x10 replication of 2x2'
N = 10000
t0 = Timer("c = repmat(a,10,10)","import numpy; repmat=numpy.repmat; a=numpy.ones((2,2))")
print 'repmat: \t', t0.timeit(number=N)
t1 = Timer("c = _reparray_alt(a,(10,10))","import numpy; from __main__ import _reparray_alt; a=numpy.ones((2,2))")
print 'reparray v0:\t', t1.timeit(number=N)
t2 = Timer("c = reparray(a,(10,10))","import numpy; from __main__ import reparray; a=numpy.ones((2,2))")
print 'reparray v1:\t', t2.timeit(number=N)
print '2x2 replication of 100x100'
N = 1000
t0 = Timer("c = repmat(a,2,2)","import numpy; repmat=numpy.repmat; a=numpy.ones((100,100))")
print 'repmat: \t', t0.timeit(number=N)
t1 = Timer("c = _reparray_alt(a,(2,2))","import numpy; from __main__ import _reparray_alt; a=numpy.ones((100,100))")
print 'reparray v0:\t', t1.timeit(number=N)
t2 = Timer("c = reparray(a,(2,2))","import numpy; from __main__ import reparray; a=numpy.ones((100,100))")
print 'reparray v1:\t', t2.timeit(number=N)
print '10x10 replication of 50x50'
N = 100
t0 = Timer("c = repmat(a,10,10)","import numpy; repmat=numpy.repmat; a=numpy.ones((50,50))")
print 'repmat: \t', t0.timeit(number=N)
t1 = Timer("c = _reparray_alt(a,(10,10))","import numpy; from __main__ import _reparray_alt; a=numpy.ones((50,50))")
print 'reparray v0:\t', t1.timeit(number=N)
t2 = Timer("c = reparray(a,(10,10))","import numpy; from __main__ import reparray; a=numpy.ones((50,50))")
print 'reparray v1:\t', t2.timeit(number=N)
print '10x1 replication of 100x100'
N = 100
t0 = Timer("c = repmat(a,10,1)","import numpy; repmat=numpy.repmat; a=numpy.ones((100,100))")
print 'repmat: \t', t0.timeit(number=N)
t1 = Timer("c = _reparray_alt(a,(10,1))","import numpy; from __main__ import _reparray_alt; a=numpy.ones((100,100))")
print 'reparray v0:\t', t1.timeit(number=N)
t2 = Timer("c = reparray(a,(10,1))","import numpy; from __main__ import reparray; a=numpy.ones((100,100))")
print 'reparray v1:\t', t2.timeit(number=N)
print '1x10 replication of 100x100'
N = 100
t0 = Timer("c = repmat(a,1,10)","import numpy; repmat=numpy.repmat; a=numpy.ones((100,100))")
print 'repmat: \t', t0.timeit(number=N)
t1 = Timer("c = _reparray_alt(a,(1,10))","import numpy; from __main__ import _reparray_alt; a=numpy.ones((100,100))")
print 'reparray v0:\t', t1.timeit(number=N)
t2 = Timer("c = reparray(a,(1,10))","import numpy; from __main__ import reparray; a=numpy.ones((100,100))")
print 'reparray v1:\t', t2.timeit(number=N)
if __name__ == '__main__':
_test()
speed_compare()-------------------------------------------------------------------------
Take Surveys. Earn Cash. Influence the Future of IT
Join SourceForge.net's Techsay panel and you'll get the chance to share your
opinions on IT & business topics through brief surveys -- and earn cash
http://www.techsay.com/default.php?page=join.php&p=sourceforge&CID=DEVDEV
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion
