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

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

[Numpy-discussion] inconsistent use of axis= keyword argument?

[Vidar G. <vid...@37...>]

(i've been updating the cross reference of MATLAB synonymous
commands in Numeric Python to NumPy. I've kept Numeric/numarray
alternatives in the source XML, but omitted it in the PDF outputs.
see, http://37mm.no/download/matlab-python-xref.pdf.
feedback is highly appreciated.)


as i was working on this, i started wondering why 

a.max(0), a.min(0), a.ptp(0), a.flatten(0), ...

does not allow the axis=0 keyword argument used with
the exact same meaning for:

m.mean(axis=0), m.sum(axis=0), ...



and i also wonder why concatenate can't be used to stack 1-d
arrays on top of each other, returning a 2-d array?
axis relates to the number of axes in the original array(s)?

n [3]: v = arange(9)

In [7]: concatenate((v,v))
Out[7]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 8])

In [8]: concatenate((v,v),axis=0)
Out[8]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 8])

In [15]: concatenate((v,v)).reshape(2,-1)
Out[15]: 
array([[0, 1, 2, 3, 4, 5, 6, 7, 8],
       [0, 1, 2, 3, 4, 5, 6, 7, 8]])


In [5]: m = v.reshape(3,-1)

In [10]: concatenate((m,m))
Out[10]: 
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8],
       [0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])

In [11]: concatenate((m,m), axis=0)
Out[11]: 
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8],
       [0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])

In [12]: concatenate((m,m), axis=1)
Out[12]: 
array([[0, 1, 2, 0, 1, 2],
       [3, 4, 5, 3, 4, 5],
       [6, 7, 8, 6, 7, 8]])

[Numpy-discussion] A question

[Colin J. W. <cj...@sy...>]

 Bar is a sub-class of ArrayType (ndarray) and bar is an instance of Bar
[Dbg]>>> self.bar
Bar([ 5,  1,  1, 11,  3,  7, 14,  0,  5,  2,  4, 12,  9, 10,  4, 12], 
dtype=uint16)
[Dbg]>>> z= self.bar
[Dbg]>>> z
1                                                                         
<< Is this expected?
[Dbg]>>> type(self.bar)
<class '__main__.Bar'>
[Dbg]>>> self.bar.__class__.base
<attribute 'base' of 'numpy.bigndarray' objects>
[Dbg]>>>

Colin W.

Re: [Numpy-discussion] Success with SVN on Cygwin

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

Jose Gomez-Dans wrote:

>Hi!
>A few days ago, I asked about compiling NumPy on Cygwin. Travis carried out some
>modifications, and with last night's SVN, I can happily report that it now
>compiles and works. The tests produced no errors, so it's all good :)
>  
>
That's good news. I wish our unit test coverage was wide enough that 
this actually meant that all is good :-)  But, it's a good start.

Thanks are really do the cygwin ports people who already had a patch 
(though they didn't let us know about it --- I found it using Google).   
I just incorporated their patch into the build tree.

I'm glad to know that it works.

-Travis

Re: [Numpy-discussion] algorithm, optimization, or other problem?

[Brian B. <bb...@br...>]

Nadav Horesh wrote:
> It is slower.
> 
> I did a little study on this issue since I got into the issue of 
> algorithms that can not be easily vectorized (like this one).
> On my PC an outer loop step took initially 17.3 seconds, and some 
> optimization brought it down to ~11 seconds.  The dot product consumed 
> about 1/3 of the time. I estimate that objects creation/destruction 
> consumes most of the cpu time. It seems that this way comes nowhere near 
> cmex speed. I suspect that maybe blitz/boost may bridge the gap.
> 

yeah, I realized that pure python would be too slow, because I ran into the exact 
same problem with matlab scripting.  these time-dependent loops are really a mess 
when it comes to speed optimization.

After posting to the Pyrex list, someone pointed out that my loop variables had not 
been declared as c datatypes.  so, I had loops like:

for it from 0 <= it <= 1000:
	for i from 0 <= i <= 100:

		(stuff)

and the "it" and "i" were being treated, due to my oversight, as python variables. 
for speed, you need to have all the variables in the loop as c datatypes.  just 
putting a line in like:

cdef int it,i

increases the speed from 8 seconds per block to 0.2 seconds per block, which is 
comparable to the mex.

I learned that I have to be a bit more careful!  :)


			thanks,

				Brian Blais


