numpy-discussion

[Numpy-discussion] repmat equivalent?

[Albert S. <fu...@gm...>]

Hello all

I recently started using NumPy and one function that I am really
missing from MATLAB/Octave is repmat. This function is very useful for
implementing algorithms as matrix multiplications instead of for
loops.

Here's my first attempt at repmat for 1d and 2d (with some
optimization by Stefan van der Walt):

def repmat(a, m, n):
    if a.ndim =3D=3D 1:
        a =3D array([a])
    (origrows, origcols) =3D a.shape
    rows =3D origrows * m
    cols =3D origcols * n
    b =3D a.reshape(1,a.size).repeat(m, 0).reshape(rows, origcols).repeat(n=
, 0)
    return b.reshape(rows, cols)

print repmat(array([[1,2],[3,4]]), 2, 3)

produces:

[[1 2 1 2 1 2]
 [3 4 3 4 3 4]
 [1 2 1 2 1 2]
 [3 4 3 4 3 4]]

which is the same as in MATLAB.

There are various issues with my function that I don't quite know how to so=
lve:

- How to handle scalar inputs (probably need asarray here)
- How to handle more than 2 dimensions

More than 2 dimensions is tricky, since NumPy and MATLAB don't seem to
agree on how more-dimensional data is organised? As such, I don't know
what a NumPy user would expect repmat to do with more than 2
dimensions.

Here are some test cases that the current repmat should pass, but doesn't:

a =3D repmat(1, 1, 1)
assert_equal(a, 1)
a =3D repmat(array([1]), 1, 1)
assert_array_equal(a, array([1]))
a =3D repmat(array([1,2]), 2, 3)
assert_array_equal(a, array([[1,2,1,2,1,2], [1,2,1,2,1,2]]))
a =3D repmat(array([[1,2],[3,4]]), 2, 3)
assert_array_equal(a, array([[1,2,1,2,1,2], [3,4,3,4,3,4],
[1,2,1,2,1,2], [3,4,3,4,3,4]]))

Any suggestions on how do repmat in NumPy would be appreciated.

Regards

Albert

Re: [Numpy-discussion] repmat equivalent?

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

On Thu, Feb 23, 2006 at 01:31:47PM +0200, Albert Strasheim wrote:
> More than 2 dimensions is tricky, since NumPy and MATLAB don't seem to
> agree on how more-dimensional data is organised? As such, I don't know
> what a NumPy user would expect repmat to do with more than 2
> dimensions.

To expand on this, here is what I see when I create (M,N,3) matrices
in both octave and numpy.  I expect to see an MxN matrix stacked 3
high:

octave
------

octave:1> zeros(2,2,3))
ans =3D

ans(:,:,1) =3D

  0  0
  0  0

ans(:,:,2) =3D

  0  0
  0  0

ans(:,:,3) =3D

  0  0
  0  0

numpy
-----

In [19]: zeros((2,3,3))
Out[19]:
array([[[0, 0, 0],
        [0, 0, 0],
        [0, 0, 0]],

       [[0, 0, 0],
        [0, 0, 0],
        [0, 0, 0]]])


There is nothing wrong with numpy's array -- but the output generated
seems counter-intuitive.

St=E9fan

Re: [Numpy-discussion] repmat equivalent?

[Albert S. <fu...@gm...>]

Hello all

On 2/23/06, Stefan van der Walt <st...@su...> wrote:
> On Thu, Feb 23, 2006 at 01:31:47PM +0200, Albert Strasheim wrote:
> > More than 2 dimensions is tricky, since NumPy and MATLAB don't seem to
> > agree on how more-dimensional data is organised? As such, I don't know
> > what a NumPy user would expect repmat to do with more than 2
> > dimensions.
>
> To expand on this, here is what I see when I create (M,N,3) matrices
> in both octave and numpy.  I expect to see an MxN matrix stacked 3
> high:

<snip>

There are other (unexpected, for me at least) differences between
MATLAB/Octave and NumPy too. For a 3D array in MATLAB, only indexing
on the last dimension yields a 2D array, where NumPy always returns a
2D array.

