[SciPy-user] repmat

[Alan G Isaac]

On Tue, 06 Dec 2005, Aarre Laakso apparently wrote: 
> Can anyone give me a clue about how I might implement a fast workalike 
> for MATLAB's repmat function? 

> http://www.mathworks.com/access/helpdesk/help/techdoc/ref/repmat.html 

> The simple use case B = repmat(A,m,n) where A is an array and m,n are 
> positive integers would be sufficient for my needs. 

> I have pasted below my own attempt, which obviously involves way too 
> much copying, because it takes a couple of minutes to replicate a 
> 50-element row vector roughly 5000 times. It seems like repeat with some 
> tricky slicing might work, but I can't figure out how to do it. Or maybe 
> there's a better way? 

> Thanks for any help you can offer. 

> def repmat(A, *args): 
>     """ Replicate and tile an array. """ 

>     result = A 
>     new_size = args 
>     for axis in range(0, len(new_size)): 
>         copy_value = result 
>         num_copies = new_size[axis] 
>         for curr_copy in range(1, num_copies): 
>             result = concatenate((result, copy_value), axis) 


What you really want is not a matrix like this---much too 
costly!---but rather an instance of a repmat class (which
would keep only one copy of the data)!

In the meantime, you can try a Kronecker product for your 
purpose. (Below.)

Cheers,
Alan Isaac



def kroneckerproduct(a,b):
    '''Compute a otimes b where otimes is the Kronecker product operator.

    
    Note: the Kronecker product is also known as the matrix direct product
    or tensor product.  It is defined as follows for 2D arrays a and b
    where shape(a)=(m,n) and shape(b)=(p,q):
    c = a otimes b  => cij = a[i,j]*b  where cij is the ij-th submatrix of c.
    So shape(c)=(m*p,n*q).

    :Parameters:
     - `a`: 2-D array
     - `b`: 2-D array
    :see: http://en.wikipedia.org/wiki/Kronecker_product
    :see: http://www.scipy.org/mailinglists/mailman?fn=scipy-user/2003-September/002138.html
    :author: Nils Wagner
    :author: Alan G. Isaac
    :date: 2004-08-12
    :copyright: public domain
    :since: 2003-09-03
    :note: Contributed to numarray in Sept 2003.

    >>> print kroneckerproduct([[1,2]],[[3],[4]])
    [[3 6]
     [4 8]]
    >>> print kroneckerproduct([[1,2]],[[3,4]])
    [ [3 4 6 8]]
    >>> print kroneckerproduct([[1],[2]],[[3],[4]])
    [[3]
     [4]
     [6]
     [8]]
    '''
    a, b = asarray(a), asarray(b)
    if not (len(shape(a))==2 and len(shape(b))==2):
            raise ValueError, 'Input must be 2D arrays.'
    if not a.iscontiguous():
        a = reshape(a, a.shape)
    if not b.iscontiguous():
        b = reshape(b, b.shape)
    o = outerproduct(a,b)
    o.shape = a.shape + b.shape
    return concatenate(concatenate(o, axis=1), axis=1)