-- 
-----------------

             bb...@br...
             http://web.bryant.edu/~bblais

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] A case for rank-0 arrays

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

Hi Sasha,

A Dissabte 18 Febrer 2006 20:49, Sasha va escriure:
> I have reviewed mailing list discussions of rank-0 arrays vs. scalars
> and I concluded that the current implementation that contains both is
> (almost) correct.  I will address the "almost" part with a concrete
> proposal at the end of this post (search for PROPOSALS if you are only
> interested in the practical part).
>
> The main criticism of supporting both scalars and rank-0 arrays is
> that it is "unpythonic" in the sense that it provides two almost
> equivalent ways to achieve the same result.  However, I am now
> convinced that this is the case where practicality beats purity.

It's a bit late, but I want to support your proposal (most of it).
I've also come to the conclusion that scalars and rank-0 arrays should
coexist. This is something that appears as a natural fact when you
have to deal regularly with general algorithms for treat objects with
different shapes. And I think you have put this very well.

Thanks,

=2D-=20
>0,0<   Francesc Altet =A0 =A0 http://www.carabos.com/
V   V   C=E1rabos Coop. V. =A0=A0Enjoy Data
 "-"

Re: [Numpy-discussion] optimizing power() for complex and real cases

[Tim H. <tim...@co...>]

I had some free time some morning, so I merged the nc_pow optimizations 
for integral powers in and committed them to the power_optimization 
branch. They could probably use more testing, but I thought someone 
might like to take a look while I'm out of town.

Also, if your looking for a way to do powers as succusive multiplies for 
fast_power or scalar_power or whatnot, starting with the algorithm for 
nc_pow would probably be a good place to start.

Hmm. Looking at this now, I realize I'm shadowing the input complex 
number 'a', with a local 'a'. Too much mindless copying from 
complexobject. That should be fixed, but I can't do it right now.

Enjoy,

-tim

Re: [Numpy-discussion] algorithm, optimization, or other problem?

[Nadav H. <Na...@Vi...>]

<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
  <meta content="text/html;charset=ISO-8859-1" http-equiv="Content-Type">
  <title></title>
</head>
<body bgcolor="#ffffff" text="#000000">
It is slower.<br>
<br>
I did a little study on this issue since I got into the issue of
algorithms that can not be easily vectorized (like this one).<br>
On my PC an outer loop step took initially 17.3 seconds, and some
optimization brought it down to ~11 seconds.&nbsp; The dot product consumed
about 1/3 of the time. I estimate that objects creation/destruction
consumes most of the cpu time. It seems that this way comes nowhere
near cmex speed. I suspect that maybe blitz/boost may bridge the gap.<br>
<br>
My optimization (did not verify numeric results)<br>
<br>
from Numeric import *<br>
from RandomArray import *<br>
import time<br>
<br>
<br>
def run():<br>
<br>
&nbsp;&nbsp;&nbsp; x=random((100,1000))&nbsp;&nbsp;&nbsp; # 1000 input vectors<br>
<br>
&nbsp;&nbsp;&nbsp; numpats=x.shape[0]<br>
&nbsp;&nbsp;&nbsp; w=random(numpats)<br>
&nbsp;&nbsp;&nbsp; th=random((1,1))<br>
&nbsp;&nbsp;&nbsp; eta=array(0.001)<br>
&nbsp;&nbsp;&nbsp; tau = array(0.01)<br>
&nbsp;&nbsp;&nbsp; old_mx=0;<br>
&nbsp;&nbsp;&nbsp; for e in range(1):<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; rnd=randint(0,numpats,250000)<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; t1=time.time()<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if 1:&nbsp; # straight python<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; for i in xrange(len(rnd)):<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; pat=rnd[i]<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; xx = x[:,pat]<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; y = array(dot(xx,w))<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; w += eta*y*(xx-y*w)<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; th += tau * (y*y-th)<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; else: # pyrex<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; dohebb(params,w,th,x,rnd)<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; print time.time()-t1<br>
<br>
if __name__ == '__main__':<br>
&nbsp;&nbsp;&nbsp; run()<br>
<br>
<br>
<br>
Bruce Southey wrote:
<blockquote
 cite="mid...@ma..."
 type="cite">
  <pre wrap="">Hi,
