Mailman 3 January 2018 - NumPy-Discussion

Invalid value encoutered : how to prevent numpy.where to do this?
by Eric Emsellem 18 Feb '23

18 Feb '23

Dear all, I have a code using lots of "numpy.where" to make some constrained calculations as in: data = arange(10) result = np.where(data == 0, 0., 1./data) # or data1 = arange(10) data2 = arange(10)+1.0 result = np.where(data1 > data2, np.sqrt(data1-data2), np.sqrt(data2-data2)) which then produces warnings like: /usr/bin/ipython:1: RuntimeWarning: invalid value encountered in sqrt or for the first example: /usr/bin/ipython:1: RuntimeWarning: divide by zero encountered in divide How do I avoid these messages to appear? I know that I could in principle use numpy.seterr. However, I do NOT want to remove these warnings for other potential divide/multiply/sqrt etc errors. Only when I am using a "where", to in fact avoid such warnings! Note that the warnings only happen once, but since I am going to release that code, I would like to avoid the user to get such messages which are irrelevant here (because I am testing, with the where, when NOT to divide by zero or take a sqrt of a negative number). thanks! Eric

5 4

Re: [Numpy-discussion] example reading binary Fortran file
by Neil Martinsen-Burrell 23 Jul '21

23 Jul '21

David Froger <david.froger.info <at> gmail.com> writes: > Hy,My question is about reading Fortran binary file (oh no this question > again...) I've posted this before, but I finally got it cleaned up for the Cookbook. For this purpose I use a subclass of file that has methods for reading unformatted Fortran data. See http://www.scipy.org/Cookbook/FortranIO/FortranFile. I'd gladly see this in numpy or scipy somewhere, but I'm not sure where it belongs. > program makeArray > implicit none > integer,parameter:: nx=10,ny=20 > real(4),dimension(nx,ny):: ux,uy,p > integer :: i,j > open(11,file='uxuyp.bin',form='unformatted') > do i = 1,nx > do j = 1,ny > ux(i,j) = real(i*j) > uy(i,j) = real(i)/real(j) > p (i,j) = real(i) + real(j) > enddo > enddo > write(11) ux,uy > write(11) p > close(11) > end program makeArray When I run the above program compiled with gfortran on my Intel Mac, I can read it back with:: >>> import numpy as np >>> from fortranfile import FortranFile >>> f=FortranFile('uxuyp.bin', endian='<') >>> uxuy = f.readReals(prec='f') # 'f' for default reals >>> len(uxuy) 400 >>> ux = np.array(uxuy[:200]).reshape((20,10)).T >>> uy = np.array(uxuy[200:]).reshape((20,10)).T >>> p = f.readReals('f').reshape((20,10)).T >>> ux array([[ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20.], [ 2., 4., 6., 8., 10., 12., 14., 16., 18., 20., 22., 24., 26., 28., 30., 32., 34., 36., 38., 40.], [ 3., 6., 9., 12., 15., 18., 21., 24., 27., 30., 33., 36., 39., 42., 45., 48., 51., 54., 57., 60.], [ 4., 8., 12., 16., 20., 24., 28., 32., 36., 40., 44., 48., 52., 56., 60., 64., 68., 72., 76., 80.], [ 5., 10., 15., 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100.], [ 6., 12., 18., 24., 30., 36., 42., 48., 54., 60., 66., 72., 78., 84., 90., 96., 102., 108., 114., 120.], [ 7.,Proxy-Connection: keep-alive Cache-Control: max-age=0 14., 21., 28., 35., 42., 49., 56., 63., 70., 77., 84., 91., 98., 105., 112., 119., 126., 133., 140.], [ 8., 16., 24., 32., 40., 48., 56., 64., 72., 80., 88., 96., 104., 112., 120., 128., 136., 144., 152., 160.], [ 9., 18., 27., 36., 45., 54., 63., 72., 81., 90., 99., 108., 117., 126., 135., 144., 153., 162., 171., 180.], [ 10., 20., 30., 40., 50., 60., 70., 80., 90., 100., 110., 120., 130., 140., 150., 160., 170., 180., 190., 200.]]) >>> uy array([[ 1. , 0.5 , 0.33333334, 0.25 , 0.2 , 0.16666667, 0.14285715, 0.125 , 0.11111111, 0.1 , 0.09090909, 0.08333334, 0.07692308, 0.07142857, 0.06666667, 0.0625 , 0.05882353, 0.05555556, 0.05263158, 0.05 ], [ 2. , 1. , 0.66666669, 0.5 , 0.40000001, 0.33333334, 0.2857143 , 0.25 , 0.22222222, 0.2 , 0.18181819, 0.16666667, 0.15384616, 0.14285715, 0.13333334, 0.125 , 0.11764706, 0.11111111, 0.10526316, 0.1 ], [ 3. , 1.5 , 1. , 0.75 , 0.60000002, 0.5 , 0.42857143, 0.375 , 0.33333334, 0.30000001, 0.27272728, 0.25 , 0.23076923, 0.21428572, 0.2 , 0.1875 , 0.17647059, 0.16666667, 0.15789473, 0.15000001], [ 4. , 2. , 1.33333337, 1. , 0.80000001, 0.66666669, 0.5714286 , 0.5 , 0.44444445, 0.40000001, 0.36363637, 0.33333334, 0.30769232, 0.2857143 , 0.26666668, 0.25 , 0.23529412, 0.22222222, 0.21052632, 0.2 ], [ 5. , 2.5 , 1.66666663, 1.25 , 1. , 0.83333331, 0.71428573, 0.625 , 0.55555558, 0.5 , 0.45454547, 0.41666666, 0.38461539, 0.35714287, 0.33333334, 0.3125 , 0.29411766, 0.27777779, 0.2631579 , 0.25 ], [ 6. , 3. , 2. , 1.5 , 1.20000005, 1. , 0.85714287, 0.75 , 0.66666669, 0.60000002, 0.54545456, 0.5 , 0.46153846, 0.42857143, 0.40000001, 0.375 , 0.35294119, 0.33333334, 0.31578946, 0.30000001], [ 7. , 3.5 , 2.33333325, 1.75 , 1.39999998, 1.16666663, 1. , 0.875 , 0.77777779, 0.69999999, 0.63636363, 0.58333331, 0.53846157, 0.5 , 0.46666667, 0.4375 , 0.41176471, 0.3888889 , 0.36842105, 0.34999999], [ 8. , 4. , 2.66666675, 2. , 1.60000002, 1.33333337, 1.14285719, 1. , 0.8888889 , 0.80000001, 0.72727275, 0.66666669, 0.61538464, 0.5714286 , 0.53333336, 0.5 , 0.47058824, 0.44444445, 0.42105263, 0.40000001], [ 9. , 4.5 , 3. , 2.25 , 1.79999995, 1.5 , 1.28571427, 1.125 , 1. , 0.89999998, 0.81818181, 0.75 , 0.69230771, 0.64285713, 0.60000002, 0.5625 , 0.52941179, 0.5 , 0.47368422, 0.44999999], [ 10. , 5. , 3.33333325, 2.5 , 2. , 1.66666663, 1.42857146, 1.25 , 1.11111116, 1. , 0.90909094, 0.83333331, 0.76923078, 0.71428573, 0.66666669, 0.625 , 0.58823532, 0.55555558, 0.52631581, 0.5 ]]) >>> p array([[ 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21.], [ 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22.], [ 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23.], [ 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24.], [ 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.], [ 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26.], [ 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27.], [ 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28.], [ 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28., 29.], [ 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28., 29., 30.]]) Note that you have to provide the shape information for ux and uy because fortran writes them together as a stream of 400 numbers. -Neil