I put some examples for the 3D case at:

http://students.ee.sun.ac.za/~albert/numscipy.html

Regards

Albert

[Numpy-discussion] Re: repmat equivalent?

[Albert S. <fu...@gm...>]

Hello all

Just to clear up any confusion:

On 2/23/06, Albert Strasheim <fu...@gm...> wrote:

> Here are some test cases that the current repmat should pass, but doesn't=
:
>
> a =3D repmat(1, 1, 1)
> assert_equal(a, 1)
> a =3D repmat(array([1]), 1, 1)
> assert_array_equal(a, array([1]))
> a =3D repmat(array([1,2]), 2, 3)
> assert_array_equal(a, array([[1,2,1,2,1,2], [1,2,1,2,1,2]]))
> a =3D repmat(array([[1,2],[3,4]]), 2, 3)
> assert_array_equal(a, array([[1,2,1,2,1,2], [3,4,3,4,3,4],
> [1,2,1,2,1,2], [3,4,3,4,3,4]]))

Only the first two tests fail. The other two pass. Presumably any test
that uses a matrix with more than 2 dimensions will also fail.

Regards

Albert

Re: [Numpy-discussion] repmat equivalent?

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

Albert Strasheim wrote:
> There are other (unexpected, for me at least) differences between
> MATLAB/Octave and NumPy too.

First: numpy is not, and was never intended to be, a MATLAB clone, 
work-alike, whatever. You should *expect* there to be differences.

 > For a 3D array in MATLAB, only indexing
> on the last dimension yields a 2D array, where NumPy always returns a
> 2D array.

I think the key here is that MATLAB's core data type is a matrix, which 
is 2-d. The ability to do 3-d arrays was added later, and it looks like 
they are still preserving the core matrix concept, so that a 3-d array 
is not really a 3-d array; it is, as someone on this thread mentioned, a 
"stack" of matrices.

In numpy, the core data type is an n-d array. That means that there is 
nothing special about 2-d vs 4-d vs whatever, except 0-d (scalars). So a 
3-d array is a cube shape, that you might want to pull a 2-d array out 
of it in any orientation. There's nothing special about which axis 
you're indexing. For that reason, it's very important that indexing any 
axis will give you the same rank array.

Here's the rule:

-- indexing reduces the rank by 1, regardless of which axis is being 
indexed.
 >>> import numpy as N
 >>> a = N.zeros((2,3,4))
 >>> a
array([[[0, 0, 0, 0],
         [0, 0, 0, 0],
         [0, 0, 0, 0]],

        [[0, 0, 0, 0],
         [0, 0, 0, 0],
         [0, 0, 0, 0]]])
 >>> a[0,:,:]
array([[0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0]])
 >>> a[:,0,:]
array([[0, 0, 0, 0],
        [0, 0, 0, 0]])
 >>> a[:,:,0]
array([[0, 0, 0],
        [0, 0, 0]])

-- slicing does not reduce the rank:
 >>> a[:,0:1,:]
array([[[0, 0, 0, 0]],

        [[0, 0, 0, 0]]])
 >>> a[:,0:1,:].shape
(2, 1, 4)

It's actually very clean, logical, and useful.

-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] repmat equivalent?

[Albert S. <fu...@gm...>]

Hello all

On 2/23/06, Christopher Barker <Chr...@no...> wrote:
> Albert Strasheim wrote:
> > There are other (unexpected, for me at least) differences between
> > MATLAB/Octave and NumPy too.
>
> First: numpy is not, and was never intended to be, a MATLAB clone,
> work-alike, whatever. You should *expect* there to be differences.

I understand this. As a new user, I'm trying to understand these difference=
s.