Actually it makes it slightly worse - given the responses on another
thread it is probably due to not pushing enough into C code.

Obviously use of blas etc will be faster but it doesn't change the
fact that removing the inner loop would be faster still.

Bruce

On 2/22/06, Nadav Horesh <a class="moz-txt-link-rfc2396E" href="mailto:na...@vi...">&lt;na...@vi...&gt;</a> wrote:
  </pre>
  <blockquote type="cite">
    <pre wrap="">You may get a significant boost by replacing the line:
  w=w+ eta * (y*x - y**2*w)
with
  w *= 1.0 - eta*y*y
  w += eta*y*x

I ran a test on a similar expression and got 5 fold speed increase.
The dot() function runs faster if you compile with dotblas.

  Nadav.


-----Original Message-----
From:   <a class="moz-txt-link-abbreviated" href="mailto:num...@li...">num...@li...</a> on behalf of Bruce Southey
Sent:   Tue 21-Feb-06 17:15
To:     Brian Blais
Cc:     <a class="moz-txt-link-abbreviated" href="mailto:pyt...@py...">pyt...@py...</a>; <a class="moz-txt-link-abbreviated" href="mailto:num...@li...">num...@li...</a>; <a class="moz-txt-link-abbreviated" href="mailto:sci...@sc...">sci...@sc...</a>
Subject:        Re: [Numpy-discussion] algorithm, optimization, or other problem?
Hi,
In the current version, note that Y is scalar so replace the squaring
(Y**2) with Y*Y as you do in the dohebb function.  On my system
without blas etc removing the squaring removes a few seconds (16.28 to
12.4). It did not seem to help factorizing Y.

Also, eta and tau are constants so define them only once as scalars
outside the loops and do the division outside the loop. It only saves
about 0.2 seconds but these add up.

The inner loop probably can be vectorized because it is just vector
operations on a matrix. You are just computing over the ith dimension
of X.  I think that you could be able to find the matrix version on
the net.

Regards
Bruce



On 2/21/06, Brian Blais <a class="moz-txt-link-rfc2396E" href="mailto:bb...@br...">&lt;bb...@br...&gt;</a> wrote:
    </pre>
    <blockquote type="cite">
      <pre wrap="">Hello,

I am trying to translate some Matlab/mex code to Python, for doing neural
simulations.  This application is definitely computing-time limited, and I need to
optimize at least one inner loop of the code, or perhaps even rethink the algorithm.
   The procedure is very simple, after initializing any variables:

1) select a random input vector, which I will call "x".  right now I have it as an
array, and I choose columns from that array randomly.  in other cases, I may need to
take an image, select a patch, and then make that a column vector.

2) calculate an output value, which is the dot product of the "x" and a weight
vector, "w", so

        y=dot(x,w)

3) modify the weight vector based on a matrix equation, like:

        w=w+ eta * (y*x - y**2*w)
               ^
               |
               +---- learning rate constant

4) repeat steps 1-3 many times

I've organized it like:

for e in 100:   # outer loop
     for i in 1000:  # inner loop
         (steps 1-3)

     display things.

so that the bulk of the computation is in the inner loop, and is amenable to
converting to a faster language.  This is my issue:

straight python, in the example posted below for 250000 inner-loop steps, takes 20
seconds for each outer-loop step.  I tried Pyrex, which should work very fast on such
a problem, takes about 8.5 seconds per outer-loop step.  The same code as a C-mex
file in matlab takes 1.5 seconds per outer-loop step.

Given the huge difference between the Pyrex and the Mex, I feel that there is
something I am doing wrong, because the C-code for both should run comparably.
Perhaps the approach is wrong?  I'm willing to take any suggestions!  I don't mind
coding some in C, but the Python API seemed a bit challenging to me.

