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

Re: [Numpy-discussion] Experience with Visit?

[Torgil S. <tor...@gm...>]

Sourceforge mailing list doesn't like me. Here's another try.

I've managed to use visit through it's cli and giving it data through
pyvtk so that this code

if __name__=='__main__':
    from numpy import *
    from numpy.random import randn
    n=10000
    t=linspace(0,1,n)+randn(n)*0.01
    a=t+t*(1+randn(n)*0.1)
    phi=randn(n)*2*pi*0.05
    scatter4( a*cos(2.0*2*pi*t+phi) , a*sin(1.5*2*pi*t+phi) , t
,1.0/sqrt(1+t**2) )

gives this result:
http://www.torgil.org/visit/scatter4.png

You can grab the source using git:  (see http://git.or.cz/ for details)
git clone http://torgil.org/pyvisit.git

I'm using the pty module in python to communicate with the visit cli
which unfortunately is only available on Unix. If you do any updates,
like an usefule API or windows support, please share (that's why I put
it in git due to it's excellent distributed development capabilities).
This  was more of a proof of concept for me.

> On 9/15/06, Rob Hetland <he...@ta...> wrote:
> >
> > I have looked at it, and I even know somebody who is using python to
> > put in a data reader for FVCOM (a finite volume ocean model).  He has
> > a prototype working.  I have no experience with the python bindings
> > personally.
> >
> > I looked into compiling it from source (for an intel mac), but it was
> > daunting so I gave up.
> >
> > I also contacted them to try to put in lower level python bindings,
> > so that we could create 3D graphics from the command line, but they
> > didn't seem too interested in that.
> >
> > Otherwise, it is a very nice workable app.
> >
> > -Rob
> >
> >
> > On Sep 15, 2006, at 11:23 AM, O'Keefe, Michael wrote:
> >
> > > I haven't tried VisIT before but thanks for the link. I also
> > > downloaded and am checking it out.
> > >
> > > Along this same line of discussion, has anyone tried OOF2 which is
> > > an FEA package that also has some strong python connections?
> > >
> > > http://www.ctcms.nist.gov/oof/oof2.html
> > >
> > > I'm working on a Windows machine and it the current code-base
> > > doesn't seem to support Windows out of the box (if at all). Looks
> > > like you can put it together for *nix or Mac OS X, though...
> > >
> > > --Michael
> > >
> > >> -----Original Message-----
> > >> From: num...@li...
> > >> [mailto:num...@li...] On
> > >> Behalf Of Robert Cimrman
> > >> Sent: Friday, September 15, 2006 8:29
> > >> To: Discussion of Numerical Python
> > >> Subject: Re: [Numpy-discussion] Experience with Visit?
> > >>
> > >> Travis Oliphant wrote:
> > >>> Has anybody had any experience with the 3-D visualization software
> > >>> VISIT?   It has Python bindings and seems to be pretty
> > >> sophisticated.
> > >>> I'm wondering why I haven't heard more about it.
> > >>>
> > >>> http://www.llnl.gov/visit/
> > >>
> > >> No reaction up to now, so...
> > >>
> > >> I have just tried the 'getting started' part and was quite impressed,
> > >> thanks for posting the link! Up to now I have used ParaView
> > >> and was very
> > >> satisfied, but the Python bindings of VisIt are a great lure.
> > >>
> > >> r.
> > >>
> > >> --------------------------------------------------------------
> > >> -----------
> > >> 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
> > >> Num...@li...
> > >> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
> > >>
> > >
> > > ----------------------------------------------------------------------
> > > ---
> > > 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
> > > Num...@li...
> > > https://lists.sourceforge.net/lists/listinfo/numpy-discussion
> >
> > ----
> > Rob Hetland, Associate Professor
> > Dept. of Oceanography, Texas A&M University
> > http://pong.tamu.edu/~rob
> > phone: 979-458-0096, fax: 979-845-6331
> >
> >
> >
> > -------------------------------------------------------------------------
> > 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
> > Num...@li...
> > https://lists.sourceforge.net/lists/listinfo/numpy-discussion
> >
>

Re: [Numpy-discussion] Problems with Numexpr and discontiguous arrays

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

Ivan Vilata i Balaguer wrote:
> It seemed that discontiguous arrays worked OK in Numexpr since r1977 or
> so, but I have come across some alignment or striding problems which can
> be seen with the following code::
>
>     import numpy
>     import numexpr
>
>     array_length = 10
>     array_descr = [('c1', numpy.int32), ('c2', numpy.uint16)]
>
>     array = numpy.empty((array_length,), dtype=array_descr)
>     for i in xrange(array_length):
>         array['c1'][i] = i
>         array['c2'][i] = 0xaaaa
>
>     print numexpr.evaluate('c1', {'c1': array['c1']})
>     print numexpr.evaluate('c1', {'c1': array['c1'].copy()})
>
> Im my computer, Pentium IV with NumPy 1.0rc1 and Numexpr r2239
> (unmodified) this gives the following result::
>
>     [      0      109226 -1431699456           2      240298 -1431699456
>            4      371370           8      633514]
>     [0 1 2 3 4 5 6 7 8 9]
>
> The test works right when ``evaluate()`` is used with 'c2' instead of
> 'c1', and also when 'c2' also measures 32 bits and fields are aligned.
> Maybe the ``memsteps`` value is not getting used somewhere.  Any ideas
> on this?
>   
I suspect that there are some assumptions that the element separation  
is an integral multiple of the element size. I certainly didn't have 
record arrays in mind when I was working on the striding stuff, so it 
wouldn't surprise me. This should be fixed: preferably to do the right 
thing and at a minimum to cleanly raise an exception rather than 
spitting out garbage. I don't know that I'll have time to mess with it 
soon though.

-tim

Re: [Numpy-discussion] speeding up y[:, i] for y.shape = (big number, small number)

[Bruce S. <bso...@gm...>]

Hi,
I think what you are after is the multivariate normal distribution.
Assuming I have it correctly, it is clearer to see (and probably more
accurate to compute) in the log form as:

-(N/2)*log(2*PI) - 0.5*log(determinant of V) - 0.5*(transpose of
(x-mu))*inverse(V)*(x-mu)

where N is the number of observations, PI is math constant, V is the
known variance co-variance matrix, x is vector of values,  mu is the
known mean.

If so, then you can vectorize the density calculation. Of course if V
has a simple structure (it is unclear from your code on this) and can
be factored out, then the only 'hard' part to compute is the product
of the vector transpose with the vector. If mu and va are scalars for
each calculation, then you can factor them out for so the only
multidimensional calculations are squares and sums of each frame.

Regards
Bruce




On 10/4/06, David Cournapeau <da...@ar...> wrote:
> Tim Hochberg wrote:
> >>
> >>
> >> I guess the problem is coming from the fact that y being C order, y[:,
> >> i] needs accessing data in a non 'linear' way. Is there a way to speed
> >> this up ? I did something like this:
> >>
> >>   y   = N.zeros((K, n))
> >>     for i in range(K):
> >>         y[i] = gauss_den(data, mu[i, :], va[i, :])
> >>     return y.T
> >>
> >> which works, but I don't like it very much.
> > Why not?
> >
> Mainly because using those transpose do not really reflect the
> intention, and this does not seem natural.
> >>  Isn't there any other way
> > That depends on the details of gauss_den.
> >
> > A general note: for this kind of microoptimization puzzle, it's much
> > easier to help if you can post a self contained example, preferably
> > something fairly simple that still illustrates the speed issue, that we
> > can experiment with.
> >
> Here we are (the difference may not seem that much between the two
> multiple_ga, but in reality, _diag_gauss_den is an internal function
> which is done in C, and is much faster... By writing this example, I've
> just realized that the function _diag_gauss_den may be slow for exactly
> the same reasons):
>
>
> #! /usr/bin/env python
> # Last Change: Wed Oct 04 07:00 PM 2006 J
>
> import numpy as N
> from numpy.random import randn
>
> def _diag_gauss_den(x, mu, va):
>     """ This function is the actual implementation
>     of gaussian pdf in scalar case. It assumes all args
>     are conformant, so it should not be used directly
>
>     Call gauss_den instead"""
>     # Diagonal matrix case
>     d       = mu.size
>     inva    = 1/va[0]
>     fac     = (2*N.pi) ** (-d/2.0) * N.sqrt(inva)
>     y       =  (x[:,0] - mu[0]) ** 2 * inva * -0.5
>     for i in range(1, d):
>         inva    = 1/va[i]
>         fac     *= N.sqrt(inva)
>         y       += (x[:,i] - mu[i]) ** 2 * inva * -0.5
>     y   = fac * N.exp(y)
>
> def multiple_gauss_den1(data, mu, va):
>     """Helper function to generate several Gaussian
>     pdf (different parameters) from the same data: unoptimized version"""
>     K   = mu.shape[0]
>     n   = data.shape[0]
>     d   = data.shape[1]
>
>     y   = N.zeros((n, K))
>     for i in range(K):
>         y[:, i] = _diag_gauss_den(data, mu[i, :], va[i, :])
>
>     return y
>
> def multiple_gauss_den2(data, mu, va):
>     """Helper function to generate several Gaussian
>     pdf (different parameters) from the same data: optimized version"""
>     K   = mu.shape[0]
>     n   = data.shape[0]
>     d   = data.shape[1]
>
>     y   = N.zeros((K, n))
>     for i in range(K):
>         y[i] = _diag_gauss_den(data, mu[i, :], va[i, :])
>
>     return y.T
>
> def bench():
>     #===========================================
>     # GMM of 30 comp, 15 dimension, 1e4 frames
>     #===========================================
>     d       = 15
>     k       = 30
>     nframes = 1e4
>     niter   = 10
>     mode    = 'diag'
>
>     mu      = randn(k, d)
>     va      = randn(k, d) ** 2
>     X       = randn(nframes, d)
>
>     print "============================================================="
>     print "(%d dim, %d components) GMM with %d iterations, for %d frames" \
>             % (d, k, niter, nframes)
>
>     for i in range(niter):
>         y1  = multiple_gauss_den1(X, mu, va)
>
>     for i in range(niter):
>         y2  = multiple_gauss_den2(X, mu, va)
>
>     se  = N.sum(y1-y2)
>
>     print se
>
> if __name__ == '__main__':
>     import hotshot, hotshot.stats
>     profile_file    = 'foo.prof'
>     prof    = hotshot.Profile(profile_file, lineevents=1)
>     prof.runcall(bench)
>     p = hotshot.stats.load(profile_file)
>     print p.sort_stats('cumulative').print_stats(20)
>     prof.close()
>
> I am a bit puzzled by all those C vs F storage, though. In Matlab, where
> the storage was always F as far as I know, I have never encountered such
> differences (eg between y(:, i) and y(:, i)); I don't know if this is
> because I am doing it badly, or because matlab is much more clever than
> numpy at handling those cases, or if that is the price to pay for the
> added flexibility of numpy...
>
> David
>
> -------------------------------------------------------------------------
> 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
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>

Re: [Numpy-discussion] Vectorizing code, for loops, and all that

[Bruce S. <bso...@gm...>]

Hi,

On 10/2/06, Albert Strasheim <fu...@gm...> wrote:
> I recently started looking at David Cournapeau's PyEM package,
> specifically his implementation of the K-Means algorithm. He
> implemented part of this algorithm with in pure Python version and
> also provided a Pyrex alternative that is significantly faster (about
> 10 times with the data I tested with). I tried to figure out why this
> is so.

Just of of curosity, have you seen the implementation at BioPython
(http://biopython.org)?

Michiel de Hoon and other wrote 'The C Clustering Library', see the
following document for more details (linked from documentation page -
http://biopython.org/wiki/Documentation):
http://biopython.org/DIST/docs/cluster/cluster.pdf

Regards
Bruce

[Numpy-discussion] Problems with Numexpr and discontiguous arrays

[Ivan V. i B. <iv...@ca...>]

It seemed that discontiguous arrays worked OK in Numexpr since r1977 or
so, but I have come across some alignment or striding problems which can
be seen with the following code::

    import numpy
    import numexpr

    array_length =3D 10
    array_descr =3D [('c1', numpy.int32), ('c2', numpy.uint16)]

    array =3D numpy.empty((array_length,), dtype=3Darray_descr)
    for i in xrange(array_length):
        array['c1'][i] =3D i
        array['c2'][i] =3D 0xaaaa

    print numexpr.evaluate('c1', {'c1': array['c1']})
    print numexpr.evaluate('c1', {'c1': array['c1'].copy()})

Im my computer, Pentium IV with NumPy 1.0rc1 and Numexpr r2239
(unmodified) this gives the following result::

    [      0      109226 -1431699456           2      240298 -1431699456
           4      371370           8      633514]
    [0 1 2 3 4 5 6 7 8 9]

The test works right when ``evaluate()`` is used with 'c2' instead of
'c1', and also when 'c2' also measures 32 bits and fields are aligned.
Maybe the ``memsteps`` value is not getting used somewhere.  Any ideas
on this?

::

	Ivan Vilata i Balaguer   >qo<   http://www.carabos.com/
	       C=C3=A1rabos Coop. V.  V  V   Enjoy Data
	                          ""

Re: [Numpy-discussion] speeding up y[:, i] for y.shape = (big number, small number)

[David C. <da...@ar...>]

David Cournapeau wrote:
> Here we are (the difference may not seem that much between the two 
> multiple_ga, but in reality, _diag_gauss_den is an internal function 
> which is done in C, and is much faster... By writing this example, I've 
> just realized that the function _diag_gauss_den may be slow for exactly 
> the same reasons):
>
>   
I checked my assumption about the _diag_gauss_den function being slow 
because of the same problem, and this is indeed the case:

If I replace X = randn(nframes, d) by X = randn(d, nframes); X   = X.T, 
the function is mode than twice as fast ! This seems way to much of a 
difference to me....

David

Re: [Numpy-discussion] speeding up y[:, i] for y.shape = (big number, small number)

[David C. <da...@ar...>]

Tim Hochberg wrote:
>>
>>
>> I guess the problem is coming from the fact that y being C order, y[:, 
>> i] needs accessing data in a non 'linear' way. Is there a way to speed 
>> this up ? I did something like this:
>>
>>   y   = N.zeros((K, n))
>>     for i in range(K):
>>         y[i] = gauss_den(data, mu[i, :], va[i, :])
>>     return y.T
>>
>> which works, but I don't like it very much.
> Why not?
>
Mainly because using those transpose do not really reflect the 
intention, and this does not seem natural.
>>  Isn't there any other way 
> That depends on the details of gauss_den.
>
> A general note: for this kind of microoptimization puzzle, it's much 
> easier to help if you can post a self contained example, preferably 
> something fairly simple that still illustrates the speed issue, that we 
> can experiment with.
>
Here we are (the difference may not seem that much between the two 
multiple_ga, but in reality, _diag_gauss_den is an internal function 
which is done in C, and is much faster... By writing this example, I've 
just realized that the function _diag_gauss_den may be slow for exactly 
the same reasons):


#! /usr/bin/env python
# Last Change: Wed Oct 04 07:00 PM 2006 J

import numpy as N
from numpy.random import randn

def _diag_gauss_den(x, mu, va):
    """ This function is the actual implementation
    of gaussian pdf in scalar case. It assumes all args
    are conformant, so it should not be used directly
   
    Call gauss_den instead"""
    # Diagonal matrix case
    d       = mu.size
    inva    = 1/va[0]
    fac     = (2*N.pi) ** (-d/2.0) * N.sqrt(inva)
    y       =  (x[:,0] - mu[0]) ** 2 * inva * -0.5
    for i in range(1, d):
        inva    = 1/va[i]
        fac     *= N.sqrt(inva)
        y       += (x[:,i] - mu[i]) ** 2 * inva * -0.5
    y   = fac * N.exp(y)

def multiple_gauss_den1(data, mu, va):
    """Helper function to generate several Gaussian
    pdf (different parameters) from the same data: unoptimized version"""
    K   = mu.shape[0]
    n   = data.shape[0]
    d   = data.shape[1]
   
    y   = N.zeros((n, K))
    for i in range(K):
        y[:, i] = _diag_gauss_den(data, mu[i, :], va[i, :])

    return y

def multiple_gauss_den2(data, mu, va):
    """Helper function to generate several Gaussian
    pdf (different parameters) from the same data: optimized version"""
    K   = mu.shape[0]
    n   = data.shape[0]
    d   = data.shape[1]
   
    y   = N.zeros((K, n))
    for i in range(K):
        y[i] = _diag_gauss_den(data, mu[i, :], va[i, :])
   
    return y.T

def bench():
    #===========================================
    # GMM of 30 comp, 15 dimension, 1e4 frames
    #===========================================
    d       = 15
    k       = 30
    nframes = 1e4
    niter   = 10
    mode    = 'diag'

    mu      = randn(k, d)
    va      = randn(k, d) ** 2
    X       = randn(nframes, d)

    print "============================================================="
    print "(%d dim, %d components) GMM with %d iterations, for %d frames" \
            % (d, k, niter, nframes)

    for i in range(niter):
        y1  = multiple_gauss_den1(X, mu, va)

    for i in range(niter):
        y2  = multiple_gauss_den2(X, mu, va)

    se  = N.sum(y1-y2)

    print se

if __name__ == '__main__':
    import hotshot, hotshot.stats
    profile_file    = 'foo.prof'
    prof    = hotshot.Profile(profile_file, lineevents=1)
    prof.runcall(bench)
    p = hotshot.stats.load(profile_file)
    print p.sort_stats('cumulative').print_stats(20)
    prof.close()

I am a bit puzzled by all those C vs F storage, though. In Matlab, where 
the storage was always F as far as I know, I have never encountered such 
differences (eg between y(:, i) and y(:, i)); I don't know if this is 
because I am doing it badly, or because matlab is much more clever than 
numpy at handling those cases, or if that is the price to pay for the 
added flexibility of numpy...

David

Re: [Numpy-discussion] "Manually" broadcasting arrays in Python

[A. M. A. <per...@gm...>]

On 04/10/06, Robert Kern <rob...@gm...> wrote:

> Yeah, a segfault would be problematic. Otherwise, it works and is in fact faster
> than what I wrote.

It's a bit tricky to trigger but I think it was fixed in 1.0rc1 (in
changeset 3125, in fact:
http://projects.scipy.org/scipy/numpy/changeset/3125 )

Would it be useful for me to contribute the tiny script I wrote to
trigger it as a regression test?

A. M. Archibald

Re: [Numpy-discussion] "Manually" broadcasting arrays in Python

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

A. M. Archibald wrote:
> On 03/10/06, Robert Kern <rob...@gm...> wrote:
>> Has anyone implemented an easier or more efficient way to broadcast arrays to a
>> common shape at the Python level? I was hoping that the broadcast iterator would
>> actually provide the broadcasted arrays, but it does not.
> 
> How about vectorize(lambda *args: args)? Almost works, only it
> sometimes segfaults with more than two arguments... but that's clearly
> a numpy bug.

Yeah, a segfault would be problematic. Otherwise, it works and is in fact faster 
than what I wrote.

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
  that is made terrible by our own mad attempt to interpret it as though it had
  an underlying truth."
   -- Umberto Eco

Re: [Numpy-discussion] Vectorizing code, for loops, and all that

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

Hello all

> -----Original Message-----
> From: Travis Oliphant [mailto:oli...@ee...]
> Sent: 03 October 2006 03:51
> To: Discussion of Numerical Python; Albert Strasheim
> Subject: Re: ***[Possible UCE]*** Re: [Numpy-discussion] Vectorizing code,
> for loops, and all that
> 
> <snip>

> >Meanwhile, I can confirm that the NOBUFFER_UFUNCLOOP case in
> >PyUFunc_GenericFunction is getting exercised in the slower case. Here's
> some
> >info on what's happening, courtesy of Rational Quantify:
> >
> >case NOBUFFER_UFUNCLOOP:
> >while (loop->index < loop->size) {
> >  for (i=0; i<self->nargs; i++)
> >    loop->bufptr[i] = loop->iters[i]->dataptr; [1]
> >
> >  loop->function((char **)loop->bufptr, &(loop->bufcnt),
> >                 loop->steps, loop->funcdata); [2]
> >  UFUNC_CHECK_ERROR(loop);
> >
> >  for (i=0; i<self->nargs; i++) {
> >      PyArray_ITER_NEXT(loop->iters[i]); [3]
> >  }
> >  loop->index++;
> >}
> >break;
> >
> >[1] 12.97% of function time
> >[2] 8.65% of functiont ime
> >[3] 62.14% of function time
> >
> >If statistics from elsewhere in the code would be helpful, let me know,
> and
> >I'll see if I can convince Quantify to cough it up.
> >
> Please run the same test but using
> 
> x1 = N.random.rand(39,2000)
> x2 = N.random.rand(39,64,1)
> 
> z1 = x1[:,N.newaxis,:] - x2

Very similar results to what I had previously:

[1] 10.88%
[2] 7.25%
[3] 68.25%

Cheers,

Albert

Re: [Numpy-discussion] Vectorizing code, for loops, and all that

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

Er, is this the PyEM in question?
It doesn't have much of a web presence...
   http://www.ar.media.kyoto-u.ac.jp/members/david/

--bb

On 10/3/06, David Cournapeau <da...@ar...> wrote:
> Albert Strasheim wrote:
> > Hello all
> >
> > I recently started looking at David Cournapeau's PyEM package, specifically
> > his implementation of the K-Means algorithm. He implemented part of this
> > algorithm with in pure Python version and also provided a Pyrex alternative
> > that is significantly faster (about 10 times with the data I tested with). I
> > tried to figure out why this is so.

Re: [Numpy-discussion] "Manually" broadcasting arrays in Python

[A. M. A. <per...@gm...>]

On 03/10/06, Robert Kern <rob...@gm...> wrote:
> Has anyone implemented an easier or more efficient way to broadcast arrays to a
> common shape at the Python level? I was hoping that the broadcast iterator would
> actually provide the broadcasted arrays, but it does not.

How about vectorize(lambda *args: args)? Almost works, only it
sometimes segfaults with more than two arguments... but that's clearly
a numpy bug.

A. M. Archibald

[Numpy-discussion] Vectorizing code, for loops, and all that

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

Hello all

I recently started looking at David Cournapeau's PyEM package,
specifically his implementation of the K-Means algorithm. He
implemented part of this algorithm with in pure Python version and
also provided a Pyrex alternative that is significantly faster (about
10 times with the data I tested with). I tried to figure out why this
is so.

The job of this part of the algorithm is pretty simple: from a bunch
of cluster means (the codebook) find the nearest cluster mean for each
data point. The first attempt at implementing this algorithm might use
two for loops, one over the data points and one over the cluster
means, computing a Euclidean distance at each step. Slightly better is
to use one loop and some broadcasting.

By randomly trying some code, I arrived at the following one-liner:

N.sum((data[...,N.newaxis,...]-code)**2, 2).argmin(axis=1)

where

data = N.random.randn(2000, 39)
nclusters = 64
code = data[0:nclusters,:]

This code was still slow, so I mangled it a bit so that I could
profile it using cProfile under Python 2.5 (use profile if you have an
older version of Python):

def kmean():
    data = N.random.randn(2000, 39)
    nclusters = 64
    code = data[0:nclusters,:]
    for i in xrange(10):
        def a(): return data[...,N.newaxis,...]
        z = a()
        def b(): return z-code
        z = b()
        def c(): return z**2
        z = c()
        def d(): return N.sum(z, 2)
        z = d()
        def e(): return z.argmin(axis=1)
        z = e()

if __name__ == '__main__':
    import cProfile
    cProfile.run("kmean()")

Profiler output:

In [570]: %run -t -e kmean2.py
         84 function calls in 8.567 CPU seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    8.567    8.567 <string>:1(<module>)
       10    0.000    0.000    0.835    0.084 fromnumeric.py:375(sum)
       10    0.000    0.000    0.000    0.000 kmean2.py:10(a)
       10    6.837    0.684    6.837    0.684 kmean2.py:12(b)
       10    0.623    0.062    0.623    0.062 kmean2.py:14(c)
       10    0.000    0.000    0.835    0.084 kmean2.py:16(d)
       10    0.080    0.008    0.080    0.008 kmean2.py:18(e)
        1    0.180    0.180    8.567    8.567 kmean2.py:5(kmean2)
       10    0.000    0.000    0.000    0.000 {isinstance}
        1    0.000    0.000    0.000    0.000 {method 'disable' of
'_lsprof.Profiler' objects}
        1    0.012    0.012    0.012    0.012 {method 'randn' of
'mtrand.RandomState' objects}
       10    0.835    0.083    0.835    0.083 {method 'sum' of
'numpy.ndarray' objects}

It seems the b part of the computation is taking most of the time. I
tried to simulate this separately:

In [571]: x1 = N.random.randn(2000,39)

In [572]: y1 = N.random.randn(64,39)

In [574]: %timeit z1=x1[...,N.newaxis,...]-y1 10 loops, best of 3: 703
ms per loop

In [575]: z1.shape
Out[575]: (2000, 64, 39)

As far as I can figure, this operation is doing 2000*64*39
subtractions. Doing this straight up yields the following:

In [576]: x2 = N.random.randn(2000,64,39)

In [577]: y2 = N.random.randn(2000,64,39)

In [578]: %timeit z2 = x2-y2
10 loops, best of 3: 108 ms per loop

Does anybody have any ideas on why this is so much faster? Hopefully I
didn't mess up somewhere...

Thanks!

Regards,

Albert

P.S. I'm used NumPy 1.0rc1 with Python 2.5 on Windows for these tests.
I get similar results with 1.0.dev3239 with Python 2.4.3 on Linux.

Re: [Numpy-discussion] Vectorizing code, for loops, and all that

[A. M. A. <per...@gm...>]

On 03/10/06, Tim Hochberg <tim...@ie...> wrote:

> I had an idea regarding which axis to operate on first. Rather than
> operate on strictly the longest axis or strictly the innermost axis, a
> hybrid approach could be used. We would operate on the longest axis that
> would not result in the inner loop overflowing the cache. The idea is
> minimize the loop overhead as we do now by choosing the largest axis,
> while at the same time attempting to maintain cache friendliness.

If elements are smaller than cache lines (usually at least eight
bytes, I think), we might end up pulling many times as many bytes into
the cache as we actually need if we don't loop along axes with small
strides first.

Can BLAS be used for some of these operations?

A. M. Archibald

A. M. Archibald

[Numpy-discussion] "Manually" broadcasting arrays in Python

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

Has anyone implemented an easier or more efficient way to broadcast arrays to a 
common shape at the Python level? I was hoping that the broadcast iterator would 
actually provide the broadcasted arrays, but it does not.

I've attached my best pure-Python effort.

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
  that is made terrible by our own mad attempt to interpret it as though it had
  an underlying truth."
   -- Umberto Eco

Re: [Numpy-discussion] return type diffences of indexed arrays with Intel C++ compiler (Python 2.5)

[Gael V. <gae...@no...>]

On Tue, Oct 03, 2006 at 11:53:22AM +0200, Lars Bittrich wrote:
> On Monday 02 October 2006 23:53, Travis Oliphant wrote:

> > This is a Python 2.5 issue (the new __index__ method) was incorrectly
> > implemented and allowing a 1-d array to be interpreted as an index.

> > This should be fixed in SVN.

> Ok, thank you. I tried that but the version of python-svn is 2.6

No, Travis was talking of the numpy svn, not the python svn.

Gael

Re: [Numpy-discussion] return type diffences of indexed arrays with Intel C++ compiler (Python 2.5)

[Lars B. <lar...@go...>]

On Monday 02 October 2006 23:53, Travis Oliphant wrote:

> This is a Python 2.5 issue (the new __index__ method) was incorrectly
> implemented and allowing a 1-d array to be interpreted as an index.
>
> This should be fixed in SVN.

Ok, thank you. I tried that but the version of python-svn is 2.6 and there are 
much more complications like:

File "[...]/lib/python2.6/site-packages/numpy/core/numeric.py", line 284
     as = a.shape
      ^
SyntaxError: invalid syntax

File "[...]/lib/python2.6/site-packages/numpy/f2py/crackfortran.py", line 1647
    as=b['args']
     ^
SyntaxError: invalid syntax

I should try Python 2.4 for more stable results.

P.S. Sorry for spaming you. I had problems with my mail address settings. I 
did not expect the mails to arrive after a few days. I thought that they were 
rejected.

Re: [Numpy-discussion] Vectorizing code, for loops, and all that

[A. M. A. <per...@gm...>]

On 02/10/06, Travis Oliphant <oli...@ee...> wrote:

> Perhaps those inner 1-d loops could be optimized (using prefetch or
> something) to reduce the number of cache misses on the inner
> computation, and the concept of looping over the largest dimension
> (instead of the last dimension) should be re-considered.

Cache control seems to be the main factor deciding the speed of many
algorithms. Prefectching could make a huge difference, particularly on
NUMA machines (like a dual opteron). I think GCC has a moderately
portable way to request it (though it may be only in beta versions as
yet).

More generally, all the tricks that ATLAS uses to accelerate BLAS
routines would (in principle) be applicable here. The implementation
would be extremely difficult, though, even if all the basic loops
could be expressed in a few primitives.

A. M. Archibald

Re: [Numpy-discussion] Complex numbers response

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

Kenny Ortmann wrote:
>> There may be a better way, but::
>>
>>    alltrue(isreal(x))
>>
>> Would work. As would:
>>
>>    not sometrue(x.imag)
>>
>> In the above test you are already negating the test, so you could just
>> drop the not.
>>     
>>>  and if so is
>>> there a way to extract the(a + ib) because the absolute value of a 
>>> complex
>>> number is like the pythagorean therom on a and b?
>>>       
>
> what i was looking for in the above is that when I looked up the 
> abs(complexnum) i saw that the abs() of a complex num is (a^2 + b^2)^.5 so i 
> was asking for a way to do this because i was assuming that i was going to 
> have to parse through the array in a for loop and if i found a complex 
> number then get the a and b to do the pythagorean theorom but it seems i 
> will not have to do this, as for the matlab code see below
>   
OK. Yeah, you shouldn't have to do element by element calculations here.

>> I'm not entirely sure what you are looking for here, but x.imag and
>> x.real will give you the real and imaginary parts. abs(x) will return
>> the magnitude of x whether x is real or complex. x.conj() will return
>> complex conjugate.
>>
>> I'm somewhat suspicious of that matlab code. The code given is
>> discontinuous as you cross the negative real axis. Does the result
>> subsequently get squared or something? I'm guessing that either the
>> matlab code is doing extra work, or there are some hidden assumptions
>> (all values are in the positive real half-plane). Or some such. In
>> either case, you'd probably be OK just skipping the check for realness
>> and always taking the absolute of array. I can't say for sure without
>> more context though.
>>
>> -tim
>>     
>
> this is taken out of a program i am converting from matlab to python, so I 
> can not explain all of it, it would take forever, BUT
> COEmix, the input, is equal to a different arrays 0 and 2 rows added 
> together.  So with that being said, it can not be empty, and must be actual 
> numbers.  The reason I am "scared" to remove the isreal statement that i 
> asked about is because the array before hand has values that are taken from 
> square roots.  I do not know if these could ever end up being complex and I 
> am meeting with the head of this project tomorrow to discuss it, but you 
> asked so i thought i would share.
>   

Just to beat a deceased equine a bit more: what I'm concerned about here 
is that the code as it stands will result in a positive number for 
(-x+jy), not matter how small y is, while resulting in a negative number 
for -x. That discontinuity seems unlikely to be what is desired.  
Actually, it's worse than that, if there are any complex numbers in the 
array, all the numbers will be converted to positive real numbers. So, 
for example:

    [-1, 1, -2, 2, 5] => [-1, 1, -2, 2, 5]

while:

    [-1, 1, -2, 2, 3+4j] => [1, 1, 2, 2, 5]

Now it's quite possible that the way COEmix is generated constrains it's 
value in some way. It's possible that it's constrained to the positive 
half plane, or perhaps more likely that it's ideally a positive real 
number that may stray slightly off the real axis, so to clean it up the 
absolute value is being taken.

Since, I can't think of a case where you actually want the discontinuity 
described above, I'm guessing that the code would be just fine as
   COEmix = abs(COEmix)
That is, unconditionally take the absolute value. However, that *is* 
just a guess. The answer lies in how COEmix is originally generated and 
perhaps in how it is used as well.

Regards,

-tim

[... CODE...]

Re: [Numpy-discussion] Complex numbers response

[Kenny O. <kor...@id...>]

> There may be a better way, but::
>
>    alltrue(isreal(x))
>
> Would work. As would:
>
>    not sometrue(x.imag)
>
> In the above test you are already negating the test, so you could just
> drop the not.
>>  and if so is
>> there a way to extract the(a + ib) because the absolute value of a 
>> complex
>> number is like the pythagorean therom on a and b?

what i was looking for in the above is that when I looked up the 
abs(complexnum) i saw that the abs() of a complex num is (a^2 + b^2)^.5 so i 
was asking for a way to do this because i was assuming that i was going to 
have to parse through the array in a for loop and if i found a complex 
number then get the a and b to do the pythagorean theorom but it seems i 
will not have to do this, as for the matlab code see below

> I'm not entirely sure what you are looking for here, but x.imag and
> x.real will give you the real and imaginary parts. abs(x) will return
> the magnitude of x whether x is real or complex. x.conj() will return
> complex conjugate.
>
> I'm somewhat suspicious of that matlab code. The code given is
> discontinuous as you cross the negative real axis. Does the result
> subsequently get squared or something? I'm guessing that either the
> matlab code is doing extra work, or there are some hidden assumptions
> (all values are in the positive real half-plane). Or some such. In
> either case, you'd probably be OK just skipping the check for realness
> and always taking the absolute of array. I can't say for sure without
> more context though.
>
> -tim

this is taken out of a program i am converting from matlab to python, so I 
can not explain all of it, it would take forever, BUT
COEmix, the input, is equal to a different arrays 0 and 2 rows added 
together.  So with that being said, it can not be empty, and must be actual 
numbers.  The reason I am "scared" to remove the isreal statement that i 
asked about is because the array before hand has values that are taken from 
square roots.  I do not know if these could ever end up being complex and I 
am meeting with the head of this project tomorrow to discuss it, but you 
asked so i thought i would share.

function [pospeakind,negpeakind]=peakdetect(COEmix)
if ~nargin
   COEmix=input('enter COEmix vector or return for empty outputs\n');
   if isempty(COEmix)
      pospeakind=[];
      negpeakind=[];
      return
   end
end
sizsig=size(COEmix);
while isempty(COEmix)|~isnumeric(COEmix)|~all(all(isfinite(COEmix)))...
      |length(sizsig)>2|min(sizsig)~=1
   COEmix=input(['COEmix is empty, nonnumeric, nonfinite, or 
nonvector:\nenter '...
         'finite vector or return for empty outputs\n']);
   if isempty(COEmix)
      pospeakind=[];
      negpeakind=[];
      return
   end
   sizsig=size(COEmix);
end
if ~isreal(COEmix)
   COEmix=abs(COEmix);
end
if ~any(COEmix-COEmix(1))
   pospeakind=[];
   negpeakind=[];
   disp('constant COEmix graph suppressed')
   return
end 

Re: [Numpy-discussion] complex numbers

[David L G. <Dav...@no...>]

Thanks Tim,

DG

Tim Hochberg wrote:
> David L Goldsmith wrote:
>   
>> PS:  The Python built in function (i.e., you don't even need numpy for 
>> this) abs(x) is "vectorized" (i.e., accepts a (nested) sequence, incl. 
>> numpy array, argument) and overloaded to give the modulus (i.e., 
>> Pythagorean "length") of a complex number when such is its argument.
>>   
>>     
> This isn't quite right. The built in abs function looks for the special 
> method __abs__. So, abs(x) is equivalent to x.__abs__(). Arrays supply 
> an appropriate __abs__ method, lists do not. For example:
>
>  >>> l = [1,2,3,4]
>  >>> abs(l)
> Traceback (most recent call last):
>   File "<stdin>", line 1, in ?
> TypeError: bad operand type for abs()
>  >>> a = numpy.arange(5)
>  >>> abs(a)
> array([0, 1, 2, 3, 4])
>
>
> -tim
>   
>> DG
>>
>> Kenny Ortmann wrote:
>>   
>>     
>>> excuse my laziness for not looking this up, I googled it but could not find 
>>> a solution.
>>> matlab has a
>>> isreal(array)
>>> where if the array is full of real numbers the value returned is 1.
>>> I'm translating matlab code and ran across
>>>
>>> if ~isreal(array)
>>>     array = abs(array)
>>>
>>> Is there a way to check to see if a number is real or complex? and if so is 
>>> there a way to extract the(a + ib) because the absolute value of a complex 
>>> number is like the pythagorean therom on a and b?
>>>
>>> thanks for your help,
>>> Kenny
>>>
>>>
>>>
>>> -------------------------------------------------------------------------
>>> 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
>>> Num...@li...
>>> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>>>   
>>>     
>>>       
>>   
>>     
>
>
>
> -------------------------------------------------------------------------
> 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
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>   


-- 
HMRD/ORR/NOS/NOAA <http://response.restoration.noaa.gov/emergencyresponse/>

Re: [Numpy-discussion] complex numbers

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

David L Goldsmith wrote:
> PS:  The Python built in function (i.e., you don't even need numpy for 
> this) abs(x) is "vectorized" (i.e., accepts a (nested) sequence, incl. 
> numpy array, argument) and overloaded to give the modulus (i.e., 
> Pythagorean "length") of a complex number when such is its argument.
>   
This isn't quite right. The built in abs function looks for the special 
method __abs__. So, abs(x) is equivalent to x.__abs__(). Arrays supply 
an appropriate __abs__ method, lists do not. For example:

 >>> l = [1,2,3,4]
 >>> abs(l)
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
TypeError: bad operand type for abs()
 >>> a = numpy.arange(5)
 >>> abs(a)
array([0, 1, 2, 3, 4])


-tim
> DG
>
> Kenny Ortmann wrote:
>   
>> excuse my laziness for not looking this up, I googled it but could not find 
>> a solution.
>> matlab has a
>> isreal(array)
>> where if the array is full of real numbers the value returned is 1.
>> I'm translating matlab code and ran across
>>
>> if ~isreal(array)
>>     array = abs(array)
>>
>> Is there a way to check to see if a number is real or complex? and if so is 
>> there a way to extract the(a + ib) because the absolute value of a complex 
>> number is like the pythagorean therom on a and b?
>>
>> thanks for your help,
>> Kenny
>>
>>
>>
>> -------------------------------------------------------------------------
>> 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
>> Num...@li...
>> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>>   
>>     
>
>
>   

Re: [Numpy-discussion] complex numbers

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

Kenny Ortmann wrote:
> excuse my laziness for not looking this up, I googled it but could not find 
> a solution.
> matlab has a
> isreal(array)
> where if the array is full of real numbers the value returned is 1.
> I'm translating matlab code and ran across
>
> if ~isreal(array)
>     array = abs(array)
>
> Is there a way to check to see if a number is real or complex?
There may be a better way, but::

    alltrue(isreal(x))

Would work. As would:

    not sometrue(x.imag)

In the above test you are already negating the test, so you could just 
drop the not.
>  and if so is 
> there a way to extract the(a + ib) because the absolute value of a complex 
> number is like the pythagorean therom on a and b?
>   
I'm not entirely sure what you are looking for here, but x.imag and 
x.real will give you the real and imaginary parts. abs(x) will return 
the magnitude of x whether x is real or complex. x.conj() will return 
complex conjugate.

I'm somewhat suspicious of that matlab code. The code given is 
discontinuous as you cross the negative real axis. Does the result 
subsequently get squared or something? I'm guessing that either the 
matlab code is doing extra work, or there are some hidden assumptions 
(all values are in the positive real half-plane). Or some such. In 
either case, you'd probably be OK just skipping the check for realness 
and always taking the absolute of array. I can't say for sure without 
more context though.

-tim

Re: [Numpy-discussion] complex numbers

[David L G. <Dav...@no...>]

PS:  The Python built in function (i.e., you don't even need numpy for 
this) abs(x) is "vectorized" (i.e., accepts a (nested) sequence, incl. 
numpy array, argument) and overloaded to give the modulus (i.e., 
Pythagorean "length") of a complex number when such is its argument.

DG

Kenny Ortmann wrote:
> excuse my laziness for not looking this up, I googled it but could not find 
> a solution.
> matlab has a
> isreal(array)
> where if the array is full of real numbers the value returned is 1.
> I'm translating matlab code and ran across
>
> if ~isreal(array)
>     array = abs(array)
>
> Is there a way to check to see if a number is real or complex? and if so is 
> there a way to extract the(a + ib) because the absolute value of a complex 
> number is like the pythagorean therom on a and b?
>
> thanks for your help,
> Kenny
>
>
>
> -------------------------------------------------------------------------
> 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
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>   


-- 
HMRD/ORR/NOS/NOAA <http://response.restoration.noaa.gov/emergencyresponse/>

Re: [Numpy-discussion] complex numbers

[David G. <Dav...@no...>]

numpy.isreal(a) is a "top level" function (i.e., not a class member 
function) in numpy; here's its help doc:

Help on function isreal in module numpy.lib.type_check:

isreal(x)
    Return a boolean array where elements are True if that element
    is real (has zero imaginary part)
   
    For scalars, return a boolean.

As for checking for pure imaginaries (i.e., real part == 0), one could use:

    x.imag == -j*x

or

    numpy.isreal(i*x)

or simply

    (x.real == 0.0) and (x.imag != 0.0)

DG

Kenny Ortmann wrote:
> excuse my laziness for not looking this up, I googled it but could not find 
> a solution.
> matlab has a
> isreal(array)
> where if the array is full of real numbers the value returned is 1.
> I'm translating matlab code and ran across
>
> if ~isreal(array)
>     array = abs(array)
>
> Is there a way to check to see if a number is real or complex? and if so is 
> there a way to extract the(a + ib) because the absolute value of a complex 
> number is like the pythagorean therom on a and b?
>
> thanks for your help,
> Kenny
>
>
>
> -------------------------------------------------------------------------
> 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
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>