10 10

Extending C with Python
by Jialin Liu 07 Feb '18

07 Feb '18

Hello, I'm extending C with python (which is opposite way of what people usually do, extending python with C), I'm currently stuck in passing a C array to python layer, could anyone plz advise? I have a C buffer in my C code and want to pass it to a python function. In the C code, I have: npy_intp dims [2]; > dims[0] = 10; > dims[1] = 20; > import_array(); > npy_intp m=2; > PyObject * py_dims = PyArray_SimpleNewFromData(1, &m, NPY_INT16 ,(void > *)dims ); // I also tried NPY_INT > PyObject_CallMethod(pInstance, method_name, "O", py_dims); In the Python code, I want to just print that array: def f(self, dims): print ("np array:%d,%d"%(dims[0],dims[1])) But it only prints the first number correctly, i.e., dims[0]. The second number is always 0. Best, Jialin LBNL/NERSC

5 7

documentation, sound and music [was: doc? music through mathematical relations between LPCM samples and musical elements/characteristics]
by Renato Fabbri 01 Feb '18

01 Feb '18

Dear Scipy-ers, If you think I should split the message so that things get more clear... But the things are: 1) How to document numpy-heavy projects? 2) How to better make these contributions available to the numpy/scipy community? Directions will be greatly appreciated. I suspect that this info is all gathered somewhere I did not find. PS. I apologize for the cross-post, this is a message sent to numpy list but without any replies, so I found fit to send the message to both lists. Thanks in advance, Best, R. ---------- Forwarded message ---------- From: Renato Fabbri <renato.fabbri(a)gmail.com> Date: Tue, Jan 30, 2018 at 8:32 PM Subject: doc? music through mathematical relations between LPCM samples and musical elements/characteristics To: Discussion of Numerical Python <numpy-discussion(a)python.org> the half-shape suite: archive.org/details/ShapeSuite was completely synthesized using psychophysical relations for each resulting 16bit 44kHz samples: https://arxiv.org/abs/1412.6853 I am thinking about the ways in which to make the documentation at least to: * mass (music and audio in sample sequences): the framework. It might be considered a toolbox, but it is not a package: https://github.com/ttm/mass/ * music: a package for making music, including the routines used to make the half-shape suite: https://github.com/ttm/music == It seems reasonable at first to: * upload the package to PyPI (mass is not a package), music is there (beta but there). * Make available some summary of the file tree, including some automated code documentation. As I follow numpy's convention (or try to), and the package and framework are particularly useful for proficient numpy users, I used Sphinx. The very preliminary version is: https://pythonmusic.github.io/ * Make a nice readthedocs documentation. "music" project name was taken so I made: http://music-documentation.readthedocs.io/en/latest/ "mass" project name was also taken so I made: http://musicmass.readthedocs.io/en/latest/ And I found these URLs clumsy. Is there a standard or would you go with musicpackage.read.. and massframework.read... ? More importantly: * would you use readthedocs for the sphinx/doxigen output? * would you use readthedocs.org or a gitbook would be better or...? Should I contact the scipy community to make available a scikit or integrate it to numpy in any way beyond using and citing it appropriately?. BTW. I am using vimwiki with some constant attempt to organize and track my I/O: https://arxiv.org/abs/1712.06933 So maybe a unified solution for all these documentation instances using a wiki structure seems reasonable at the moment. Maybe upload the generated html from Sphinx and from Vimwiki to readthedocs...? Best, R. -- Renato Fabbri GNU/Linux User #479299 labmacambira.sourceforge.net -- Renato Fabbri GNU/Linux User #479299 labmacambira.sourceforge.net