One note: I am using the Numeric package, not numpy, only because I want to be able
to use the Enthought version for Windows.  I develop on Linux, and haven't had a
chance to see if I can compile numpy using the Enthought Python for Windows.

If there is anything else anyone needs to know, I'll post it.  I put the main script,
and a dohebb.pyx code below.


                        thanks!

                                        Brian Blais

--
-----------------

             <a class="moz-txt-link-abbreviated" href="mailto:bb...@br...">bb...@br...</a>
             <a class="moz-txt-link-freetext" href="http://web.bryant.edu/~bblais">http://web.bryant.edu/~bblais</a>




# Main script:

from dohebb import *
import pylab as p
from Numeric import *
from RandomArray import *
import time

x=random((100,1000))    # 1000 input vectors

numpats=x.shape[0]
w=random((numpats,1));

th=random((1,1))

params={}
params['eta']=0.001;
params['tau']=100.0;
old_mx=0;
for e in range(100):

     rnd=randint(0,numpats,250000)
     t1=time.time()
     if 0:  # straight python
         for i in range(len(rnd)):
             pat=rnd[i]
             xx=reshape(x[:,pat],(1,-1))
             y=matrixmultiply(xx,w)
             w=w+params['eta']*(y*transpose(xx)-y**2*w);
             th=th+(1.0/params['tau'])*(y**2-th);
     else: # pyrex
         dohebb(params,w,th,x,rnd)
     print time.time()-t1


p.plot(w,'o-')
p.xlabel('weights')
p.show()


#=============================================

#         dohebb.pyx

cdef extern from "Numeric/arrayobject.h":

   struct PyArray_Descr:
     int type_num, elsize
     char type

   ctypedef class Numeric.ArrayType [object PyArrayObject]:
     cdef char *data
     cdef int nd
     cdef int *dimensions, *strides
     cdef object base
     cdef PyArray_Descr *descr
     cdef int flags


def dohebb(params,ArrayType w,ArrayType th,ArrayType X,ArrayType rnd):


     cdef int num_iterations
     cdef int num_inputs
     cdef int offset
     cdef double *wp,*xp,*thp
     cdef int *rndp
     cdef double eta,tau

     eta=params['eta']  # learning rate
     tau=params['tau']  # used for variance estimate

     cdef double y
     num_iterations=rnd.dimensions[0]
     num_inputs=w.dimensions[0]

     # get the pointers
     wp=&lt;double *&gt;w.data
     xp=&lt;double *&gt;X.data
     rndp=&lt;int *&gt;rnd.data
     thp=&lt;double *&gt;th.data

     for it from 0 &lt;= it &lt; num_iterations:

         offset=rndp[it]*num_inputs

         # calculate the output
         y=0.0
         for i from 0 &lt;= i &lt; num_inputs:
             y=y+wp[i]*xp[i+offset]

         # change in the weights
         for i from 0 &lt;= i &lt; num_inputs:
             wp[i]=wp[i]+eta*(y*xp[i+offset] - y*y*wp[i])

         # estimate the variance
         thp[0]=thp[0]+(1.0/tau)*(y**2-thp[0])







-------------------------------------------------------
This SF.net email is sponsored by: Splunk Inc. Do you grep through log files
for problems?  Stop!  Download the new AJAX search engine that makes
searching your log files as easy as surfing the  web.  DOWNLOAD SPLUNK!
<a class="moz-txt-link-freetext" href="http://sel.as-us.falkag.net/sel?cmd=lnk&kid=103432&bid=230486&dat=121642">http://sel.as-us.falkag.net/sel?cmd=lnk&amp;kid=103432&amp;bid=230486&amp;dat=121642</a>
_______________________________________________
Numpy-discussion mailing list
<a class="moz-txt-link-abbreviated" href="mailto:Num...@li...">Num...@li...</a>
<a class="moz-txt-link-freetext" href="https://lists.sourceforge.net/lists/listinfo/numpy-discussion">https://lists.sourceforge.net/lists/listinfo/numpy-discussion</a>

      </pre>
    </blockquote>
    <pre wrap="">
