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

Re: [Numpy-discussion] Biulding Numarray with atlas.

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

Todd Miller wrote:
> I think this is a SciPy FAQ and following their advice is how I've
> worked around it in the past myself:
> 
> http://www.scipy.org/documentation/mailman?fn=scipy-dev/2001-November/000105.html

I'd seem that note before, and it wasn't very helpful, but then I 
noticed, in the Numeric customize.py:

  Example: Some Linux distributions have ATLAS, but not LAPACK (if you're
     liblapack file is smaller than a megabyte, this is probably you)
     Use our own f2c'd lapack libraries, but use ATLAS for BLAS.

My linux distro (FC4) has a full lapack, but it's not optimized, in 
fact, it's  slower than lapack-lite. However, I got atlas from the atlas 
site, and it doesn't come with a full lapack, which would explain the 
missing symbols. Using:

#  Using ATLAS blas in /usr/local with lapack-lite
if 1:
     use_system_lapack = 0
     use_system_blas = 1
     lapack_library_dirs = ['/usr/local/lib/atlas']
     lapack_libraries = ['cblas', 'f77blas', 'atlas', 'g2c']

It compiles, works, and is about twice as fast as using the the 
Numeric-supplied blas.

I still would like to find a full, optimized lapack for Fedora Core 4. 
On my last machine, I got about a six times speed-up using the atlas 
that came with Gentoo.

But for the moment, I'm OK.

Also, I saw something in the atlas docs about how to merge the atlas 
lapack with another full lapack, so maybe I can get that to work.

Now on to see if I can do something similar with numarray.

-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] Biulding Numarray with atlas.

[Todd M. <jm...@st...>]

On Thu, 2005-07-21 at 14:23, Chris Barker wrote:
> Nadav Horesh wrote:
> > Whenever you need a patch to compile numarray with lapack/atlas, you
> > should check if the linear_algebra module works. 
> 
> Of course. And test I have, and no, it doesn't work:
> 
> ImportError:
> /usr/lib/python2.4/site-packages/numarray/linear_algebra/lapack_lite2.so:
> undefined symbol: dgesdd_
> 
> I'm guessing this is a symbol mangling mismatch between what
> lapack_lite2 is expecting and what my atlas install used (i.e. too many
> or not enough underscores) For what it's worth, I get the same error
> with Numeric.

I think this is a SciPy FAQ and following their advice is how I've
worked around it in the past myself:

http://www.scipy.org/documentation/mailman?fn=scipy-dev/2001-November/000105.html

Re: [Numpy-discussion] Biulding Numarray with atlas.

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

Nadav Horesh wrote:
> Whenever you need a patch to compile numarray with lapack/atlas, you
> should check if the linear_algebra module works. 

Of course. And test I have, and no, it doesn't work:

ImportError:
/usr/lib/python2.4/site-packages/numarray/linear_algebra/lapack_lite2.so:
undefined symbol: dgesdd_

I'm guessing this is a symbol mangling mismatch between what
lapack_lite2 is expecting and what my atlas install used (i.e. too many
or not enough underscores) For what it's worth, I get the same error
with Numeric.


which brings me back to my first question:

Where should I get a good lapack/blas for Fedora core 4?

There are lapack and blas rpms, but there is no indication of whether
they are optimized or not.

-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] Biulding Numarray with atlas.

[Todd M. <jm...@st...>]

On Wed, 2005-07-20 at 17:56, Chris Barker wrote:
> Hi all,
> 
> I'm working on building numarray with atlas on a Fedora core 4 system. I 
> installed atlas by downloading a binary from the atlas sourceforge site, 
> and copying stuff into /usr/local/...
> 
> Then in cfc_packages.py, I added:
> 
> USE_LAPACK = True
> 
> at the top. That doesn't seem like where I should do it, but it works.
> 
> then I got:
> 
> Src/_dotblas.c:15:19: error: cblas.h: No such file or directory
> 
> So I looked and saw:
> 
>      else:
>          lapack_dirs = ['/usr/local/lib/atlas']
>          lapack_libs = ['lapack', 'cblas', 'f77blas', 'atlas', 'g2c', 'm']
> 
> which makes sense, we haven't added the include directory that cblas.h 
> is in.
> 
> So I added:
> lapack_include_dirs += ['/usr/local/include/atlas/']
> 
> and now it seems to work.
> 
> Perhaps an include line should be added to the stock cfc_packages.py

Sounds good to me.  I added this.

Thanks,
Todd

[Numpy-discussion] Biulding Numarray with atlas.

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

Hi all,

I'm working on building numarray with atlas on a Fedora core 4 system. I 
installed atlas by downloading a binary from the atlas sourceforge site, 
and copying stuff into /usr/local/...

Then in cfc_packages.py, I added:

USE_LAPACK = True

at the top. That doesn't seem like where I should do it, but it works.

then I got:

Src/_dotblas.c:15:19: error: cblas.h: No such file or directory

So I looked and saw:

     else:
         lapack_dirs = ['/usr/local/lib/atlas']
         lapack_libs = ['lapack', 'cblas', 'f77blas', 'atlas', 'g2c', 'm']

which makes sense, we haven't added the include directory that cblas.h 
is in.

So I added:
lapack_include_dirs += ['/usr/local/include/atlas/']

and now it seems to work.

Perhaps an include line should be added to the stock cfc_packages.py

-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] LAPACK/BLAS for Fedora core 4 i386 ?

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

Darren Dale wrote:
> I can't comment on the Fedora question. With Gentoo, you should not use the 
> "atlas" package, because it comes with an incomplete LAPACK, as you pointed 
> out. 
> 
> Instead, Gentoo-ers should use the "blas-atlas" and "lapack-atlas" packages.

Well, I've moved away from Gentoo, but I  just may go back. However, I'm 
pretty sure I was using the lapack-atlas package, and while it was 
described as "complete", it did not, in fact, have all of lapack. Maybe 
I was doing something wrong, however.

-thanks, Chris

Re: [Numpy-discussion] LAPACK/BLAS for Fedora core 4 i386 ?

[Darren D. <dd...@co...>]

Hi Chris,

On Wednesday 20 July 2005 02:46 pm, Chris Barker wrote:
> Hi all,
>
> Can anyone recommend the best way to get a native LAPACK installed on
> Fedora core 4?
>
> I'm quite surprised that I can't find an atlas rpm with yum.
>
> While I'm at it, when I was using Gentoo, it came with a nice atlas
> blas/lapack. Unfortunately, the atlas lapack does not include all of
> lapack. It has everything Numeric/numarray needs, but not some stuff i
> need for another project (Banded Matrix solvers). Does anyone have a
> suggestion for how to add the addional lapack stuff I need, while still
> using atlas stuff where possible?

I can't comment on the Fedora question. With Gentoo, you should not use the 
"atlas" package, because it comes with an incomplete LAPACK, as you pointed 
out. 

Instead, Gentoo-ers should use the "blas-atlas" and "lapack-atlas" packages.

-- 

Darren

[Numpy-discussion] LAPACK/BLAS for Fedora core 4 i386 ?

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

Hi all,

Can anyone recommend the best way to get a native LAPACK installed on 
Fedora core 4?

I'm quite surprised that I can't find an atlas rpm with yum.

While I'm at it, when I was using Gentoo, it came with a nice atlas 
blas/lapack. Unfortunately, the atlas lapack does not include all of 
lapack. It has everything Numeric/numarray needs, but not some stuff i 
need for another project (Banded Matrix solvers). Does anyone have a 
suggestion for how to add the addional lapack stuff I need, while still 
using atlas stuff where possible?

Frankly, I don't understand why atlas doesn't just include all of 
lapack, using generic versions for anything they haven't optimized, it 
would be a lot easier to get one stop shopping.

-Chris

Re: [Numpy-discussion] unexpected behavior of na.where (no short-circuiting)

[Warren F. <fo...@sl...>]

Well, that's the way Python works - function arguments are evaluated
before the function is called.

Even if there were some way to pass the arguments in as expressions, and
evaluate them on the fly, you'd need some sort of lazy evaluation, and
some pretty deep parsing of the expression - in this case, you could calc
1./a[x,y,z,whatever] for the appropriate values of [x,y,z,whatever], but
what if, instead of 1./a, you have divide.outer(b,c)?  Or
videoCamera.getNextFrame()?  I don't think the general problem of figuring
out how to iterate over the indices of subexpressions in the argument to
get the appropriate index into the evaluated argument is soluble, even
when there are indexable subexpressions.

In this case, you can pre-mask the array:
>>> a = na.arange(4)-2
a[a==0] = 1./999
>>> na.where(a != 0,1./a, 999)

Warren Focke

On Tue, 19 Jul 2005, Sebastian Haase wrote:

> Hi,
> any thoughts on this ? I really think it's counter intuitive ?
>
> Thanks,
> Sebastian Haase
>
> On Tuesday 05 July 2005 17:34, Sebastian Haase wrote:
> > Hi,
> >
> > I was very surprised when I got this warning:
> > >>> a = na.arange(4)-2
> > >>> na.where(a != 0,1./a, 999)
> >
> > Warning: Encountered divide by zero(s)  in divide
> > [  -0.5   -1.   999.     1. ]
> >
> > Then I realized that this is generally referred to as (not) "short
> > circuiting" (e.g. in the case of the '?:'-C-operator, the middle part never
> > gets evaluated at all if the first part evals to 0 )
> >
> > Especially annoying was this because (for debugging) I had set this
> >
> > error-mode:
> > >>> na.Error.setMode(dividebyzero="error")
> >
> > My questions are:
> >  a) Did other people encounter this problem ?
> >  b) What is the general feeling about this actually being a "problem" ?
> >  c) Could this (at all possible) be implemented differently ?
> >
> > Thanks,
> > Sebastian Haase
> >
> >
> >
> >
> >
> >
> >
> > -------------------------------------------------------
> > SF.Net email is sponsored by: Discover Easy Linux Migration Strategies
> > from IBM. Find simple to follow Roadmaps, straightforward articles,
> > informative Webcasts and more! Get everything you need to get up to
> > speed, fast. http://ads.osdn.com/?ad_id=7477&alloc_id=16492&op=click
> > _______________________________________________
> > Numpy-discussion mailing list
> > Num...@li...
> > https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>
>
> -------------------------------------------------------
> SF.Net email is sponsored by: Discover Easy Linux Migration Strategies
> from IBM. Find simple to follow Roadmaps, straightforward articles,
> informative Webcasts and more! Get everything you need to get up to
> speed, fast. http://ads.osdn.com/?ad_id=7477&alloc_id=16492&op=click
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>