>  > For a 3D array in MATLAB, only indexing
> > on the last dimension yields a 2D array, where NumPy always returns a
> > 2D array.
>
> I think the key here is that MATLAB's core data type is a matrix, which
> is 2-d. The ability to do 3-d arrays was added later, and it looks like
> they are still preserving the core matrix concept, so that a 3-d array
> is not really a 3-d array; it is, as someone on this thread mentioned, a
> "stack" of matrices.
>
> In numpy, the core data type is an n-d array. That means that there is
> nothing special about 2-d vs 4-d vs whatever, except 0-d (scalars). So a
> 3-d array is a cube shape, that you might want to pull a 2-d array out
> of it in any orientation. There's nothing special about which axis
> you're indexing. For that reason, it's very important that indexing any
> axis will give you the same rank array.
>
> Here's the rule:
>
> -- indexing reduces the rank by 1, regardless of which axis is being
> indexed.

<snip>

Thanks for your comments. These cleared up a few questions I had about
NumPy's design. However, I'm still wondering how the average NumPy
user would expect repmat implemented for NumPy to behave with arrays
with more than 2 dimensions.

I would like to clear this up, since I think that a good repmat
function is an essential tool for implementing algorithms that use
matrix multiplication instead of for loops to perform operations
(hopefully with a significant speed increase). If there is another way
of accomplishing this, I would love to know.

Regards

Albert

Re: [Numpy-discussion] repmat equivalent?

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

Albert Strasheim wrote:

>Hello all
>
>I recently started using NumPy and one function that I am really
>missing from MATLAB/Octave is repmat. This function is very useful for
>implementing algorithms as matrix multiplications instead of for
>loops.
>  
>
There is a function in scipy.linalg called kron that could be brought 
over which can do a repmat.


In file: /usr/lib/python2.4/site-packages/scipy/linalg/basic.py

def kron(a,b):
    """kronecker product of a and b

    Kronecker product of two matrices is block matrix
    [[ a[ 0 ,0]*b, a[ 0 ,1]*b, ... , a[ 0 ,n-1]*b  ],
     [ ...                                   ...   ],
     [ a[m-1,0]*b, a[m-1,1]*b, ... , a[m-1,n-1]*b  ]]
    """
    if not a.flags['CONTIGUOUS']:
        a = reshape(a, a.shape)
    if not b.flags['CONTIGUOUS']:
        b = reshape(b, b.shape)
    o = outerproduct(a,b)
    o=o.reshape(a.shape + b.shape)
    return concatenate(concatenate(o, axis=1), axis=1)


Thus,

kron(ones((2,3)), arr)

 >>> sl.kron(ones((2,3)),arr)
array([[1, 2, 1, 2, 1, 2],
       [3, 4, 3, 4, 3, 4],
       [1, 2, 1, 2, 1, 2],
       [3, 4, 3, 4, 3, 4]])

gives you the equivalent of

repmat(arr, 2,3)

We could bring this over from scipy into numpy as it is simple enough.

It has a multidimensional extension (i.e. you can pass in a and b as 
higher dimensional arrays),  But, don't ask me to explain it to you 
because I can't without further study....

-Travis

Re: [Numpy-discussion] repmat equivalent?

[Albert S. <fu...@gm...>]

Hello all

On 2/23/06, Travis Oliphant <oli...@ie...> wrote:
> Albert Strasheim wrote:
>
> >Hello all
> >
> >I recently started using NumPy and one function that I am really
> >missing from MATLAB/Octave is repmat. This function is very useful for
> >implementing algorithms as matrix multiplications instead of for
> >loops.
> >
> >
> There is a function in scipy.linalg called kron that could be brought
> over which can do a repmat.

<snip>

> Thus,
>
> kron(ones((2,3)), arr)
>
>  >>> sl.kron(ones((2,3)),arr)
> array([[1, 2, 1, 2, 1, 2],
>        [3, 4, 3, 4, 3, 4],
>        [1, 2, 1, 2, 1, 2],
>        [3, 4, 3, 4, 3, 4]])
>
> gives you the equivalent of
>
> repmat(arr, 2,3)

Thanks! Merging this into numpy would be much appreciated. Stefan van
der Walt did some benchmarks and this approach seems faster than
anything we managed for 2D arrays.

However, I'm a bit concerned about the ones(n,m) that is needed by
this implementation. It seems to me that this would double the memory
requirements of this repmat function, which is fine when working with
small matrices, but could be a problem with larger ones.