-------------------------------------------------------
This SF.net email is sponsored by: Splunk Inc. Do you grep through log files
for problems?  Stop!  Download the new AJAX search engine that makes
searching your log files as easy as surfing the  web.  DOWNLOAD SPLUNK!
<a class="moz-txt-link-freetext" href="http://sel.as-us.falkag.net/sel?cmd=k&kid3432&bid#0486&dat1642">http://sel.as-us.falkag.net/sel?cmd=k&amp;kid3432&amp;bid#0486&amp;dat1642</a>
_______________________________________________
Numpy-discussion mailing list
<a class="moz-txt-link-abbreviated" href="mailto:Num...@li...">Num...@li...</a>
<a class="moz-txt-link-freetext" href="https://lists.sourceforge.net/lists/listinfo/numpy-discussion">https://lists.sourceforge.net/lists/listinfo/numpy-discussion</a>




    </pre>
  </blockquote>
  <pre wrap=""><!---->
  </pre>
</blockquote>
</body>
</html>

Re: [Numpy-discussion] Timings for various round functions

[Bill B. <wb...@gm...>]

On 2/23/06, Travis Oliphant <oli...@ie...> wrote:
>Sasha wrote:

>I think we have agreed that C99 functions are good candidates to become
>ufuncs.  The only problem is figuring out what to do on platforms that
>don't define them.

>For example, we could define a separate module of C99 functions that is
>only available on certain platforms.


Might this be some help?  http://mega-nerd.com/FPcast/

--Bill

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?

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

Hello

The problem is that repeat is not the same as repmat.

For example:

>> repmat([1 2; 3 4], 2, 1)
ans =3D
     1     2
     3     4
     1     2
     3     4

In [12]: repeat(array([[1, 2],[3,4]]), 2)
Out[12]:
array([[1, 2],
       [1, 2],
       [3, 4],
       [3, 4]])

How can I use the repeat function as is to accomplish this?

Thanks.

Regards

Albert


On 2/23/06, Nadav Horesh <na...@vi...> wrote:
> You should really use the "repeat" function.
>
>    Nadav.

<snip>

RE: [Numpy-discussion] repmat equivalent?

[Nadav H. <na...@vi...>]

You should really use the "repeat" function.

   Nadav.


-----Original Message-----
From:	num...@li... on behalf of Albert =
Strasheim
Sent:	Thu 23-Feb-06 13:31
To:	num...@li...
Cc:=09
Subject:	[Numpy-discussion] repmat equivalent?
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 =
solve:

- 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


-------------------------------------------------------
This SF.Net email is sponsored by xPML, a groundbreaking scripting =
language
that extends applications into web and mobile media. Attend the live =
webcast
and join the prime developer group breaking into this new coding =
territory!
http://sel.as-us.falkag.net/sel?cmd=3Dk&kid=110944&bid$1720&dat=121642
_______________________________________________
Numpy-discussion mailing list
Num...@li...
https://lists.sourceforge.net/lists/listinfo/numpy-discussion

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

[Numpy-discussion] system_info problem

[Robert C. <cim...@nt...>]

Hi,

I am trying to (finally) move UMFPACK out of the sandbox to the scipy 
proper and so I need checking for its libraries via system_info.py of 
numpy distutils.

I have added a new section to site.cfg:

[umfpack]
library_dirs = 
/home/share/software/packages/UMFPACK/UMFPACK/Lib:/home/share/software/packages/UMFPACK/AMD/Lib
include_dirs = /home/share/software/packages/UMFPACK/UMFPACK/Include
umfpack_libs = umfpack, amd

the names of the libraries are libumfpack.a, libamd.a - they are 
correctly found in 'system_info._check_libs()' by 
'self._lib_list(lib_dir, libs, exts)', but then the function fails, 
since 'len(found_libs) == len(libs)' which is wrong. Can some 
numpy.distutils expert help me? Below is the new umfpack_info class I 
have written using the blas_info class as template.

yours clueless,
r.

