A random.normal function with stdev as array
Hi, I am trying to optimize a code where I derive random numbers many times and having an array of values for the stdev parameter. I wish to have an efficient way of doing something like: ################## stdev = array([1.1,1.2,1.0,2.2]) result = numpy.zeros(stdev.shape, Float) for i in range(len(stdev)) : result[i] = numpy.random.normal(0, stdev[i]) ################## In my case, stdev can in fact be an array of a few millions floats... so I really need to optimize things. Any hint on how to code this efficiently ? And in general, where could I find tips for optimizing a code where I unfortunately have too many loops such as "for i in range(Nbody) : " with Nbody being > 10^6 ? thanks! Eric
Eric Emsellem wrote:
Hi,
I am trying to optimize a code where I derive random numbers many times and having an array of values for the stdev parameter.
I wish to have an efficient way of doing something like: ################## stdev = array([1.1,1.2,1.0,2.2]) result = numpy.zeros(stdev.shape, Float) for i in range(len(stdev)) : result[i] = numpy.random.normal(0, stdev[i]) ##################
You can use the fact that the standard deviation of a normal distribution is a scale parameter. You can get random normal deviates of varying standard deviation by multiplying a standard normal deviate by the desired standard deviation (how's that for confusing terminology, eh?). result = numpy.random.standard_normal(stdev.shape) * stdev
In my case, stdev can in fact be an array of a few millions floats... so I really need to optimize things.
Any hint on how to code this efficiently ?
And in general, where could I find tips for optimizing a code where I unfortunately have too many loops such as "for i in range(Nbody) : " with Nbody being > 10^6 ?
Tim Hochberg recently made this list: """ 0. Think about your algorithm. 1. Vectorize your inner loop. 2. Eliminate temporaries 3. Ask for help 4. Recode in C. 5. Accept that your code will never be fast. Step zero should probably be repeated after every other step ;) """ That's probably the best general advice. To get better advice, we would need to know the specifics of the problem. -- Robert Kern robert.kern@gmail.com "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
Hi Eric, In the past , I've done things like ###### normdist = lambda x: numpy.random.normal(0,x) vecnormal = numpy.vectorize(normdist) stdev = numpy.array([1.1,1.2,1.0,2.2]) result = vecnormal(stdev) ###### This works fine for up to 10k elements for stdev for some reason. Any larger then that and i get a Bus error on my PPC mac and a segfault on my x86 linux box. I'm running numpy 0.9.7.2325 on both machines. Perhaps for larger inputs, you could break up your loop into smaller vectorized chunks. Regards, John On Wed, Apr 05, 2006 at 03:32:06PM +0200, Eric Emsellem wrote:
Hi,
I am trying to optimize a code where I derive random numbers many times and having an array of values for the stdev parameter.
I wish to have an efficient way of doing something like: ################## stdev = array([1.1,1.2,1.0,2.2]) result = numpy.zeros(stdev.shape, Float) for i in range(len(stdev)) : result[i] = numpy.random.normal(0, stdev[i]) ##################
In my case, stdev can in fact be an array of a few millions floats... so I really need to optimize things.
Any hint on how to code this efficiently ?
And in general, where could I find tips for optimizing a code where I unfortunately have too many loops such as "for i in range(Nbody) : " with Nbody being > 10^6 ?
thanks! Eric
------------------------------------------------------- This SF.Net email is sponsored by xPML, a groundbreaking scripting language that extends applications into web and mobile media. Attend the live webcast and join the prime developer group breaking into this new coding territory! http://sel.as-us.falkag.net/sel?cmd=lnk&kid=110944&bid=241720&dat=121642 _______________________________________________ Numpy-discussion mailing list Numpy-discussion@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/numpy-discussion
-- If liberty and equality, as is thought by some are chiefly to be found in democracy, they will be best attained when all persons alike share in the government to the utmost. -- Aristotle, Politics
John Byrnes wrote:
Hi Eric,
In the past , I've done things like
###### normdist = lambda x: numpy.random.normal(0,x) vecnormal = numpy.vectorize(normdist)
stdev = numpy.array([1.1,1.2,1.0,2.2]) result = vecnormal(stdev)
######
This works fine for up to 10k elements for stdev for some reason. Any larger then that and i get a Bus error on my PPC mac and a segfault on my x86 linux box.
This needs to be tracked down. It looks like some-kind of error is not being caught correctly. You should not get a segfault. Could you provide a stack-trace when the problem occurs? One issue is that vectorize is using object arrays under the covers which is consuming roughly 2x the memory than you may think. An object array is created and the function is called for every element. This object array is then converted to a number type after the fact. The segfault should be tracked down in any case. -Travis
On Thu, Apr 06, 2006 at 03:36:59AM -0600, Travis Oliphant wrote:
John Byrnes wrote:
Hi Eric,
In the past , I've done things like
###### normdist = lambda x: numpy.random.normal(0,x) vecnormal = numpy.vectorize(normdist)
stdev = numpy.array([1.1,1.2,1.0,2.2]) result = vecnormal(stdev)
######
This works fine for up to 10k elements for stdev for some reason. Any larger then that and i get a Bus error on my PPC mac and a segfault on my x86 linux box.
This needs to be tracked down. It looks like some-kind of error is not being caught correctly. You should not get a segfault. Could you provide a stack-trace when the problem occurs?
One issue is that vectorize is using object arrays under the covers which is consuming roughly 2x the memory than you may think. An object array is created and the function is called for every element. This object array is then converted to a number type after the fact.
The segfault should be tracked down in any case.
-Travis
Hi Travis, Here is a backtrace from gdb on my mac. John #0 0x00470b88 in log1pl () #1 0x00000000 in ?? () Cannot access memory at address 0x0 Cannot access memory at address 0x0 #2 0x004708ec in log1pl () #3 0x1000c348 in PyObject_Call (func=0x4, arg=0x4, kw=0x15fb) at /Users/bob/src/Python-2.4.1/Objects/abstract.c:1751 #4 0x1007ce34 in ext_do_call (func=0x1, pp_stack=0xbfffed90, flags=211904, na=8656012, nk=1194304) at /Users/bob/src/Python-2.4.1/Python/ceval.c:3824 #5 0x1007a230 in PyEval_EvalFrame (f=0x848410) at /Users/bob/src/Python-2.4.1/Python/ceval.c:2203 #6 0x1007b284 in PyEval_EvalCodeEx (co=0x2, globals=0x4, locals=0x1, args=0x3, argcount=1049072, kws=0x841150, kwcount=1, defs=0x8411fc, defcount=0, closure=0x0) at /Users/bob/src/Python-2.4.1/Python/ceval.c:2730 #7 0x10026274 in function_call (func=0x880bb0, arg=0x1001f0, kw=0x848410) at /Users/bob/src/Python-2.4.1/Objects/funcobject.c:548 #8 0x1000c348 in PyObject_Call (func=0x4, arg=0x4, kw=0x15fb) at /Users/bob/src/Python-2.4.1/Objects/abstract.c:1751 #9 0x10015a88 in instancemethod_call (func=0x52eef0, arg=0x54a170, kw=0x0) at /Users/bob/src/Python-2.4.1/Objects/classobject.c:2431 #10 0x1000c348 in PyObject_Call (func=0x4, arg=0x4, kw=0x15fb) at /Users/bob/src/Python-2.4.1/Objects/abstract.c:1751 #11 0x10059358 in slot_tp_call (self=0x53e4f0, args=0x5b310, kwds=0x0) at /Users/bob/src/Python-2.4.1/Objects/typeobject.c:4526 #12 0x1000c348 in PyObject_Call (func=0x4, arg=0x4, kw=0x15fb) at /Users/bob/src/Python-2.4.1/Objects/abstract.c:1751 #13 0x1007c9e4 in do_call (func=0x53e4f0, pp_stack=0x53e4f0, na=0, nk=8655844) at /Users/bob/src/Python-2.4.1/Python/ceval.c:3755 #14 0x1007c6dc in call_function (pp_stack=0x0, oparg=4) at /Users/bob/src/Python-2.4.1/Python/ceval.c:3570 #15 0x1007a140 in PyEval_EvalFrame (f=0x10e200) at /Users/bob/src/Python-2.4.1/Python/ceval.c:2163 #16 0x1007c83c in fast_function (func=0x4, pp_stack=0x10e360, n=268927488, na=268755664, nk=1) at /Users/bob/src/Python-2.4.1/Python/ceval.c:3629 #17 0x1007c6c4 in call_function (pp_stack=0xbffff5bc, oparg=4) at /Users/bob/src/Python-2.4.1/Python/ceval.c:3568 #18 0x1007a140 in PyEval_EvalFrame (f=0x10e030) at /Users/bob/src/Python-2.4.1/Python/ceval.c:2163 #19 0x1007b284 in PyEval_EvalCodeEx (co=0x0, globals=0x4, locals=0x1, args=0x10078200, argcount=1049072, kws=0x841150, kwcount=1, defs=0x8411fc, defcount=0, closure=0x0) at /Users/bob/src/Python-2.4.1/Python/ceval.c:2730 #20 0x1007e678 in PyEval_EvalCode (co=0x4, globals=0x4, locals=0x15fb) at /Users/bob/src/Python-2.4.1/Python/ceval.c:484 #21 0x100b2ee0 in run_node (n=0x10078200, filename=0x4 <Address 0x4 out of bounds>, globals=0x0, locals=0x10e180, flags=0x2) at /Users/bob/src/Python-2.4.1/Python/pythonrun.c:1265 #22 0x100b23b0 in PyRun_InteractiveOneFlags (fp=0x54a1a5, filename=0x56ca0 "", flags=0x10e030) at /Users/bob/src/Python-2.4.1/Python/pythonrun.c:762 #23 0x100b2190 in PyRun_InteractiveLoopFlags (fp=0x56b94, filename=0xd440 "", flags=0x100f21b8) at /Users/bob/src/Python-2.4.1/Python/pythonrun.c:695 #24 0x100b3bb0 in PyRun_AnyFileExFlags (fp=0xa0001554, filename=0x100f36ac "<stdin>", closeit=0, flags=0xbffff934) at /Users/bob/src/Python-2.4.1/Python/pythonrun.c:658 #25 0x100bf640 in Py_Main (argc=269413412, argv=0x20000000) at /Users/bob/src/Python-2.4.1/Modules/main.c:484 #26 0x000018d0 in start () #27 0x8fe1a278 in __dyld__dyld_start ()
Hi, Can you provide more details on what you are doing, especially how you are using this? The one item that is not directly part of Tim's list is that some times you need to reorder your loops (perhaps this is part of "Think about your algorithm"?). Loop swapping is very common to improve performance. However, it usually requires a very clear head or someone else to do it. Also, you can might need to break loops into pieces where you repeat the same tasks and computations over and over. The other aspect is to do some algebra on the calculations as the stdev is essentially a constant so depending on how you use it you can factor it out further. Again it all depends on what you are actually doing with these numbers.
From a different view, you need to be very careful with your (pseudo)random number generator with that many samples. These have a tendency to repeat so your random number stream is no longer random. See the Wikipedia entry: http://en.wikipedia.org/wiki/Pseudorandom_number_generator
If I recall correctly, the Python random number generator is a Mersenne twister but ranlib is not and so prone to the mentioned problems. I do not know if SciPy adds any other generators. Finally I would also cheat by reducing the stdev values because in many cases you will not see a real difference between a normal with mean zero and variance 1.0 and a normal with mean zero and variance 1.1 (especially if you are doing more than comparing distributions so there are more sources of 'error') unless you have a really large number of samples. Regards Bruce On 4/5/06, Eric Emsellem <emsellem@obs.univ-lyon1.fr> wrote:
Hi,
I am trying to optimize a code where I derive random numbers many times and having an array of values for the stdev parameter.
I wish to have an efficient way of doing something like: ################## stdev = array([1.1,1.2,1.0,2.2]) result = numpy.zeros(stdev.shape, Float) for i in range(len(stdev)) : result[i] = numpy.random.normal(0, stdev[i]) ##################
In my case, stdev can in fact be an array of a few millions floats... so I really need to optimize things.
Any hint on how to code this efficiently ?
And in general, where could I find tips for optimizing a code where I unfortunately have too many loops such as "for i in range(Nbody) : " with Nbody being > 10^6 ?
thanks! Eric
------------------------------------------------------- This SF.Net email is sponsored by xPML, a groundbreaking scripting language that extends applications into web and mobile media. Attend the live webcast and join the prime developer group breaking into this new coding territory! http://sel.as-us.falkag.net/sel?cmd=lnk&kid=110944&bid=241720&dat=121642 _______________________________________________ Numpy-discussion mailing list Numpy-discussion@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/numpy-discussion
Bruce Southey wrote:
From a different view, you need to be very careful with your (pseudo)random number generator with that many samples. These have a tendency to repeat so your random number stream is no longer random. See the Wikipedia entry: http://en.wikipedia.org/wiki/Pseudorandom_number_generator
If I recall correctly, the Python random number generator is a Mersenne twister but ranlib is not and so prone to the mentioned problems. I do not know if SciPy adds any other generators.
numpy.random uses the Mersenne Twister. RANLIB is dead! Long live MT19937! -- Robert Kern robert.kern@gmail.com "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
participants (5)
-
Bruce Southey -
Eric Emsellem -
John Byrnes -
Robert Kern -
Travis Oliphant