Re: [Numpy-discussion] unexpected behavior of na.where (no short-circuiting)

[Sebastian H. <ha...@ms...>]

Hi,
any thoughts on this ? I really think it's counter intuitive ?

Thanks,
Sebastian Haase

On Tuesday 05 July 2005 17:34, Sebastian Haase wrote:
> Hi,
>
> I was very surprised when I got this warning:
> >>> a = na.arange(4)-2
> >>> na.where(a != 0,1./a, 999)
>
> Warning: Encountered divide by zero(s)  in divide
> [  -0.5   -1.   999.     1. ]
>
> Then I realized that this is generally referred to as (not) "short
> circuiting" (e.g. in the case of the '?:'-C-operator, the middle part never
> gets evaluated at all if the first part evals to 0 )
>
> Especially annoying was this because (for debugging) I had set this
>
> error-mode:
> >>> na.Error.setMode(dividebyzero="error")
>
> My questions are:
>  a) Did other people encounter this problem ?
>  b) What is the general feeling about this actually being a "problem" ?
>  c) Could this (at all possible) be implemented differently ?
>
> Thanks,
> Sebastian Haase
>
>
>
>
>
>
>
> -------------------------------------------------------
> SF.Net email is sponsored by: Discover Easy Linux Migration Strategies
> from IBM. Find simple to follow Roadmaps, straightforward articles,
> informative Webcasts and more! Get everything you need to get up to
> speed, fast. http://ads.osdn.com/?ad_id=7477&alloc_id=16492&op=click
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion

Re: [Numpy-discussion] Reversing elements of array

[Robert K. <rk...@uc...>]

Egg wrote:
> Hi all,
> 
> I would like to reverse the elements of an array (or vector).  If it were a
> regular Python list, I could say L.reverse().  If it were in Matlab (ick), I
> could say U(4:-1:1).  Can anyone give me a hint on doing this in Numeric (not
> Numarray)?

In [2]: U = arange(10)

In [3]: U[::-1]
Out[3]: array([9, 8, 7, 6, 5, 4, 3, 2, 1, 0])

-- 
Robert Kern
rk...@uc...

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

Re: [Numpy-discussion] Reversing elements of array

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

A Monday 18 July 2005 20:00, Egg va escriure:
> I would like to reverse the elements of an array (or vector).  If it were=
 a
> regular Python list, I could say L.reverse().  If it were in Matlab (ick),
> I could say U(4:-1:1).  Can anyone give me a hint on doing this in Numeric
> (not Numarray)?

>>> import Numeric
>>> a=3DNumeric.arange(10)
>>> a[::-1]
array([9, 8, 7, 6, 5, 4, 3, 2, 1, 0])

(the same works for python lists)

Cheers,

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

[Numpy-discussion] Reversing elements of array

[Egg <web...@gm...>]

Hi all,

I would like to reverse the elements of an array (or vector).  If it were a
regular Python list, I could say L.reverse().  If it were in Matlab (ick), I
could say U(4:-1:1).  Can anyone give me a hint on doing this in Numeric (not
Numarray)?

Thx

Re: [Numpy-discussion] optimizations

[Thomas G. <gr...@gr...>]

Ok,
i'm answering myself partly since no-one else did

> I realized that the ufuncs in numarrays are not as fast as possible 
> because they are coded very straightforward. Is there any particular 
> reason why this is the case, like cleanness of code, confidence in the 
> compiler or because the code was automatically generated?
>

sorry for disregarding the manual in this respect... clearly the 
codelets are generated

> I would like to contribute assembly SIMD codelets (SSE and Altivec), i 
> have been successfully using in my projects for quite some time. Is it 
> feasible to submit patches, so that these go into the main numarray 
> distribution, or should i rather implement this as separate ufuncs?
>

well then, is it an option to provide hand-written ufunc 
implementations, circumventing the code generation?
Is the function nomenclature stable enough that this would work for a 
longer time?

thanks,
Thomas

Re: [Numpy-discussion] problem with poisson generators

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

Hi,
This seems more and more like some coding problem.  The fact that
everything you have tried seems to work especially replacing the
ranlib generator with R's (using the standalone library  just avoids
the calls to R).  I am not an expert in this but it seems that your
code may not be correctly calling things correctly.

Perhaps you should find the types of the variables like Lpos_esp as
this really should be an C int.

Also, if it remains a problem I wouldn't mind seeing the values and
types of the variables  used to compute Lpos_esp.

Regards
Bruce