PS: I do this because I prefer having the umfpack installed separately. 
It will be used, if present, to replace the default superLU-based sparse 
solver. Moving its sources under scipy/Lib/sparse would solve this 
issue, but Tim Davis recently changed the license of UMFPACK to GPL, and 
so the last version available for the direct inclusion is 4.4. (4.6 is 
the current one). Opinions are of course welcome.

--
class umfpack_info(system_info):
     section = 'umfpack'
     dir_env_var = 'UMFPACK'
     _lib_names = ['umfpack', 'amd']
     includes = 'umfpack.h'
     notfounderror = UmfpackNotFoundError

     def calc_info(self):
         info = {}

         lib_dirs = self.get_lib_dirs()
         print lib_dirs
         umfpack_libs = self.get_libs('umfpack_libs', self._lib_names)
         print umfpack_libs
         for d in lib_dirs:
             libs = self.check_libs(d,umfpack_libs,[])
             print d, libs
             if libs is not None:
                 dict_append(info,libraries=libs)
                 break
         else:
             return

         include_dirs = self.get_include_dirs()
         print include_dirs
         h = (self.combine_paths(lib_dirs+include_dirs,includes) or 
[None])[0]
         if h:
             h = os.path.dirname(h)
             dict_append(info,include_dirs=[h])

         print info

         self.set_info(**info)

[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] Multiple inheritance from ndarray

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

On Wed, Feb 22, 2006 at 08:58:28PM -0700, Travis Oliphant wrote:
> Here's an outline of what you need to do.  This is, of course,=20
> untested....  For example, I don't really know what actImage is.
>=20
> from numpy import ndarray, array
>=20
> class Image(ndarray, actImage):
>    def __new__(subtype, *args)
>        act1 =3D actImage.__new__(actImage, *args)
>        actImage.__init__(act1, *args)
>        arr =3D array(act1.getArray(), 'd', copy=3DFalse)
>        self =3D arr.view(subtype)
>        # you might need to copy attributes from act1 over to self here.=
..
>        return self

This is probably the right place to use super, i.e.:

def __new__(subtype, *args):
	act1 =3D super(Image, subtype).__new__(subtype, *args)
	...

def __init__(self, *args):
	super(Image, self).__init__(*args)

The attached script shows how multiple inheritance runs through
different classes.

St=E9fan

[Numpy-discussion] Success with SVN on Cygwin

[Jose Gomez-D. <jos...@gm...>]

Hi!
A few days ago, I asked about compiling NumPy on Cygwin. Travis carried out some
modifications, and with last night's SVN, I can happily report that it now
compiles and works. The tests produced no errors, so it's all good :)

Many thanks to all, and to Travis specially, for his really fast response.

Many thanks!
Jose

Re: [Numpy-discussion] Re: Why floor and ceil change the type of the array?

[Sasha <nd...@ma...>]

On 2/23/06, Travis Oliphant <oli...@ie...> wrote:
> How do you propose to determine if the output fits into the data-type?
> Are you proposing to have different output rules for different
> functions.  Sheer madness...  The rules now are (relatively) simple and
> easy to program to.
>
I did not propose that.  I just mentioned the output to argue that the
rule to use minimal floating point type that can represent the input
is an arbitrary one and no better than cast all integers to doubles.=20
"Sheer madness...", however is too strong a characterization.  Note
that python (so far) changes the type of the result depending on its
value in some cases:

>>> type(2**30)
<type 'int'>
>>> type(2**32)
<type 'long'>

This is probably unacceptable to numpy for performance reasons, but it
is not madness.

Try to explain the following to someone who is used to python arithmetics:

>>> 2*array(2**62)+array(2*2**62)
0L

> Hardly a good example.  Are you also concerned about the following?
>
>  >>> exp(1000)
> inf
>
>  >>> exp(array(1000,'g'))
> 1.97007111401704699387e+434

No, but I think it is because I am conditioned by C.  To me exp() is a
double-valued function that happened to work for ints with the help of
an implicit cast.  You may object that this is so because C does not
allow function overloading, but C++ does overload exp so that
exp((float)1) is float, and exp((long double)1) is long doubles but
exp((short)1), exp((char)1) and exp((long long)1) are all double.