Any thoughts?

Regards

Albert

Re: [Numpy-discussion] repmat equivalent?

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

On Thu, Feb 23, 2006 at 11:21:39PM +0200, Albert Strasheim wrote:
> > Thus,
> >
> > kron(ones((2,3)), arr)
> >
> >  >>> sl.kron(ones((2,3)),arr)
> > array([[1, 2, 1, 2, 1, 2],
> >        [3, 4, 3, 4, 3, 4],
> >        [1, 2, 1, 2, 1, 2],
> >        [3, 4, 3, 4, 3, 4]])
> >
> > gives you the equivalent of
> >
> > repmat(arr, 2,3)
>=20
> Thanks! Merging this into numpy would be much appreciated. Stefan van
> der Walt did some benchmarks and this approach seems faster than
> anything we managed for 2D arrays.

My benchmark was wrong -- this function is not as fast as the version
Albert previously proposed.  Below follows the benchmark of seven
possible repmat functions:

-------------------------------------------------------------------------=
--
0 : 1.09316706657 (Albert)
1 : 6.15612506866 (Stefan)
2 : 5.21671295166 (Stefan)
3 : 2.78160500526 (Stefan)
4 : 1.20426011086 (Albert Optimised)
5 : 11.0923781395 (Travis)
6 : 3.47499799728 (Alex)
-------------------------------------------------------------------------=
--
0 : 1.17543005943
1 : 6.03165698051
2 : 5.7597899437
3 : 2.40381717682
4 : 1.09497308731
5 : 11.6657807827
6 : 7.11567497253
-------------------------------------------------------------------------=
--
0 : 2.03999996185
1 : 9.87535595894
2 : 8.86893296242
3 : 4.56993699074
4 : 2.02298903465
5 : 22.8858327866
6 : 10.7882151604
-------------------------------------------------------------------------=
--

I attach the code.

St=E9fan

Re[2]: [Numpy-discussion] repmat equivalent?

[Alan G I. <ai...@am...>]

On Fri, 24 Feb 2006, Stefan van der Walt apparently wrote: 
> Below follows the benchmark of seven possible repmat 
> functions: 

Probably I a misunderstanding something here.  But I thought 
the idea of repmat was that it used a single copy of the 
data to represent multiple copies in a matrix, and all these 
functions seem ultimately to use multiple copies of the 
data.  If that is right, then a repmat should be a
subclass of matrix. I think ...

Cheers,
Alan Isaac

Re: [Numpy-discussion] repmat equivalent?

[Albert S. <fu...@gm...>]

Hello

On 2/23/06, Travis Oliphant <oli...@ie...> wrote:
> Albert Strasheim wrote:
>
> >Hello all
> >
> >I recently started using NumPy and one function that I am really
> >missing from MATLAB/Octave is repmat. This function is very useful for
> >implementing algorithms as matrix multiplications instead of for
> >loops.
> >
> >
> There is a function in scipy.linalg called kron that could be brought
> over which can do a repmat.

I quickly tried a few of my test cases with following implementation of rep=
mat:

from numpy import asarray
from scipy.linalg import kron

def repmat(a, m, n):
   a =3D asarray(a)
   return kron(ones((m, n)), a)

This test:

a =3D repmat(1, 1, 1)
assert_equal(a, 1)

fails with:

ValueError: 0-d arrays can't be concatenated

and this test:

a =3D repmat(array([1,2]), 2, 3)
assert_array_equal(a, array([[1,2,1,2,1,2], [1,2,1,2,1,2]]))

fails with:

AssertionError: Arrays are not equal (shapes (12,), (2, 6) mismatch)

Regards

Albert

Re: [Numpy-discussion] repmat equivalent?

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

Travis Oliphant schrieb:

> There is a function in scipy.linalg called kron that could be brought
> over which can do a repmat.

> 
> def kron(a,b):
>    """kronecker product of a and b

That would be great, I have been missing it in numpy already! (Because
scipy is rather big, I'd like to avoid depending on it for such things.)

So please do bring it over.
Thanks,
Sven