2 2

custom Welch method for power spectral density
by Seb 31 Jan '18

31 Jan '18

Hello, I'm trying to compute a power spectral density of a signal, using the Welch method, in the broad sense; i.e. splitting the signal into segments for deriving smoother spectra. This is well implemented in scipy.signal.welch. However, I'd like to use exponentially increasing (power 2) segment length to dampen increasing variance in spectra at higher frequencies. Before hacking the scipy.signal.spectral module for this, I'd appreciate any tips on available packages/modules that allow for this kind of binning scheme, or other suggestions. Thanks, -- Seb

5 9

Setting custom dtypes and 1.14
by Chris Barker 31 Jan '18

31 Jan '18

Hi all, I'm pretty sure this is the same thing as recently discussed on this list about 1.14, but to confirm: I had failures in my code with an upgrade for 1.14 -- turns out it was a single line in a single test fixture, so no big deal, but a regression just the same, with no deprecation warning. I was essentially doing this: In [*48*]: dt Out[*48*]: dtype([('time', '<i8'), ('value', [('u', '<f8'), ('v', '<f8')])], align=True) In [*49*]: uv Out[*49*]: array([[1., 1.], [1., 1.], [1., 1.], [1., 1.]]) In [*50*]: time Out[*50*]: array([1, 1, 1, 1]) In [*51*]: full = np.array(zip(time, uv), dtype=dt) --------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-51-ed726f71dd4a> in <module>() ----> 1 full = np.array(zip(time, uv), dtype=dt) ValueError: setting an array element with a sequence. It took some poking, but the solution was to do: full = np.array(zip(time, (tuple(w) *for* w *in* uv)), dtype=dt) That is, convert the values to nested tuples, rather than an array in a tuple, or a list in a tuple. As I said, my problem is solved, but to confirm: 1) This is a known change with good reason? 2) My solution was the best (only) one -- the only way to set a nested dtype like that is with tuples? If so, then I think we should: A) improve the error message. "ValueError: setting an array element with a sequence." Is not really clear -- I spent a while trying to figure out how I could set a nested dtype like that without a sequence? and I was actually using a ndarray, so it wasn't even a generic sequence. And a tuple is a sequence, too... I had a vague recollection that in some circumstances, numpy treats tuples and lists (and arrays) differently (fancy indexing??), so I tried the tuple thing and that worked. But I've been around numpy a long time -- that could have been very very confusing to many people. So could the message be changed to something like: "ValueError: setting an array element with a generic sequence. Only the tuple type can be used in this context." or something like that -- I'm not sure where else this same error message might pop up, so that could be totally inappropriate. 2) maybe add a .totuple()method to ndarray, much like the .tolist() method? that would have been handy here. -Chris -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/OR&R (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker(a)noaa.gov

7 34

doc? music through mathematical relations between LPCM samples and musical elements/characteristics
by Renato Fabbri 31 Jan '18

31 Jan '18

the half-shape suite: archive.org/details/ShapeSuite was completely synthesized using psychophysical relations for each resulting 16bit 44kHz samples: https://arxiv.org/abs/1412.6853 I am thinking about the ways in which to make the documentation at least to: * mass (music and audio in sample sequences): the framework. It might be considered a toolbox, but it is not a package: https://github.com/ttm/mass/ * music: a package for making music, including the routines used to make the half-shape suite: https://github.com/ttm/music == It seems reasonable at first to: * upload the package to PyPI (mass is not a package), music is there (beta but there). * Make available some summary of the file tree, including some automated code documentation. As I follow numpy's convention (or try to), and the package and framework are particularly useful for proficient numpy users, I used Sphinx. The very preliminary version is: https://pythonmusic.github.io/ * Make a nice readthedocs documentation. "music" project name was taken so I made: http://music-documentation.readthedocs.io/en/latest/ "mass" project name was also taken so I made: http://musicmass.readthedocs.io/en/latest/ And I found these URLs clumsy. Is there a standard or would you go with musicpackage.read.. and massframework.read... ? More importantly: * would you use readthedocs for the sphinx/doxigen output? * would you use readthedocs.org or a gitbook would be better or...? Should I contact the scipy community to make available a scikit or integrate it to numpy in any way beyond using and citing it appropriately?. BTW. I am using vimwiki with some constant attempt to organize and track my I/O: https://arxiv.org/abs/1712.06933 So maybe a unified solution for all these documentation instances using a wiki structure seems reasonable at the moment. Maybe upload the generated html from Sphinx and from Vimwiki to readthedocs...? Best, R. -- Renato Fabbri GNU/Linux User #479299 labmacambira.sourceforge.net

1 0

Optimisation of matrix multiplication
by Andrew Nelson 30 Jan '18

30 Jan '18

Hi all, I have a matrix multiplication that I'd like to optimize. I have a matrix `a` (dtype=complex) with shape (N, M, 2, 2). I'd like to do the following multiplication: a[:, 0] @ a[:, 1] @ ... @ a[:, M-1] where the first dimension, N, is element wise (and hopefully vectorisable) and M>=2. So for each row do M-1 matrix multiplications of 2x2 matrices. The output should have shape (N, 2, 2). What would be the best way of going about this? -- _____________________________________ Dr. Andrew Nelson _____________________________________

1 1

Re: [Numpy-discussion] Moving NumPy's PRNG Forward
by Kevin Sheppard 30 Jan '18

30 Jan '18

I agree with pretty much everything you wrote Robert. I didn't have quote the right frame but the generic class that takes a low-level core PRNG sounds like the right design, and this should make user-generated distributions easier to develop. I was thinking along these lines inspired by the SpiPy changes that use a LowLevelCallable, e.g., https://ilovesymposia.com/2017/03/12/scipys-new-lowlevelcallable-is-a-game-… This might also allow users to extend the core PRNGs using something like Numba JIT classes as an alternative. Another area that needs though is how to correctly spawn in Multiprocess application. This might be most easily addressed by providing a guide on a good way rather than the arbitrary way used now. Kevin On Sat, Jan 27, 2018 at 5:03 PM <numpy-discussion-request(a)python.org> wrote: > Send NumPy-Discussion mailing list submissions to > numpy-discussion(a)python.org > > To subscribe or unsubscribe via the World Wide Web, visit > https://mail.python.org/mailman/listinfo/numpy-discussion > or, via email, send a message with subject or body 'help' to > numpy-discussion-request(a)python.org > > You can reach the person managing the list at > numpy-discussion-owner(a)python.org > > When replying, please edit your Subject line so it is more specific > than "Re: Contents of NumPy-Discussion digest..." > > > Today's Topics: > > 1. Re: Moving NumPy's PRNG Forward (Robert Kern) > 2. Re: Using np.frombuffer and cffi.buffer on array of C structs > (problem with struct member padding) (Joe) > > > ---------------------------------------------------------------------- > > Message: 1 > Date: Sat, 27 Jan 2018 09:28:54 +0900 > From: Robert Kern <robert.kern(a)gmail.com> > To: Discussion of Numerical Python <numpy-discussion(a)python.org> > Subject: Re: [Numpy-discussion] Moving NumPy's PRNG Forward > Message-ID: > <CAF6FJitFLv3U7gHkYBCFW69A5BbVe=HzBAm5oxYVcXGbBdMU= > w(a)mail.gmail.com> > Content-Type: text/plain; charset="utf-8" > > On Sat, Jan 27, 2018 at 1:14 AM, Kevin Sheppard < > kevin.k.sheppard(a)gmail.com> > wrote: > > > > I am a firm believer that the current situation is not sustainable. > There are a lot of improvements that can practically be incorporated. > While many of these are performance related, there are also improvements in > accuracy over some ranges of parameters that cannot be incorporated. I also > think that perfect stream reproducibility is a bit of a myth across > versions since this really would require identical OS, compiler and > possibly CPU for some of the generators that produce floats. > > > > I believe there is a case for separating the random generator from core > NumPy. Some points that favor becoming a subproject: > > > > 1. It is a pure consumer of NumPy API. Other parts of the API do no > depend on random. > > 2. A stand alone package could be installed along side many different > version of core NumPy which would reduce the pressure on freezing the > stream. > > Removing numpy.random (or freezing it as deprecated legacy while all PRNG > development moves elsewhere) is probably a non-starter. It's too used for > us not to provide something. That said, we can (and ought to) make it much > easier for external packages to provide PRNG capabilities (core PRNGs and > distributions) that interoperate with the core functionality that numpy > provides. I'm also happy to place a high barrier on adding more > distributions to numpy.random once that is in place. > > Specifically, core uniform PRNGs should have a small common C API that > distribution functions can use. This might just be a struct with an opaque > `void*` state pointer and then 2 function pointers for drawing a uint64 > (whole range) and a double in [0,1) from the state. It's important to > expose our core uniform PRNGs as a C API because there has been a desire to > interoperate at that level, using the same PRNG state inside C or Fortran > or GPU code. If that's in place, then people can write new efficient > distribution functions in C that use this small C API agnostic to the core > PRNG algorithm. It also makes it easy to implement new core PRNGs that the > distribution functions provided by numpy.random can use. > > > In terms of what is needed, I think that the underlying PRNG should be > swappable. The will provide a simple mechanism to allow certain types of > advancement while easily providing backward compat. In the current design > this is very hard and requires compiling many nearly identical copies of > RandomState. In pseudocode something like > > > > standard_normal(prng) > > > > where prng is a basic class that retains the PRNG state and has a small > set of core random number generators that belong to the underlying PRNG -- > probably something like int32, int64, double, and possibly int53. I am not > advocating explicitly passing the PRNG as an argument, but having > generators which can take any suitable PRNG would add a lot of flexibility > in terms of taking advantage of improvements in the underlying PRNGs (see, > e.g., xoroshiro128/xorshift1024). The "small" core PRNG would have > responsibility over state and streams. The remainder of the module would > transform the underlying PRNG into the required distributions. > > (edit: after writing the following verbiage, I realize it can be summed up > with more respect to your suggestion: yes, we should do this design, but we > don't need to and shouldn't give up on a class with distribution methods.) > > Once the core PRNG C API is in place, I don't think we necessarily need to > move away from a class structure per se, though it becomes an option. We > just separate the core PRNG object from the distribution-providing class. > We don't need to make copies of the distribution-providing class just to > use a new core PRNG. I'm coming around to Nathaniel's suggestion for the > constructor API (though not the distribution-versioning, for reasons I can > get into later). We have a couple of core uniform PRNG classes like > `MT19937` and `PCG128`. Those have a tiny API, and probably don't have a > lot of unnecessary code clones between them. Their constructors can be > different depending on the different ways they can be instantiated, > depending on the PRNG's features. I'm not sure that they'll have any common > methods besides `__getstate__/__setstate__` and probably a `copy()`. They > will expose their C API as a Python-opaque attribute. They can have > whatever algorithm-dependent methods they need (e.g. to support jumpahead). > I might not even expose to Python the uint64 and U(0,1) double sampling > methods, but maybe so. > > Then we have a single `Distributions` class that provides all of the > distributions that we want to support in numpy.random (i.e. what we > currently have on `RandomState` and whatever passes our higher bar in the > future). It takes one of the core PRNG instances as an argument to the > constructor (nominally, at least; we can design factory functions to make > this more convenient). > > prng = Distributions(PCG128(seed)) > x = prng.normal(mean, std) > > If someone wants to write a WELL512 core PRNG, they can just implement that > object and pass it to `Distributions()`. The `Distributions` code doesn't > need to be copied, nor do we need to much around with `__new__` tricks in > Cython. > > Why do this instead of distribution functions? Well, I have a damn lot of > code that is expecting an object with at least the broad outlines of the > `RandomState` interface. I'm not going to update that code to use functions > instead. And if I'm not, no one is. There isn't a good transition path > since the PRNG object needs to thread through all of the code, whether it's > my code that I'm writing greenfield or library code that I don't control. > That said, people writing new distributions outside of numpy.random would > be writing functions, not trying to add to the `Distributions` class, but > that's fine. > > It also allows the possibility for the `Distributions` class to be stateful > if we want to do things like caching the next Box-Muller variate and not > force that onto the core PRNG state like I currently do. Though I'd rather > just drop Box-Muller, and that's not a common pattern outside of > Box-Muller. But it's a possibility. > > -- > Robert Kern >

1 0

Re: [Numpy-discussion] Moving NumPy's PRNG Forward
by Kevin Sheppard 30 Jan '18

30 Jan '18

I am a firm believer that the current situation is not sustainable. There are a lot of improvements that can practically be incorporated. While many of these are performance related, there are also improvements in accuracy over some ranges of parameters that cannot be incorporated. I also think that perfect stream reproducibility is a bit of a myth across versions since this really would require identical OS, compiler and possibly CPU for some of the generators that produce floats. I believe there is a case for separating the random generator from core NumPy. Some points that favor becoming a subproject: 1. It is a pure consumer of NumPy API. Other parts of the API do no depend on random. 2. A stand alone package could be installed along side many different version of core NumPy which would reduce the pressure on freezing the stream. In terms of what is needed, I think that the underlying PRNG should be swappable. The will provide a simple mechanism to allow certain types of advancement while easily providing backward compat. In the current design this is very hard and requires compiling many nearly identical copies of RandomState. In pseudocode something like standard_normal(prng) where prng is a basic class that retains the PRNG state and has a small set of core random number generators that belong to the underlying PRNG -- probably something like int32, int64, double, and possibly int53. I am not advocating explicitly passing the PRNG as an argument, but having generators which can take any suitable PRNG would add a lot of flexibility in terms of taking advantage of improvements in the underlying PRNGs (see, e.g., xoroshiro128/xorshift1024). The "small" core PRNG would have responsibility over state and streams. The remainder of the module would transform the underlying PRNG into the required distributions. This would also simplify making improvements, since old versions could be saved or improved versions could be added to the API. For example, from numpy.random import standard_normal, prng # Preferred versions standard_normal(prng) # Ziggurat from numpy.random.legacy import standard_normal_bm, mt19937 # legacy generators standard_normal_bm(mt19937) # Box-Muller Kevin

3 3

NumPy-Discussion January 2018