Both numpy and C++ made an arbitrary design choice.  I find C++ choice
simpler and more natural, but I can live with the numpy choice once
I've learned what it is.

Re: [Numpy-discussion] Thoughts on an ndarray super-class

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

Sasha wrote:

>On 2/23/06, Travis Oliphant <oli...@ie...> wrote:
>  
>
>>...
>>I have been thinking, however, of replacing it with a super-class that
>>does not define the dimensions or strides.
>>
>>    
>>
>Having a simple 1-d array in numpy would be great.  In an ideal world
>I would rather see a 1-d array implemented in C together with a set of
>array operations that is rich enough to allow trivial implementation
>of ndarray in pure python.
>  
>
You do realize that this is essentially numarray, right?  And your dream 
of *rich enough* 1-d operations to allow *trivial* implementation may be 
a bit far-fetched, but I'm all for dreaming.

>When you say "does not define the dimensions or strides" do you refer
>to python interface or to C struct?  I thought python did not allow to
>add data members to object structs in subclasses.
>  
>
The C-struct.  

Yes, you can add data-members to object structs in sub-classes.  Every 
single Python Object does it.  The standard Python Object just defines 
PyObject_HEAD or PyObject_VAR_HEAD.  

This is actually the essence of inheritance in C and it is why 
subclasses written in C must have compatible memory layouts.  The first 
part of the C-structure must be identical, but you can add to it all you 
want.

It all comes down to:  Can I cast to the base-type C-struct and have 
everything still work out when I dereference a particular field? 

This will be true if  PyArrayObject is

struct {
PyBaseArrayObject
int nd
intp *dimensions
intp *strides
}

I suppose we could change the

int nd

to

intp nd

and place it in the PyBaseArrayObject where it would be used as a length

But, I don't really like that...

>>In other words, the default array would be just a block of memory.  The
>>standard array would inherit from the default and add dimension and
>>strides pointers.
>>
>>    
>>
>If python lets you do it, how will that block of memory know its size?
>
>
>  
>
It won't of course by itself unless you add an additional size field.  
Thus, I'm not really sure whether it's a good idea or not.   I don't 
like the idea of adding more and more fields to the basic C-struct that 
has been around for 10 years unless we have a good reason.   The other 
issue is that the data-pointer doesn't always refer to memory that the 
ndarray has allocated, so it's actually incorrect to think of the 
ndarray as both the block of memory and the dimensioned indexing. 

The memory pointer is just that (a memory pointer).  We are currently 
allowing ndarray's to create their own memory but that could easily 
change so that they always use some other object to allocate memory.

In short, I don't see how to really do it so that the base object is 
actually useable.

Re: [Numpy-discussion] Timings for various round functions

[Sasha <nd...@ma...>]

On 2/23/06, Travis Oliphant <oli...@ie...> wrote:
> ...
> I think we have agreed that C99 functions are good candidates to become
> ufuncs.  The only problem is figuring out what to do on platforms that
> don't define them.
>

