Hi, What is the seed used? It is really important that you set the seed? Did you build Python and numarray from source? Can you reduce your code to a few lines that have the problem? 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 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. Regards Bruce On 7/13/05, Flavio Coelho 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 :

...

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

On Wed, 2005-07-13 at 12:13 -0300, Flavio Coelho wrote:

2005/7/13, Bruce Southey : 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 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 : > > 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

On Fri, 2005-07-15 at 21:30 -0300, Flavio Coelho wrote:

2005/7/15, Todd Miller :

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 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 : > > > 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

2005/7/13, David M. Cooke :

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, 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?

This looks like bug #1123145 in Numeric:

http://sourceforge.net/tracker/index.php?func=detail&aid=1123145&group_id=1369&atid=101369

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 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) --

|>|\/|<

/--------------------------------------------------------------------------\ |David M. Cooke http://arbutus.physics.mcmaster.ca/dmc/ |cookedm@physics.mcmaster.ca http://physics.mcmaster.ca

-- 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

2005/7/14, Robert Kern :

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?

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

Robert Kern rkern@ucsd.edu

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

------------------------------------------------------- This SF.Net http://SF.Net email is sponsored by the 'Do More With Dual!' 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 HP, AMD, and NVIDIA. To register visit http://www.hp.com/go/dualwebinar _______________________________________________ Numpy-discussion mailing list Numpy-discussion@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/numpy-discussion

-- 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