On 7/15/05, Flavio Coelho <fcc...@gm...> wrote:
>=20
>=20
> 2005/7/15, Todd Miller <jm...@st...>:
>=20
>  Todd,
> =20
>  the  function  within which the problem happens runs fine by itself, as
> well as the simplified version you sent me.
> =20
>  I am running the code on a AMD64 (though all my software and OS is compi=
led
> to 686 - Gentoo Linux). Is there any issue I should be aware of,  regardi=
ng
> this specific cpu? I can try to run the software on another machine...
> =20
>=20
> > Having this code is a good start to solving the problem.  I think the
> > next step is to simplify your example to make it runnable and provide=
=20
> > known inputs for all the parameters which lead to a failure for you.
> >=20
> > Being really literal (spelling out every single damn thing) cuts down o=
n
> > speculation and inefficiency on our end.
> >=20
> > It would also be good to express the condition (or it's inverse) you=20
> > think is an error as an assertion,  so something like this might be wha=
t
> > we need:
> >=20
> > from numarray.random_array import poisson
> >=20
> > E, I, S =3D (0,0,0)
> > beta,alpha,e,r,delta,B,w,p =3D (0.001,1,1,0.5,1,0,0,0)
> > theta =3D 0
> > npass =3D 0
> > N =3D 100 # total guess here for me
> > Lpos_esp =3D beta*S*((I+theta)/(N+npass))**alpha #Number of
> new cases
> > Lpos =3D poisson(Lpos_esp)
> > assert not (theta =3D=3D 0 and Lpos_esp =3D=3D 0 and Lpos > 0)
> >=20
> > The problem is,  the above works for me.  Make it sane and get it to
> > expose the error for you and I'll look into this some more.
> >=20
> > Regards,
> > Todd
> >=20
> > >
> > > I wonder if called by itself  it would trigger the problem... The=20
> > > commented Lines Is where I catch the errors: when poisson returns a
> > > greater than zero number, when called with mean 0.
> > >
> > >
> > >
> > >         Ranlib uses static floats so it may relate to numerical=20
> > >         precision (as
> > >         well as the random uniform random number generator). Really
> > >         the only
> > >         way is for you to debug the ranlib.c file
> > >
> > > I'll see what I can do.=20
> > >
> > >
> > >         Note that R provides a standalone math library in C that
> > >         includes the
> > >         random number generators so you could you those instead. SWIG
> > >         handles
> > >         it rather well.
> > >
> > >
> > > I think thats what Rpy already does...is it not?
> > >
> > > thanks Bruce,
> > >
> > > Fl=E1vio
> > >
> > >
> > >         Regards
> > >         Bruce
> > >=20
> > >         On 7/13/05, Flavio Coelho <fcc...@gm...> wrote:
> > >         > Well,
> > >         >   I am definetly glad that someone has also stumbled onto=
=20
> > >         the same problem.
> > >         >
> > >         >  But it is nevertheless intriguing, that you can run poisso=
n
> > >         a million times
> > >         > with mean  zero or negative(it assumes zero mean inthis=20
> > >         case) without any
> > >         > problems by itself. But when I call it within my code, the
> > >         rate of error is
> > >         > very high (I would say it returns a wrong result every time=
,=20
> > >         but I haven't
> > >         > checked every answer).
> > >         >
> > >         >  Meanwhile, my solution will be:
> > >         >
> > >         >  import rpy
> > >         >=20
> > >         >  n =3D rpy.r.rpois(n,mean)
> > >         >
> > >         >  I don't feel I can trust poisson  while this "funny"
> > >         behavior is still
> > >         > there...
> > >         >  If someone has any Idea of how I could trace this bug
> > >         please tell me, and
> > >         > I'll hunt it down. At least I can reproduce it in a very=20
> > >         specific context.
> > >         >
> > >         >  thanks,
> > >         >
> > >         >  Fl=E1vio
> > >         >
> > >         > 2005/7/12, Sebastian Haase < ha...@ms...>:
> > >         > > Hi Flavio!
> > >         > > I had reported this long time ago and this list (about
> > >         numarray).
> > >         > > Somehow this got more or less dismissed. If I recall=20
> > >         correctly the
> > >         > argument
> > >         > > was that nobody could reproduce it (I ran this on Debian
> > >         Woody , py2.2,
> > >         > (with
> > >         > > CVS numarray at the time).=20
> > >         > >
> > >         > > I ended up writting my own wrapper(s):
> > >         > > def poissonArr(shape=3Ddefshape, mean=3D1):
> > >         > >         from numarray import random_array as ra
> > >         > >         if mean =3D=3D 0:
> > >         > >                 return zeroArrF(shape)
> > >         > >         elif mean < 0:
> > >         > >                 raise "poisson not defined for mean < 0"
> > >         > >         else:
> > >         > >                 return ra.poisson(mean, shape).astype
> > >         (na.UInt16)
> > >         > >
> > >         > > def poissonize(arr):
> > >         > >         from numarray import random_array as ra
> > >         > >         return na.where(arr<=3D0, 0, ra.poisson(arr)).ast=
ype
> > >         ( na.UInt16)
> > >         > >
> > >         > > (I use the astype( na.UInt16) because of some OpenGL code=
)
> > >         > >
> > >         > > Just last week had this problem on a windows98 computer=
=20
> > >         (python2.4).
> > >         > >
> > >         > > This should get sorted out ...
> > >         > >
> > >         > > Thanks for reporting this problem.
> > >         > > Sebastian Haase=20
> > >         > >
> > >         > >
> > >         > >
> > >         > >
> > >         > > On Tuesday 12 July 2005 13:32, Flavio Coelho wrote:
> > >         > > > Hi,
> > >         > > >
> > >         > > > I am having problems with the poisson random number
> > >         generators of both
> > >         > > > Numarray and Numeric.
> > >         > > > I can't replicate it when calling the function from the=
=20
> > >         python cosonle,
> > >         > but
> > >         > > > it has consistently malfunctioned when used within one
> > >         of my scripts.
> > >         > > >
> > >         > > > What happens is that it frequently return a value=20
> > >         greater than zero when
> > >         > > > called with zero mean: poisson( 0.0)
> > >         > > >
> > >         > > > Unfortunately My program is too big to send attached bu=
t=20
> > >         I have
> > >         > confirmed
> > >         > > > the malfunction by printing both the mean and the resul=
t
> > >         whenever it
> > >         > spits
> > >         > > > out a wrong result.=20
> > >         > > >
> > >         > > > This is very weird indeed, I have run poisson millions
> > >         of times by itsel
> > >         > on
> > >         > > > the python console, without any problems...=20
> > >         > > >
> > >         > > > I hope it is some stupid mistake, but when I replace th=
e
> > >         poisson
> > >         > function
> > >         > > > call within my program by the R equivalent command=20
> > >         (rpois) via the rpy
> > >         > > > wrapper, everything works just fine...
> > >         > > >
> > >         > > > any Ideas?
> > >         > >
> > >         >=20
> > >         >
> > >         >
> > >         > --
> > >         >
> > >         > Fl=E1vio Code=E7o Coelho
> > >         > registered Linux user # 386432
> > >         > ---------------------------=20
> > >         > Great spirits have always found violent opposition from
> > >         mediocrities. The
> > >         > latter cannot understand it when a man does not
> > >         thoughtlessly submit to
> > >         > hereditary prejudices but honestly and courageously uses hi=
s
> > >         intelligence.
> > >         >         Albert Einstein
> > >         >
> > >
> > >
> > >
> > > --
> > > Fl=E1vio Code=E7o Coelho=20
> > > registered Linux user # 386432
> > > ---------------------------
> > > Great spirits have always found violent opposition from mediocrities.
> > > The
> > > latter cannot understand it when a man does not thoughtlessly submit=
=20
> > > to
> > > hereditary prejudices but honestly and courageously uses his
> > > intelligence.
> > >         Albert Einstein
> >=20
> >=20
>=20
>=20
>=20
> --=20
>=20
> Fl=E1vio Code=E7o Coelho
> registered Linux user # 386432=20
> ---------------------------
> Great spirits have always found violent opposition from mediocrities. The
> latter cannot understand it when a man does not thoughtlessly submit to
> hereditary prejudices but honestly and courageously uses his intelligence=
.=20
>         Albert Einstein=20
>

Re: [Numpy-discussion] problem with poisson generators

[Todd M. <jm...@st...>]

On Fri, 2005-07-15 at 21:30 -0300, Flavio Coelho wrote:
> 
> 
> 2005/7/15, Todd Miller <jm...@st...>:
> 
> Todd,
> 
> the  function  within which the problem happens runs fine by itself,
> as well as the simplified version you sent me.

Maybe you could stick in an assert rather than the print and then use
pdb.pm() to find the actual parameters leading to the failure.  (The
parameters may or may not matter, i.e. the problem may be hidden state
somewhere your program accumulates in some complicated way,  but maybe
we'll get lucky...)

> I am running the code on a AMD64 (though all my software and OS is
> compiled to 686 - Gentoo Linux). Is there any issue I should be aware
> of,  regarding this specific cpu? 

Other than the fact that 64-bit addressing functionality isn't where we
want it to be,  there are no remaining issues I'm aware of for AMD64.
But the addressing is currently hamstrung by 32-bit limits in Python's
sequence protocol so you should be aware that even with your AMD64 you
may run out of head room.

> I can try to run the software on another machine...

I wouldn't bother:  our goal should be to replicate the problem simply
and that suggests to me maintaining identical initial conditions.  BTW,
I'm using AMD64 too so that won't be another variable.

Regards,
Todd

> 
>         Having this code is a good start to solving the problem.  I
>         think the
>         next step is to simplify your example to make it runnable and
>         provide 
>         known inputs for all the parameters which lead to a failure
>         for you.
>         
>         Being really literal (spelling out every single damn thing)
>         cuts down on
>         speculation and inefficiency on our end.
>         
>         It would also be good to express the condition (or it's
>         inverse) you 
>         think is an error as an assertion,  so something like this
>         might be what
>         we need:
>         
>         from numarray.random_array import poisson
>         
>         E, I, S = (0,0,0)
>         beta,alpha,e,r,delta,B,w,p = (0.001,1,1,0.5,1,0,0,0)
>         theta = 0
>         npass = 0
>         N = 100 # total guess here for me
>         Lpos_esp = beta*S*((I+theta)/(N+npass))**alpha #Number of new
>         cases
>         Lpos = poisson(Lpos_esp)
>         assert not (theta == 0 and Lpos_esp == 0 and Lpos > 0)
>         
>         The problem is,  the above works for me.  Make it sane and get
>         it to
>         expose the error for you and I'll look into this some more.
>         
>         Regards,
>         Todd
>         
>         >
>         > I wonder if called by itself  it would trigger the
>         problem... The 
>         > commented Lines Is where I catch the errors: when poisson
>         returns a
>         > greater than zero number, when called with mean 0.
>         >
>         >
>         >
>         >         Ranlib uses static floats so it may relate to
>         numerical 
>         >         precision (as
>         >         well as the random uniform random number generator).
>         Really
>         >         the only
>         >         way is for you to debug the ranlib.c file
>         >
>         > I'll see what I can do. 
>         >
>         >
>         >         Note that R provides a standalone math library in C
>         that
>         >         includes the
>         >         random number generators so you could you those
>         instead. SWIG
>         >         handles
>         >         it rather well.
>         >
>         >
>         > I think thats what Rpy already does...is it not?
>         >
>         > thanks Bruce,
>         >
>         > Flávio
>         >
>         >
>         >         Regards
>         >         Bruce
>         > 
>         >         On 7/13/05, Flavio Coelho <fcc...@gm...>
>         wrote:
>         >         > Well,
>         >         >   I am definetly glad that someone has also
>         stumbled onto 
>         >         the same problem.
>         >         >
>         >         >  But it is nevertheless intriguing, that you can
>         run poisson
>         >         a million times
>         >         > with mean  zero or negative(it assumes zero mean
>         inthis 
>         >         case) without any
>         >         > problems by itself. But when I call it within my
>         code, the
>         >         rate of error is
>         >         > very high (I would say it returns a wrong result
>         every time, 
>         >         but I haven't
>         >         > checked every answer).
>         >         >
>         >         >  Meanwhile, my solution will be:
>         >         >
>         >         >  import rpy
>         >         > 
>         >         >  n = rpy.r.rpois(n,mean)
>         >         >
>         >         >  I don't feel I can trust poisson  while this
>         "funny"
>         >         behavior is still
>         >         > there...
>         >         >  If someone has any Idea of how I could trace this
>         bug
>         >         please tell me, and
>         >         > I'll hunt it down. At least I can reproduce it in
>         a very 
>         >         specific context.
>         >         >
>         >         >  thanks,
>         >         >
>         >         >  Flávio
>         >         >
>         >         > 2005/7/12, Sebastian Haase <ha...@ms...>:
>         >         > > Hi Flavio!
>         >         > > I had reported this long time ago and this list
>         (about
>         >         numarray).
>         >         > > Somehow this got more or less dismissed. If I
>         recall 
>         >         correctly the
>         >         > argument
>         >         > > was that nobody could reproduce it (I ran this
>         on Debian
>         >         Woody , py2.2,
>         >         > (with
>         >         > > CVS numarray at the time). 
>         >         > >
>         >         > > I ended up writting my own wrapper(s):
>         >         > > def poissonArr(shape=defshape, mean=1):
>         >         > >         from numarray import random_array as ra
>         >         > >         if mean == 0:
>         >         > >                 return zeroArrF(shape)
>         >         > >         elif mean < 0:
>         >         > >                 raise "poisson not defined for
>         mean < 0"
>         >         > >         else:
>         >         > >                 return ra.poisson(mean,
>         shape).astype
>         >         (na.UInt16)
>         >         > >
>         >         > > def poissonize(arr):
>         >         > >         from numarray import random_array as ra
>         >         > >         return na.where(arr<=0, 0, ra.poisson
>         (arr)).astype
>         >         ( na.UInt16)
>         >         > >
>         >         > > (I use the astype( na.UInt16) because of some
>         OpenGL code)
>         >         > >
>         >         > > Just last week had this problem on a windows98
>         computer 
>         >         (python2.4).
>         >         > >
>         >         > > This should get sorted out ...
>         >         > >
>         >         > > Thanks for reporting this problem.
>         >         > > Sebastian Haase 
>         >         > >
>         >         > >
>         >         > >
>         >         > >
>         >         > > On Tuesday 12 July 2005 13:32, Flavio Coelho
>         wrote:
>         >         > > > Hi,
>         >         > > >
>         >         > > > I am having problems with the poisson random
>         number
>         >         generators of both
>         >         > > > Numarray and Numeric.
>         >         > > > I can't replicate it when calling the function
>         from the 
>         >         python cosonle,
>         >         > but
>         >         > > > it has consistently malfunctioned when used
>         within one
>         >         of my scripts.
>         >         > > >
>         >         > > > What happens is that it frequently return a
>         value 
>         >         greater than zero when
>         >         > > > called with zero mean: poisson( 0.0)
>         >         > > >
>         >         > > > Unfortunately My program is too big to send
>         attached but 
>         >         I have
>         >         > confirmed
>         >         > > > the malfunction by printing both the mean and
>         the result
>         >         whenever it
>         >         > spits
>         >         > > > out a wrong result. 
>         >         > > >
>         >         > > > This is very weird indeed, I have run poisson
>         millions
>         >         of times by itsel
>         >         > on
>         >         > > > the python console, without any problems... 
>         >         > > >
>         >         > > > I hope it is some stupid mistake, but when I
>         replace the
>         >         poisson
>         >         > function
>         >         > > > call within my program by the R equivalent
>         command 
>         >         (rpois) via the rpy
>         >         > > > wrapper, everything works just fine...
>         >         > > >
>         >         > > > any Ideas?
>         >         > >
>         >         > 
>         >         >
>         >         >
>         >         > --
>         >         >
>         >         > Flávio Codeço Coelho
>         >         > registered Linux user # 386432
>         >         > --------------------------- 
>         >         > Great spirits have always found violent opposition
>         from
>         >         mediocrities. The
>         >         > latter cannot understand it when a man does not
>         >         thoughtlessly submit to
>         >         > hereditary prejudices but honestly and
>         courageously uses his
>         >         intelligence.
>         >         >         Albert Einstein
>         >         >
>         >
>         >
>         >
>         > --
>         > Flávio Codeço Coelho 
>         > registered Linux user # 386432
>         > ---------------------------
>         > Great spirits have always found violent opposition from
>         mediocrities.
>         > The
>         > latter cannot understand it when a man does not
>         thoughtlessly submit 
>         > to
>         > hereditary prejudices but honestly and courageously uses his
>         > intelligence.
>         >         Albert Einstein
>         
> 
> 
> 
> -- 
> Flávio Codeço Coelho
> registered Linux user # 386432 
> ---------------------------
> Great spirits have always found violent opposition from mediocrities.
> The
> latter cannot understand it when a man does not thoughtlessly submit
> to
> hereditary prejudices but honestly and courageously uses his
> intelligence. 
>         Albert Einstein 

Re: [Numpy-discussion] problem with poisson generators

[Flavio C. <fcc...@gm...>]

2005/7/15, Todd Miller <jm...@st...>:
Todd,

the function within which the problem happens runs fine by itself, as well=
=20
as the simplified version you sent me.

I am running the code on a AMD64 (though all my software and OS is compiled=
=20
to 686 - Gentoo Linux). Is there any issue I should be aware of, regarding=
=20
this specific cpu? I can try to run the software on another machine...

Having this code is a good start to solving the problem. I think the
> next step is to simplify your example to make it runnable and provide
> known inputs for all the parameters which lead to a failure for you.
>=20
> Being really literal (spelling out every single damn thing) cuts down on
> speculation and inefficiency on our end.
>=20
> It would also be good to express the condition (or it's inverse) you
> think is an error as an assertion, so something like this might be what
> we need:
>=20
> from numarray.random_array import poisson
>=20
> E, I, S =3D (0,0,0)
> beta,alpha,e,r,delta,B,w,p =3D (0.001,1,1,0.5,1,0,0,0)
> theta =3D 0
> npass =3D 0
> N =3D 100 # total guess here for me
> Lpos_esp =3D beta*S*((I+theta)/(N+npass))**alpha #Number of new cases
> Lpos =3D poisson(Lpos_esp)
> assert not (theta =3D=3D 0 and Lpos_esp =3D=3D 0 and Lpos > 0)
>=20
> The problem is, the above works for me. Make it sane and get it to
> expose the error for you and I'll look into this some more.
>=20
> Regards,
> Todd
>=20
> >
> > I wonder if called by itself it would trigger the problem... The
> > commented Lines Is where I catch the errors: when poisson returns a
> > greater than zero number, when called with mean 0.
> >
> >
> >
> > Ranlib uses static floats so it may relate to numerical
> > precision (as
> > well as the random uniform random number generator). Really
> > the only
> > way is for you to debug the ranlib.c file
> >
> > I'll see what I can do.
> >
> >
> > Note that R provides a standalone math library in C that
> > includes the
> > random number generators so you could you those instead. SWIG
> > handles
> > it rather well.
> >
> >
> > I think thats what Rpy already does...is it not?
> >
> > thanks Bruce,
> >
> > Fl=E1vio
> >
> >
> > Regards
> > Bruce
> >
> > On 7/13/05, Flavio Coelho <fcc...@gm...> wrote:
> > > Well,
> > > I am definetly glad that someone has also stumbled onto
> > the same problem.
> > >
> > > But it is nevertheless intriguing, that you can run poisson
> > a million times
> > > with mean zero or negative(it assumes zero mean inthis
> > case) without any
> > > problems by itself. But when I call it within my code, the
> > rate of error is
> > > very high (I would say it returns a wrong result every time,
> > but I haven't
> > > checked every answer).
> > >
> > > Meanwhile, my solution will be:
> > >
> > > import rpy
> > >
> > > n =3D rpy.r.rpois(n,mean)
> > >
> > > I don't feel I can trust poisson while this "funny"
> > behavior is still
> > > there...
> > > If someone has any Idea of how I could trace this bug
> > please tell me, and
> > > I'll hunt it down. At least I can reproduce it in a very
> > specific context.
> > >
> > > thanks,
> > >
> > > Fl=E1vio
> > >
> > > 2005/7/12, Sebastian Haase <ha...@ms...>:
> > > > Hi Flavio!
> > > > I had reported this long time ago and this list (about
> > numarray).
> > > > Somehow this got more or less dismissed. If I recall
> > correctly the
> > > argument
> > > > was that nobody could reproduce it (I ran this on Debian
> > Woody , py2.2,
> > > (with
> > > > CVS numarray at the time).
> > > >
> > > > I ended up writting my own wrapper(s):
> > > > def poissonArr(shape=3Ddefshape, mean=3D1):
> > > > from numarray import random_array as ra
> > > > if mean =3D=3D 0:
> > > > return zeroArrF(shape)
> > > > elif mean < 0:
> > > > raise "poisson not defined for mean < 0"
> > > > else:
> > > > return ra.poisson(mean, shape).astype
> > (na.UInt16)
> > > >
> > > > def poissonize(arr):
> > > > from numarray import random_array as ra
> > > > return na.where(arr<=3D0, 0, ra.poisson(arr)).astype
> > ( na.UInt16)
> > > >
> > > > (I use the astype( na.UInt16) because of some OpenGL code)
> > > >
> > > > Just last week had this problem on a windows98 computer
> > (python2.4).
> > > >
> > > > This should get sorted out ...
> > > >
> > > > Thanks for reporting this problem.
> > > > Sebastian Haase
> > > >
> > > >
> > > >
> > > >
> > > > On Tuesday 12 July 2005 13:32, Flavio Coelho wrote:
> > > > > Hi,
> > > > >
> > > > > I am having problems with the poisson random number
> > generators of both
> > > > > Numarray and Numeric.
> > > > > I can't replicate it when calling the function from the
> > python cosonle,
> > > but
> > > > > it has consistently malfunctioned when used within one
> > of my scripts.
> > > > >
> > > > > What happens is that it frequently return a value
> > greater than zero when
> > > > > called with zero mean: poisson( 0.0)
> > > > >
> > > > > Unfortunately My program is too big to send attached but
> > I have
> > > confirmed
> > > > > the malfunction by printing both the mean and the result
> > whenever it
> > > spits
> > > > > out a wrong result.
> > > > >
> > > > > This is very weird indeed, I have run poisson millions
> > of times by itsel
> > > on
> > > > > the python console, without any problems...
> > > > >
> > > > > I hope it is some stupid mistake, but when I replace the
> > poisson
> > > function
> > > > > call within my program by the R equivalent command
> > (rpois) via the rpy
> > > > > wrapper, everything works just fine...
> > > > >
> > > > > any Ideas?
> > > >
> > >
> > >
> > >
> > > --
> > >
> > > Fl=E1vio Code=E7o Coelho
> > > registered Linux user # 386432
> > > ---------------------------
> > > Great spirits have always found violent opposition from
> > mediocrities. The
> > > latter cannot understand it when a man does not
> > thoughtlessly submit to
> > > hereditary prejudices but honestly and courageously uses his
> > intelligence.
> > > Albert Einstein
> > >
> >
> >
> >
> > --
> > Fl=E1vio Code=E7o Coelho
> > registered Linux user # 386432
> > ---------------------------
> > Great spirits have always found violent opposition from mediocrities.
> > The
> > latter cannot understand it when a man does not thoughtlessly submit
> > to
> > hereditary prejudices but honestly and courageously uses his
> > intelligence.
> > Albert Einstein
>=20
>=20


--=20
Fl=E1vio Code=E7o Coelho
registered Linux user # 386432
---------------------------
Great spirits have always found violent opposition from mediocrities. The
latter cannot understand it when a man does not thoughtlessly submit to
hereditary prejudices but honestly and courageously uses his intelligence.
Albert Einstein

RE: [Numpy-discussion] problem with poisson generators

[<Seb...@el...>]

Hi Flavio,

Could you give us every call to poisson or negative_binomial (ie all functi=
ons related to random numbers) preceeded with the seed ?

Adding before your declaration of stepSEIR_s code like:

randoutput =3D file("randoutput.py", "w")
old_poisson =3D poisson
def poisson(m):
    print >> randoutput, "print get_seed(), poisson(%s)"%m
    result =3D old_poisson(m)
    print >> randoutput, "# result =3D %s"%result
    return result

old_negative_binomial =3D negative_binomial
def negative_binomial(i,p):
    print >> randoutput, "print get_seed(), negative_binomial(%s,%s)"%(i,p)
    result =3D old_negative_binomial(i,p)
    print >> randoutput, "# result =3D %s"%result
    return result


will ouput every call.
We can then easily try on our machines what it gives.

Best,

Sebastien

> -----Original Message-----
> From: num...@li...
> [mailto:num...@li...]On Behalf Of Todd
> Miller
> Sent: vendredi 15 juillet 2005 12:54
> To: Flavio Coelho
> Cc: Bruce Southey; Sebastian Haase; numpy-discussion
> Subject: Re: [Numpy-discussion] problem with poisson generators
>=20
>=20
> On Wed, 2005-07-13 at 12:13 -0300, Flavio Coelho wrote:
> >=20
> >=20
> > 2005/7/13, Bruce Southey <bso...@gm...>:
> >         Hi,
> >         What is the seed used?
> >=20
> > I am not setting the seed.
> > =20
> >=20
> >         It is really important that you set the seed?
> >=20
> > No.
> >=20
> >=20
> >         Did you build Python and numarray from source?
> >=20
> >=20
> > Yes,  I  use Gentoo. I  build everything from source.=20
> >=20
> >=20
> >         Can you reduce your code to a few lines that have=20
> the problem?
> >=20
> >=20
> > Not really, poisson and binomial are the only two functions from
> > Numeric that I use but they are called inside a rather=20
> complex object
> > oriented code environment (objects within objetcs, being called
> > recursively...) Bellow is the function within which poisson=20
> is called:
> >=20
> > def=20
> stepSEIR_s(self,inits=3D(0,0,0),par=3D(0.001,1,1,0.5,1,0,0,0),theta=3D0,
> > npass=3D0,dist=3D'poisson'):
> >         """
> >         Defines an stochastic model SEIR:
> >         - inits =3D (E,I,S)
> >         - par =3D (Beta, alpha, E,r,delta,B,w,p) see docs.
> >         - theta =3D infectious individuals from neighbor sites
> >         """
> >         E,I,S =3D inits
> >         N =3D self.parentSite.totpop
> >         beta,alpha,e,r,delta,B,w,p =3D par
> >         Lpos_esp =3D beta*S*((I+theta)/(N+npass))**alpha=20
> #Number of new
> > cases
> >        =20
> >         if dist =3D=3D 'poisson':
> >             Lpos =3D poisson(Lpos_esp)
> > ##            if theta =3D=3D 0 and Lpos_esp =3D=3D 0 and Lpos > 0: =20
> > ##                print
> > Lpos,Lpos_esp,S,I,theta,N,self.parentSite.sitename
> >         elif dist =3D=3D 'negbin':
> >             prob =3D I/(I+Lpos_esp) #convertin between=20
> parameterizations
> >             Lpos =3D negative_binomial(I,prob)
> >         self.parentSite.totalcases +=3D Lpos #update number=20
> of cases   =20
> >         Epos =3D (1-e)*E + Lpos
> >         Ipos =3D e*E + (1-r)*I
> >         Spos =3D S + B - Lpos=20
> >         Rpos =3D N-(Spos+Epos+Ipos)
> >         #self.parentSite.totpop =3D Spos+Epos+Ipos+Rpos
> >         self.parentSite.incidence.append(Lpos)
> >         if not self.parentSite.infected:
> >             if Lpos > 0:
> >                 self.parentSite.infected =3D
> > self.parentSite.parentGraph.simstep
> >                 self.parentSite.parentGraph.epipath.append
> >=20
> ((self.parentSite.parentGraph.simstep,self.parentSite,self.par
> entSite.infector))
> >=20
> >        =20
> >         return [Epos,Ipos,Spos]
> >=20
>=20
> Having this code is a good start to solving the problem.  I think the
> next step is to simplify your example to make it runnable and provide
> known inputs for all the parameters which lead to a failure for you.=20
>=20
> Being really literal (spelling out every single damn thing)=20
> cuts down on
> speculation and inefficiency on our end.
>=20
> It would also be good to express the condition (or it's inverse) you
> think is an error as an assertion,  so something like this=20
> might be what
> we need:
>=20
> from numarray.random_array import poisson
>=20
> E, I, S =3D (0,0,0)
> beta,alpha,e,r,delta,B,w,p =3D (0.001,1,1,0.5,1,0,0,0)
> theta =3D 0
> npass =3D 0
> N =3D 100 # total guess here for me
> Lpos_esp =3D beta*S*((I+theta)/(N+npass))**alpha #Number of new cases
> Lpos =3D poisson(Lpos_esp)
> assert not (theta =3D=3D 0 and Lpos_esp =3D=3D 0 and Lpos > 0)
>=20
> The problem is,  the above works for me.  Make it sane and get it to
> expose the error for you and I'll look into this some more.
>=20
> Regards,
> Todd
>=20
> >=20
> > I wonder if called by itself  it would trigger the problem... The
> > commented Lines Is where I catch the errors: when poisson returns a
> > greater than zero number, when called with mean 0.
> >=20
> >=20
> >=20
> >         Ranlib uses static floats so it may relate to numerical
> >         precision (as
> >         well as the random uniform random number generator). Really
> >         the only=20
> >         way is for you to debug the ranlib.c file
> >=20
> > I'll see what I can do.=20
> >=20
> >=20
> >         Note that R provides a standalone math library in C that
> >         includes the
> >         random number generators so you could you those=20
> instead. SWIG
> >         handles=20
> >         it rather well.
> >=20
> >=20
> > I think thats what Rpy already does...is it not?=20
> >=20
> > thanks Bruce,
> >=20
> > Fl=C3=A1vio
> > =20
> >=20
> >         Regards
> >         Bruce
> >        =20
> >         On 7/13/05, Flavio Coelho <fcc...@gm...> wrote:
> >         > Well,
> >         >   I am definetly glad that someone has also stumbled onto
> >         the same problem.
> >         >
> >         >  But it is nevertheless intriguing, that you can=20
> run poisson
> >         a million times=20
> >         > with mean  zero or negative(it assumes zero mean inthis
> >         case) without any
> >         > problems by itself. But when I call it within my code, the
> >         rate of error is
> >         > very high (I would say it returns a wrong result=20
> every time,
> >         but I haven't=20
> >         > checked every answer).
> >         >
> >         >  Meanwhile, my solution will be:
> >         >
> >         >  import rpy
> >         >
> >         >  n =3D rpy.r.rpois(n,mean)
> >         >
> >         >  I don't feel I can trust poisson  while this "funny"
> >         behavior is still=20
> >         > there...
> >         >  If someone has any Idea of how I could trace this bug
> >         please tell me, and
> >         > I'll hunt it down. At least I can reproduce it in a very
> >         specific context.
> >         >
> >         >  thanks,
> >         >
> >         >  Fl=C3=A1vio=20
> >         >
> >         > 2005/7/12, Sebastian Haase <ha...@ms...>:
> >         > > Hi Flavio!
> >         > > I had reported this long time ago and this list (about
> >         numarray).
> >         > > Somehow this got more or less dismissed. If I recall
> >         correctly the=20
> >         > argument
> >         > > was that nobody could reproduce it (I ran this on Debian
> >         Woody , py2.2,
> >         > (with
> >         > > CVS numarray at the time).
> >         > >
> >         > > I ended up writting my own wrapper(s):=20
> >         > > def poissonArr(shape=3Ddefshape, mean=3D1):
> >         > >         from numarray import random_array as ra
> >         > >         if mean =3D=3D 0:
> >         > >                 return zeroArrF(shape)
> >         > >         elif mean < 0:=20
> >         > >                 raise "poisson not defined for mean < 0"
> >         > >         else:
> >         > >                 return ra.poisson(mean, shape).astype
> >         (na.UInt16)
> >         > >
> >         > > def poissonize(arr):
> >         > >         from numarray import random_array as ra
> >         > >         return na.where(arr<=3D0, 0,=20
> ra.poisson(arr)).astype
> >         ( na.UInt16)
> >         > >
> >         > > (I use the astype( na.UInt16) because of some=20
> OpenGL code)
> >         > >
> >         > > Just last week had this problem on a windows98 computer
> >         (python2.4).
> >         > >
> >         > > This should get sorted out ...=20
> >         > >
> >         > > Thanks for reporting this problem.
> >         > > Sebastian Haase
> >         > >
> >         > >
> >         > >
> >         > >
> >         > > On Tuesday 12 July 2005 13:32, Flavio Coelho wrote:
> >         > > > Hi,=20
> >         > > >
> >         > > > I am having problems with the poisson random number
> >         generators of both
> >         > > > Numarray and Numeric.
> >         > > > I can't replicate it when calling the=20
> function from the
> >         python cosonle,=20
> >         > but
> >         > > > it has consistently malfunctioned when used within one
> >         of my scripts.
> >         > > >
> >         > > > What happens is that it frequently return a value
> >         greater than zero when
> >         > > > called with zero mean: poisson( 0.0)
> >         > > >
> >         > > > Unfortunately My program is too big to send=20
> attached but
> >         I have
> >         > confirmed
> >         > > > the malfunction by printing both the mean and=20
> the result
> >         whenever it
> >         > spits=20
> >         > > > out a wrong result.
> >         > > >
> >         > > > This is very weird indeed, I have run poisson millions
> >         of times by itsel
> >         > on
> >         > > > the python console, without any problems...
> >         > > >
> >         > > > I hope it is some stupid mistake, but when I=20
> replace the
> >         poisson
> >         > function
> >         > > > call within my program by the R equivalent command
> >         (rpois) via the rpy
> >         > > > wrapper, everything works just fine...=20
> >         > > >
> >         > > > any Ideas?
> >         > >
> >         >
> >         >
> >         >
> >         > --
> >         >
> >         > Fl=C3=A1vio Code=C3=A7o Coelho
> >         > registered Linux user # 386432
> >         > ---------------------------
> >         > Great spirits have always found violent opposition from
> >         mediocrities. The=20
> >         > latter cannot understand it when a man does not
> >         thoughtlessly submit to
> >         > hereditary prejudices but honestly and=20
> courageously uses his
> >         intelligence.
> >         >         Albert Einstein
> >         >
> >=20
> >=20
> >=20
> > --=20
> > Fl=C3=A1vio Code=C3=A7o Coelho
> > registered Linux user # 386432
> > ---------------------------
> > Great spirits have always found violent opposition from=20
> mediocrities.
> > The
> > latter cannot understand it when a man does not thoughtlessly submit
> > to=20
> > hereditary prejudices but honestly and courageously uses his
> > intelligence.
> >         Albert Einstein=20
>=20
>=20
>=20
> -------------------------------------------------------
> SF.Net email is sponsored by: Discover Easy Linux Migration Strategies
> from IBM. Find simple to follow Roadmaps, straightforward articles,
> informative Webcasts and more! Get everything you need to get up to
> speed, fast. http://ads.osdn.com/?ad_id=3D7477&alloc_id=3D16492&op=3Dclick
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>=20


=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D
This message is confidential. It may also be privileged or otherwise protec=
ted by work product immunity or other legal rules. If you have received it =
by mistake please let us know by reply and then delete it from your system;=
 you should not copy it or disclose its contents to anyone. All messages se=
nt to and from Electrabel may be monitored to ensure compliance with intern=
al policies and to protect our business. Emails are not secure and cannot b=
e guaranteed to be error free as they can be intercepted, amended, lost or =
destroyed, or contain viruses. Anyone who communicates with us by email is =
taken to accept these risks.

http://www.electrabel.be/homepage/general/disclaimer_EN.asp
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D

[Numpy-discussion] ANN: PyTables 1.1 released

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

=========================
 Announcing PyTables 1.1
=========================

The PyTables development team is happy to announce the availability of
a new version of PyTables package.

On this version you will find support for a nice set of new features,
like nested datatypes, enumerated datatypes, nested iterators (for
reading only), support for native HDF5 multidimensional attributes, a
new object for dealing with compressed, non-enlargeable arrays
(CArray), bzip2 compression support and more. Many bugs has been
addressed as well.

Go to the PyTables web site for downloading the beast:
http://pytables.sourceforge.net/

Or keep reading for more info about the new features and bugs fixed.


Changes more in depth
=====================

Improvements:

- Support for nested datatypes is in place. You can now made columns
  of tables that host another columns for an unlimited depth (well,
  theoretically, in practice until the python recursive limit would be
  reached). Convenient NestedRecArray objects has been implemented as
  data containers. Cols and Description accessors has been improved so
  you can navigate on the type hierarchy very easily (natural naming
  is has been implemented for the task).

- ``Table``, ``EArray`` and ``VLArray`` objects now support enumerated
  types.  ``Array`` objects support opening existing HDF5 enumerated
  arrays.  Enumerated types are restricted sets of ``(name, value)``
  pairs.  Use the ``Enum`` class to easily define new enumerations
  that will be saved along with your data.

- Now, the HDF5 library is responsible to do data conversions when the
  datasets are written in a machine with different byte-ordering than
  the machine that reads the dataset. With this, all the data is
  converted on-the-fly and you always get native datatypes in
  memory. I think this approach to be more convenient in terms of CPU
  consumption when using these datasets. Right now, this only works
  for tables, though.

- Added support for native HDF5 multidimensional attributes. Now, you
  can load native HDF5 files that contains fully multidimensional
  attributes; these attributes will be mapped to NumArray
  objects. Also, when you save NumArray objects as attributes, they
  get saved as native HDF5 attributes (before, NumArray attributes
  where pickled).

- A brand-new class, called CArray, has been introduced. It's mainly
  like an Array class (i.e. non-enlargeable), but with compression
  capabilities enabled. The existence of CArray also allows PyTables
  to read native HDF5 chunked, non-enlargeable datasets.

- Bzip2 compressor is supported. Such a support was already in
  PyTables 1.0, but forgot to announce it.

- New LZO2 (http://www.oberhumer.com/opensource/lzo/lzonews.php)
  compressor is supported. The installer now recognizes whether LZO1
  or LZO2 is installed, and adapts automatically to it. If both are
  installed in your system, then LZO2 is chosen. LZO2 claims to be
  fully compatible (both backward and forward) with LZO1, so you
  should not experience any problem during this transition.

- The old limit of 256 columns in a table has been released. Now, you
  can have tables with any number of columns, although if you try to
  use a too high number (i.e. > 1024), you will start to consume a lot
  of system resources. You have been warned!.

- The limit in the length of column names has been released also.

- Nested iterators for reading in tables are supported now.

- A new section in tutorial about how to modify values in tables and
  arrays has been added to the User's Manual.

Backward-incompatible changes:

- None.

Bug fixes:

- VLArray now correctly updates the number of rows internal counter
  when opening an existing VLArray object. Now you can add new rows
  to existing VLA's without problems.

- Tuple flavor for VLArrays now works as intended, i.e. reading
  VLArray objects will always return tuples even in the case of
  multidimensional Atoms. Before, this operations returned a mix of
  tuples and lists.

- If a column was not able to be indexed because it has too few
  entries, then _whereInRange is called instead of
  _whereIndexed. Fixes #1203202.

- You can call now Row.append() in the middle of Table iterators without
  resetting loop counters. Fixes #1205588.

- PyTables used to give a segmentation fault when removing the last
  row out of a table with the table.removeRows() method. This is due
  to a limitation in the HDF5 library. Until this get fixed in HDF5, a
  NotImplemented error is raised when trying to do that. Address
  #1201023.

- You can safely break a loop over an iterator returned by
  Table.where(). Fixes #1234637.

- When removing a Group with hidden child groups, those are
  effectively closed now.

- Now, there is a distinction between shapes 1 and (1,) in tables. The
  former represents a scalar, and the later a 1-D array with just one
  element. That follows the numarray convention for records, and makes
  more sense as well. Before 1.1, shapes 1 and (1,) were
  represented by an scalar on disk.

Known bugs:

- Classes inheriting from IsDescription subclasses do not
  inherit columns defined in the super-class. See SF bug #1207732 for
  more info.

- Time datatypes are non-portable between big-endian and little-endian
  architectures. This is ultimately a consequence of a HDF5
  limitation. See SF bug #1234709 for more info.

Known issues:

- UCL compressor seems to work badly on MacOSX platforms. Until the
  problem would be isolated and eventually solved, UCL will not be
  compiled by default on MacOSX platforms, even if the installer finds
  it in the system. However, if you still want to get UCL support on
  MacOSX, you can use the --force-ucl flag in setup.py.


Important note for Python 2.4 and Windows users
===============================================

If you are willing to use PyTables with Python 2.4 in Windows
platforms, you will need to get the HDF5 library compiled for MSVC
7.1, aka .NET 2003.  It can be found at:
ftp://ftp.ncsa.uiuc.edu/HDF/HDF5/current/bin/windows/5-164-win-net.ZIP

Users of Python 2.3 on Windows will have to download the version of
HDF5 compiled with MSVC 6.0 available in:
ftp://ftp.ncsa.uiuc.edu/HDF/HDF5/current/bin/windows/5-164-win.ZIP


What it is
==========

**PyTables** is a package for managing hierarchical datasets and
designed to efficiently cope with extremely large amounts of data
(with support for full 64-bit file addressing).  It features an
object-oriented interface that, combined with C extensions for the
performance-critical parts of the code, makes it a very easy-to-use
tool for high performance data storage and retrieval.

Perhaps its more interesting feature is that it optimizes memory and
disk resources so that data take much less space (between a factor 3
to 5, and more if the data is compressible) than other solutions, like
for example, relational or object oriented databases.

Besides, PyTables I/O for table objects is buffered, implemented in C
and carefully tuned so that you can reach much better performance with
PyTables than with your own home-grown wrappings to the HDF5
library. PyTables sports indexing capabilities as well, allowing doing
selections in tables exceeding one billion of rows in just seconds.


Where can PyTables be applied?
==============================

PyTables is not designed to work as a relational database competitor,
but rather as a teammate.  If you want to work with large datasets of
multidimensional data (for example, for multidimensional analysis), or
just provide a categorized structure for some portions of your
cluttered RDBS, then give PyTables a try.  It works well for storing
data from data acquisition systems (DAS), simulation software, network
data monitoring systems (for example, traffic measurements of IP
packets on routers), very large XML files, or for creating a
centralized repository for system logs, to name only a few possible
uses.


What is a table?
================

A table is defined as a collection of records whose values are stored
in fixed-length fields.  All records have the same structure and all
values in each field have the same data type.  The terms
"fixed-length" and "strict data types" seem to be quite a strange
requirement for a language like Python that supports dynamic data
types, but they serve a useful function if the goal is to save very
large quantities of data (such as is generated by many scientific
applications, for example) in an efficient manner that reduces demand
on CPU time and I/O resources.


What is HDF5?
=============

For those people who know nothing about HDF5, it is a general purpose
library and file format for storing scientific data made at NCSA.
HDF5 can store two primary objects: datasets and groups.  A dataset is
essentially a multidimensional array of data elements, and a group is
a structure for organizing objects in an HDF5 file.  Using these two
basic constructs, one can create and store almost any kind of
scientific data structure, such as images, arrays of vectors, and
structured and unstructured grids.  You can also mix and match them in
HDF5 files according to your needs.


Platforms
=========

We are using Linux on top of Intel32 as the main development platform,
but PyTables should be easy to compile/install on other UNIX machines.
This package has also been successfully compiled and tested on a
FreeBSD 5.4 with Opteron64 processors, a UltraSparc platform with
Solaris 7 and Solaris 8, a SGI Origin3000 with Itanium processors
running IRIX 6.5 (using the gcc compiler), Microsoft Windows and
MacOSX (10.2 although 10.3 should work fine as well). In particular,
it has been thoroughly tested on 64-bit platforms, like Linux-64 on
top of an Intel Itanium, AMD Opteron (in 64-bit mode) or PowerPC G5
(in 64-bit mode) where all the tests pass successfully.

Regarding Windows platforms, PyTables has been tested with Windows
2000 and Windows XP (using the Microsoft Visual C compiler), but it
should also work with other flavors as well.


Web site
========

Go to the PyTables web site for more details:

http://pytables.sourceforge.net/

To know more about the company behind the PyTables development, see:

http://www.carabos.com/


Share your experience
=====================

Let us know of any bugs, suggestions, gripes, kudos, etc. you may
have.


----

  **Enjoy data!**

  -- The PyTables Team

Re: [Numpy-discussion] problem with poisson generators

[Todd M. <jm...@st...>]

On Wed, 2005-07-13 at 12:13 -0300, Flavio Coelho wrote:
> 
> 
> 2005/7/13, Bruce Southey <bso...@gm...>:
>         Hi,
>         What is the seed used?
> 
> I am not setting the seed.
>  
> 
>         It is really important that you set the seed?
> 
> No.
> 
> 
>         Did you build Python and numarray from source?
> 
> 
> Yes,  I  use Gentoo. I  build everything from source. 
> 
> 
>         Can you reduce your code to a few lines that have the problem?
> 
> 
> Not really, poisson and binomial are the only two functions from
> Numeric that I use but they are called inside a rather complex object
> oriented code environment (objects within objetcs, being called
> recursively...) Bellow is the function within which poisson is called:
> 
> def stepSEIR_s(self,inits=(0,0,0),par=(0.001,1,1,0.5,1,0,0,0),theta=0,
> npass=0,dist='poisson'):
>         """
>         Defines an stochastic model SEIR:
>         - inits = (E,I,S)
>         - par = (Beta, alpha, E,r,delta,B,w,p) see docs.
>         - theta = infectious individuals from neighbor sites
>         """
>         E,I,S = inits
>         N = self.parentSite.totpop
>         beta,alpha,e,r,delta,B,w,p = par
>         Lpos_esp = beta*S*((I+theta)/(N+npass))**alpha #Number of new
> cases
>         
>         if dist == 'poisson':
>             Lpos = poisson(Lpos_esp)
> ##            if theta == 0 and Lpos_esp == 0 and Lpos > 0:  
> ##                print
> Lpos,Lpos_esp,S,I,theta,N,self.parentSite.sitename
>         elif dist == 'negbin':
>             prob = I/(I+Lpos_esp) #convertin between parameterizations
>             Lpos = negative_binomial(I,prob)
>         self.parentSite.totalcases += Lpos #update number of cases    
>         Epos = (1-e)*E + Lpos
>         Ipos = e*E + (1-r)*I
>         Spos = S + B - Lpos 
>         Rpos = N-(Spos+Epos+Ipos)
>         #self.parentSite.totpop = Spos+Epos+Ipos+Rpos
>         self.parentSite.incidence.append(Lpos)
>         if not self.parentSite.infected:
>             if Lpos > 0:
>                 self.parentSite.infected =
> self.parentSite.parentGraph.simstep
>                 self.parentSite.parentGraph.epipath.append
> ((self.parentSite.parentGraph.simstep,self.parentSite,self.parentSite.infector))
> 
>         
>         return [Epos,Ipos,Spos]
> 

Having this code is a good start to solving the problem.  I think the
next step is to simplify your example to make it runnable and provide
known inputs for all the parameters which lead to a failure for you. 

Being really literal (spelling out every single damn thing) cuts down on
speculation and inefficiency on our end.

It would also be good to express the condition (or it's inverse) you
think is an error as an assertion,  so something like this might be what
we need:

from numarray.random_array import poisson

E, I, S = (0,0,0)
beta,alpha,e,r,delta,B,w,p = (0.001,1,1,0.5,1,0,0,0)
theta = 0
npass = 0
N = 100 # total guess here for me
Lpos_esp = beta*S*((I+theta)/(N+npass))**alpha #Number of new cases
Lpos = poisson(Lpos_esp)
assert not (theta == 0 and Lpos_esp == 0 and Lpos > 0)

The problem is,  the above works for me.  Make it sane and get it to
expose the error for you and I'll look into this some more.

Regards,
Todd

> 
> I wonder if called by itself  it would trigger the problem... The
> commented Lines Is where I catch the errors: when poisson returns a
> greater than zero number, when called with mean 0.
> 
> 
> 
>         Ranlib uses static floats so it may relate to numerical
>         precision (as
>         well as the random uniform random number generator). Really
>         the only 
>         way is for you to debug the ranlib.c file
> 
> I'll see what I can do. 
> 
> 
>         Note that R provides a standalone math library in C that
>         includes the
>         random number generators so you could you those instead. SWIG
>         handles 
>         it rather well.
> 
> 
> I think thats what Rpy already does...is it not? 
> 
> thanks Bruce,
> 
> Flávio
>  
> 
>         Regards
>         Bruce
>         
>         On 7/13/05, Flavio Coelho <fcc...@gm...> wrote:
>         > Well,
>         >   I am definetly glad that someone has also stumbled onto
>         the same problem.
>         >
>         >  But it is nevertheless intriguing, that you can run poisson
>         a million times 
>         > with mean  zero or negative(it assumes zero mean inthis
>         case) without any
>         > problems by itself. But when I call it within my code, the
>         rate of error is
>         > very high (I would say it returns a wrong result every time,
>         but I haven't 
>         > checked every answer).
>         >
>         >  Meanwhile, my solution will be:
>         >
>         >  import rpy
>         >
>         >  n = rpy.r.rpois(n,mean)
>         >
>         >  I don't feel I can trust poisson  while this "funny"
>         behavior is still 
>         > there...
>         >  If someone has any Idea of how I could trace this bug
>         please tell me, and
>         > I'll hunt it down. At least I can reproduce it in a very
>         specific context.
>         >
>         >  thanks,
>         >
>         >  Flávio 
>         >
>         > 2005/7/12, Sebastian Haase <ha...@ms...>:
>         > > Hi Flavio!
>         > > I had reported this long time ago and this list (about
>         numarray).
>         > > Somehow this got more or less dismissed. If I recall
>         correctly the 
>         > argument
>         > > was that nobody could reproduce it (I ran this on Debian
>         Woody , py2.2,
>         > (with
>         > > CVS numarray at the time).
>         > >
>         > > I ended up writting my own wrapper(s): 
>         > > def poissonArr(shape=defshape, mean=1):
>         > >         from numarray import random_array as ra
>         > >         if mean == 0:
>         > >                 return zeroArrF(shape)
>         > >         elif mean < 0: 
>         > >                 raise "poisson not defined for mean < 0"
>         > >         else:
>         > >                 return ra.poisson(mean, shape).astype
>         (na.UInt16)
>         > >
>         > > def poissonize(arr):
>         > >         from numarray import random_array as ra
>         > >         return na.where(arr<=0, 0, ra.poisson(arr)).astype
>         ( na.UInt16)
>         > >
>         > > (I use the astype( na.UInt16) because of some OpenGL code)
>         > >
>         > > Just last week had this problem on a windows98 computer
>         (python2.4).
>         > >
>         > > This should get sorted out ... 
>         > >
>         > > Thanks for reporting this problem.
>         > > Sebastian Haase
>         > >
>         > >
>         > >
>         > >
>         > > On Tuesday 12 July 2005 13:32, Flavio Coelho wrote:
>         > > > Hi, 
>         > > >
>         > > > I am having problems with the poisson random number
>         generators of both
>         > > > Numarray and Numeric.
>         > > > I can't replicate it when calling the function from the
>         python cosonle, 
>         > but
>         > > > it has consistently malfunctioned when used within one
>         of my scripts.
>         > > >
>         > > > What happens is that it frequently return a value
>         greater than zero when
>         > > > called with zero mean: poisson( 0.0)
>         > > >
>         > > > Unfortunately My program is too big to send attached but
>         I have
>         > confirmed
>         > > > the malfunction by printing both the mean and the result
>         whenever it
>         > spits 
>         > > > out a wrong result.
>         > > >
>         > > > This is very weird indeed, I have run poisson millions
>         of times by itsel
>         > on
>         > > > the python console, without any problems...
>         > > >
>         > > > I hope it is some stupid mistake, but when I replace the
>         poisson
>         > function
>         > > > call within my program by the R equivalent command
>         (rpois) via the rpy
>         > > > wrapper, everything works just fine... 
>         > > >
>         > > > any Ideas?
>         > >
>         >
>         >
>         >
>         > --
>         >
>         > Flávio Codeço Coelho
>         > registered Linux user # 386432
>         > ---------------------------
>         > Great spirits have always found violent opposition from
>         mediocrities. The 
>         > latter cannot understand it when a man does not
>         thoughtlessly submit to
>         > hereditary prejudices but honestly and courageously uses his
>         intelligence.
>         >         Albert Einstein
>         >
> 
> 
> 
> -- 
> Flávio Codeço Coelho
> registered Linux user # 386432
> ---------------------------
> Great spirits have always found violent opposition from mediocrities.
> The
> latter cannot understand it when a man does not thoughtlessly submit
> to 
> hereditary prejudices but honestly and courageously uses his
> intelligence.
>         Albert Einstein 

Re: [Numpy-discussion] problem with poisson generators

[Flavio C. <fcc...@gm...>]

2005/7/14, Robert Kern <rk...@uc...>:
>=20
> Flavio Coelho wrote:
>=20
> > I have both numeric 23.7 and numarray 1.3.1 installed and on neither
> > of them I could reproduce the bug when I called them directly from the
> > python interpreter. However they fail on mean 0.0 every time when calle=
d
> > within my code. So it appears to me that the Bug you mentioned is not
> > what is causing my problem. It seems to stem from interaction with the
> > way its being called, maybe some carryover from previous calls...
> > this is the test I ran on the interpreter:
> > [(poisson(i),i) for i in uniform(-20,20,1000) if i<=3D0]
> >
> > I also ran:
> >
> > sum(poisson(0,100000))
> >
> > they both worked flawlessly. In the first test I wanted to see if there
> > was some carry over from previous runs when called with various means
> > (which is closer to the way it is used within my code), but it returned
> > zero every time. (I don't use negative means in my code, but I wanted t=
o
> > try it here just in case)
>=20
> Are you sure that you don't have other versions of Numeric or numarray
> installed that might be getting picked up by your code?



I have just double checked. And the answer is no, there is a single version=
=20
of each of the modules.=20

--
> Robert Kern
> rk...@uc...
>=20
> "In the fields of hell where the grass grows high
> Are the graves of dreams allowed to die."
> -- Richard Harter
>=20
>=20
> -------------------------------------------------------
> This SF.Net <http://SF.Net> email is sponsored by the 'Do More With Dual!=
'=20
> webinar happening
> July 14 at 8am PDT/11am EDT. We invite you to explore the latest in dual
> core and dual graphics technology at this free one hour event hosted by=
=20
> HP,
> AMD, and NVIDIA. To register visit http://www.hp.com/go/dualwebinar
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>=20



--=20
Fl=E1vio Code=E7o Coelho
registered Linux user # 386432
---------------------------
Great spirits have always found violent opposition from mediocrities. The
latter cannot understand it when a man does not thoughtlessly submit to
hereditary prejudices but honestly and courageously uses his intelligence.
Albert Einstein

Re: [Numpy-discussion] numarray 1.3.3 build question

[Todd M. <jm...@st...>]

On Wed, 2005-07-13 at 10:09 -0700, Russell E. Owen wrote:
> I'm on MacOS X 10.3.9 and am getting mysterious warnings when building 
> numarray 1.3.3 which I've never seen before.
> 
> Building for the unix/x11 version yields:
> [dhcp-254:~/Archives\342\200\242/PythonPackages/numarray-1.3.3] rowen% 
> python setup.py config build --gencode
> Running sitecustomize.pyc
> Using EXTRA_COMPILE_ARGS = ['-Ddarwin']
> Running sitecustomize.pyc
> Running sitecustomize.pyc
> generating new API module 'libnumarray' .c  & .h
> generating new API module 'libteacup' .c  & .h
> generating new API module 'libnumeric' .c  & .h
> Using external BLAS and LAPACK
> running config
> Wrote config.h
> ...
> Packages/LinearAlgebra2/Src/lapack_litemodule.c:851: warning: function 
> declaration isn't a prototype
> Packages/LinearAlgebra2/Src/lapack_litemodule.c:848: warning: 
> `lapack_liteError' defined but not used
> gcc: -framework: linker input file unused because linking not done
> gcc: vecLib: linker input file unused because linking not done
> 
> Should I worry? 



> I haven't figured out how to run the self-tests w/out 
> installing and I'm reluctant to install w/out knowing if the warnings 
> are a problem.
> 
> also:
> - why do we need --gencode? 

You don't need it,  I do.  In general,  it's not necessary to use 
but it's there because the code generation needs to remain "live" in
order to be an asset but is also time consuming for people doing
repetitive installs to fix problems.   Prior to Python-2.4 and VC.NET,
--gen_code had meaning on win32 because VC6.0 didn't fully support
UInt64.  Everywhere else,  the CVS'ed C source just works.

> Is there some way to make that the default 
> so we can leave it off?

Sure.  I would say just don't worry about it unless you're building on
win32 from source using VC6.0.  

> - why do we need "config"? That's new to 1.3.3 isn't it?

No,  I think that's been there since 1.2.1 or so.  It's an artifact of
user submitted fixes to the linear algebra package which I accepted as-
is because they were submitted by someone who both knew enough
technically and had enough idealism to make the effort.

Regards,
Todd

Re: [Numpy-discussion] bug or feature?

[Todd M. <jm...@st...>]

On Mon, 2005-07-11 at 15:22 -0400, Perry Greenfield wrote:
> For array indexing, assignment does 
> modify the original array (as your case above illustrates), but once 
> you sliced the array-indexed array, it changed the context of the 
> indexing. I.e., (to take a simpler, 1-d array x)
> 
> x[[0,2]] = 1  # modifies x
> 
> This form results in __setitem__ being called on x
> 
> x[[0,2]][:4] = 1  # doesn't modify x
> 
> This form results in __getitem__ being called on x (and thus producing 
> a new copy of the array) for the array indexing and __setitem__ being 
> called for the slice. At least that's what I think is happening. Todd 
> or someone more recently familiar with how Python handles this can 
> correct me if I'm wrong.
> 
> Perry
> 

This is indeed what happens.  A temporary is created by the first sub-
expression x[[0,2]],   modified by the slice subexpression [:4] = 1, and
then discarded.  I think the array indexing "grammar" could be changed
to support x[[0,2],:4] as a single expression which would create the
temporary and then copy-back at the end...  but numarray's array
indexing is not currently that fancy.

Todd

Re: [Numpy-discussion] problem with poisson generators

[Robert K. <rk...@uc...>]

Flavio Coelho wrote:

> I  have  both numeric 23.7 and numarray 1.3.1 installed  and on neither 
> of them I could reproduce the bug when I called them directly from the 
> python interpreter. However they fail on mean 0.0 every time when called 
> within my code. So it appears to me that the Bug you mentioned is not 
> what is causing my problem. It seems to stem from interaction with the 
> way its being called, maybe some carryover from previous calls...
> this is the test I ran on the interpreter:
> [(poisson(i),i) for i in uniform(-20,20,1000) if i<=0]
> 
> I also ran:
> 
>  sum(poisson(0,100000))
> 
> they both worked flawlessly. In the first test I wanted to see if there 
> was some carry over from previous runs when called with various means 
> (which is closer to the way it is used within my code), but it returned 
> zero every time. (I don't use negative means in my code, but I wanted to 
> try it here just in case)

Are you sure that you don't have other versions of Numeric or numarray 
installed that might be getting picked up by your code?

-- 
Robert Kern
rk...@uc...

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

Re: [Numpy-discussion] problem with poisson generators

[Flavio C. <fcc...@gm...>]

2005/7/13, David M. Cooke <co...@ph...>:
>=20
> On Tue, Jul 12, 2005 at 05:32:25PM -0300, Flavio Coelho wrote:
> > Hi,
> >
> > I am having problems with the poisson random number generators of both
> > Numarray and Numeric.
> > I can't replicate it when calling the function from the python cosonle,=
=20
> but
> > it has consistently malfunctioned when used within one of my scripts.
> >
> > What happens is that it frequently return a value greater than zero whe=
n
> > called with zero mean: poisson(0.0)
> >
> > Unfortunately My program is too big to send attached but I have=20
> confirmed
> > the malfunction by printing both the mean and the result whenever it=20
> spits
> > out a wrong result.
> >
> > This is very weird indeed, I have run poisson millions of times by itse=
l=20
> on
> > the python console, without any problems...
> >
> > I hope it is some stupid mistake, but when I replace the poisson=20
> function
> > call within my program by the R equivalent command (rpois) via the rpy
> > wrapper, everything works just fine...
> >
> > any Ideas?
>=20
> This looks like bug #1123145 in Numeric:
>=20
>=20
> http://sourceforge.net/tracker/index.php?func=3Ddetail&aid=3D1123145&grou=
p_id=3D1369&atid=3D101369
>=20
> which was fixed a few months ago. numarray, I believe, originally took
> ranlib.c from Numeric, so it doesn't have this bug fix. Try replacing
> numarray's ranlib.c with the version from Numeric 24.0b2 (or CVS).



I have both numeric 23.7 and numarray 1.3.1 installed and on neither of the=
m=20
I could reproduce the bug when I called them directly from the python=20
interpreter. However they fail on mean 0.0 every time when called within my=
=20
code. So it appears to me that the Bug you mentioned is not what is causing=
=20
my problem. It seems to stem from interaction with the way its being called=
,=20
maybe some carryover from previous calls...
this is the test I ran on the interpreter:
[(poisson(i),i) for i in uniform(-20,20,1000) if i<=3D0]=20

I also ran:

sum(poisson(0,100000))

they both worked flawlessly. In the first test I wanted to see if there was=
=20
some carry over from previous runs when called with various means (which is=
=20
closer to the way it is used within my code), but it returned zero every=20
time. (I don't use negative means in my code, but I wanted to try it here=
=20
just in case)

--
> |>|\/|<
>=20
> /------------------------------------------------------------------------=
--\
> |David M. Cooke http://arbutus.physics.mcmaster.ca/dmc/
> |co...@ph... <http://physics.mcmaster.ca>
>=20



--=20
Fl=E1vio Code=E7o Coelho
registered Linux user # 386432
---------------------------
Great spirits have always found violent opposition from mediocrities. The
latter cannot understand it when a man does not thoughtlessly submit to
hereditary prejudices but honestly and courageously uses his intelligence.
Albert Einstein