I was going to ask this question myself, but then realized that the
answer is in the source code: for functions missing on a platform
numpy provides its own implementations.  (See for example a comment in
umathmodule "if C99 extensions not available then define dummy
functions...") I was going to just use rint instead of round and
nearbyint on platforms that dont have them.

> For example, we could define a separate module of C99 functions that is
> only available on certain platforms.

This is certainly the easiest to implement option, but we don't want
make numpy users  worry about portability of their code.

Re: [Numpy-discussion] Re: Why floor and ceil change the type of the array?

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

Sasha wrote:

>On 2/23/06, Robert Kern <rob...@gm...> wrote:
>  
>
>>>I cannot really think of any reason for the current numpy behaviour
>>>other than the consistency with transcendental functions.
>>>      
>>>
>>It's simply the easiest thing to do with the ufunc machinery.
>>    
>>
>>AFAICT, the story goes like this: sin() has two implementations, one for
>>single-precision floats and one for doubles. The ufunc machinery sees the int16
>>and picks single-precision as the smallest type of the two that can fit an int16
>>without losing precision. Naturally, you probably want the function to operate
>>in higher precision, but that's not really information that the ufunc machinery
>>knows about.
>>    
>>
>
>According to your theory long (i8) integers should cast to long doubles, but
>  
>
Robert is basically right, except there is a special case for long 
integers because long doubles are not cross platform.  The relevant code 
is PyArray_CanCastSafely.  This is basically the coercion rule  table.   
You will notice the special checks for long double placed there after it 
was noticed that on 64-bit platforms long doubles were cropping up an 
awful lot and it was decided that because long doubles are not very 
ubiquitous (for example many platforms don't distinguish between long 
double and double), we should special-case the 64-bit integer rule. 

You can read about it in the archives if you want. 

>dtype('<f8')
>
>
>Given that python's floating point object is a double, I think it
>would be natural to cast integer arguments to double for all sizes.
>
Perhaps, but that is not what is done.  I don't think it's that big a 
deal because to get "different size" integers you have to ask for them 
and then you should know that conversion to floating point is not 
necessarily a double. 

I think the only accetable direction to pursue is to raise an error and 
not do automatic upcasting if a ufunc does not have a definition for any 
of the given types.   But, this is an old behavior from Numeric, and I 
would think such changes now would rightly be considered as gratuitous 
breakage.

>  I
>would also think that in choosing the precision for a function it is
>also important that the output fits into data type.  
>
How do you propose to determine if the output fits into the data-type?  
Are you proposing to have different output rules for different 
functions.  Sheer madness...  The rules now are (relatively) simple and 
easy to program to.

>I find the
>following unfortunate:
>  
>
>>>>exp(400)
>>>>        
>>>>
>5.2214696897641443e+173
>  
>
>>>>exp(array(400,'h'))
>>>>        
>>>>
>inf
>  
>
Hardly a good example.  Are you also concerned about the following?

 >>> exp(1000)
inf

 >>> exp(array(1000,'g'))
1.97007111401704699387e+434

Re: [Numpy-discussion] Thoughts on an ndarray super-class

[Sasha <nd...@ma...>]

On 2/23/06, Travis Oliphant <oli...@ie...> wrote:
> ...
> I have been thinking, however, of replacing it with a super-class that
> does not define the dimensions or strides.
>
Having a simple 1-d array in numpy would be great.  In an ideal world
I would rather see a 1-d array implemented in C together with a set of
array operations that is rich enough to allow trivial implementation
of ndarray in pure python.

When you say "does not define the dimensions or strides" do you refer
to python interface or to C struct?  I thought python did not allow to
add data members to object structs in subclasses.


> In other words, the default array would be just a block of memory.  The
> standard array would inherit from the default and add dimension and
> strides pointers.
>
If python lets you do it, how will that block of memory know its size?

[Numpy-discussion] Re: Why floor and ceil change the type of the array?

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

Sasha wrote:
> On 2/23/06, Robert Kern <rob...@gm...> wrote:

>>AFAICT, the story goes like this: sin() has two implementations, one for
>>single-precision floats and one for doubles. The ufunc machinery sees the int16
>>and picks single-precision as the smallest type of the two that can fit an int16
>>without losing precision. Naturally, you probably want the function to operate
>>in higher precision, but that's not really information that the ufunc machinery
>>knows about.
> 
> According to your theory long (i8) integers should cast to long doubles, but
> 
>>>>sin(array(0,'i8')).dtype
> 
> dtype('<f8')

Out of curiosity, what does array(0, dtype=longdouble).dtype give you on your
platform?

The relevant function here is select_types() in ufuncobject.c.

> Given that python's floating point object is a double, I think it
> would be natural to cast integer arguments to double for all sizes.  I
> would also think that in choosing the precision for a function it is
> also important that the output fits into data type.  I find the
> following unfortunate:
> 
>>>>exp(400)
> 
> 5.2214696897641443e+173
> 
>>>>exp(array(400,'h'))
> 
> inf

I prefer consistent, predictable rules that are dependent on the input, not the
output. If I want my outputs to be double precision, I will cast appopriately.

-- 
Robert Kern
rob...@gm...

"In the fields of hell where the grass grows high
 Are the graves of dreams allowed to die."
  -- Richard Harter