Mailman 3 September 2013 - NumPy-Discussion

mixed mode arithmetic
by Neal Becker July 11, 2023

July 11, 2023

I've been browsing the numpy source. I'm wondering about mixed-mode arithmetic on arrays. I believe the way numpy handles this is that it never does mixed arithmetic, but instead converts arrays to a common type. Arguably, that might be efficient for a mix of say, double and float. Maybe not. But for a mix of complex and a scalar type (say, CDouble * Double), it's clearly suboptimal in efficiency. So, do I understand this correctly? If so, is that something we should improve?

4 6

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

Feb. 18, 2023

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 July 22, 2021

July 22, 2021

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

C-API: multidimensional array indexing?
by Johann Bauer March 31, 2016

March 31, 2016

Dear experts, is there a C-API function for numpy which implements Python's multidimensional indexing? Say, I have a 2d-array PyArrayObject * M; and an index int i; how do I extract the i-th row or column M[i,:] respectively M[:,i]? I am looking for a function which gives again a PyArrayObject * and which is a view to M (no copied data; the result should be another PyArrayObject whose data and strides points to the correct memory portion of M). I searched the API documentation, Google and mailing lists for quite a long time but didn't find anything. Can you help me? Thanks, Johann

3 3

taking a 2D uneven surface slice
by Bryan Woods Dec. 3, 2014

Dec. 3, 2014

I'm trying to do something that at first glance I think should be simple but I can't quite figure out how to do it. The problem is as follows: I have a 3D grid Values[Nx, Ny, Nz] I want to slice Values at a 2D surface in the Z dimension specified by Z_index[Nx, Ny] and return a 2D slice[Nx, Ny]. It is not as simple as Values[:,:,Z_index]. I tried this: >>> values.shape (4, 5, 6) >>> coords.shape (4, 5) >>> slice = values[:,:,coords] >>> slice.shape (4, 5, 4, 5) >>> slice = np.take(values, coords, axis=2) >>> slice.shape (4, 5, 4, 5) >>> Obviously I could create an empty 2D slice and then fill it by using np.ndenumerate to fill it point by point by selecting values[i, j, Z_index[i, j]]. This just seems too inefficient and not very pythonic.

3 3

Speedup by avoiding memory alloc twice in scalar array
by Arink Verma Jan. 8, 2014

Jan. 8, 2014

Hi, I am working on performance parity between numpy scalar/small array and python array as GSOC mentored By Charles. Currently I am looking at PyArray_Return, which allocate separate memory just for scalar return. Unlike python which allocate memory once for returning result of scalar operations; numpy calls malloc twice once for the array object itself, and a second time for the array data. These memory allocations are happening in PyArray_NewFromDescr and PyArray_Scalar. Stashing both within a single allocation would be more efficient. In, PyArray_Scalar, new struct (PyLongScalarObject) need allocation in case of scalar arrays. Instead, can we just some how convert/cast PyArrayObject to PyLongScalarObject.?? -- Arink Verma www.arinkverma.in

5 12

1.8.0rc1
by Charles R Harris Oct. 2, 2013

Oct. 2, 2013

Hi All, NumPy 1.8.0rc1 is up now on sourceforge<http://sourceforge.net/projects/numpy/files/NumPy/1.8.0rc1/>.The binary builds are included except for Python 3.3 on windows, which will arrive later. Many thanks to Ralf for the binaries, and to those who found and fixed the bugs in the last beta. Any remaining bugs are all my fault ;) I hope this will be the last release before final, so please test it thoroughly. Chuck

10 31

[ANN] MATLAB fmincon now available in Python2
by Dmitrey Sept. 30, 2013

Sept. 30, 2013

Hi all, current state of Python <-> MATLAB connection soft doesn't allow passing of function handlers, however, a walkaround has been implemented via some tricks, so now MATLAB function fmincon is available in Python-written OpenOpt and FuncDesigner frameworks (with possibility of automatic differentiation, example ). Future plans include MATLAB fsolve, ode23, ode45 (unlike scipy fsolve and ode they can handle sparse matrices), fgoalattain, maybe global optimization toolbox solvers. I intend to post the message to several forums, so to keep discussion in a single place use OpenOpt forum thread http://forum.openopt.org/viewtopic.php?id=769 Regards, D.

1 0

Re: [Numpy-discussion] Baffling error: ndarray += csc_matrix -> "ValueError: setting an array element with a sequence"
by Pauli Virtanen Sept. 27, 2013

Sept. 27, 2013

27.09.2013 19:33, Nathaniel Smith kirjoitti: [clip] > I really don't understand what arcane magic is used to make ndarray += > csc_matrix work at all, but my question is, is it going to break when > we complete the casting transition described above? It was just > supposed to catch things like int += float. This maybe clarifies it: >>> import numpy >>> import scipy.sparse >>> x = numpy.ones((2,2)) >>> y = scipy.sparse.csr_matrix(x) >>> z = x >>> z += y >>> x array([[ 1., 1.], [ 1., 1.]]) >>> z matrix([[ 2., 2.], [ 2., 2.]]) The execution flows like this: ndarray.__iadd__(arr, sparr) np.add(arr, sparr, out=???) return NotImplemented # wtf return NotImplemented Python does arr = sparr.__radd__(arr) Since Scipy master sparse matrices now have __numpy_ufunc__, but it doesn't handle out= arguments, the second step currently raises a TypeError (for Numpy master + Scipy master). And this is actually the correct thing to do, as having np.add return NotImplemented is just broken. Only ndarray.__iadd__ has the authority to return the NotImplemented. To make the in-place ops work again, it seems Numpy needs some additional fixes in its binary op machinery, before __numpy_ufunc__ business works fully as intended. Namely, the binary op routines will need to catch TypeErrors and convert them to NotImplemented. The code paths where Numpy ufuncs currently return NotImplemented could also be changed to raise TypeErrors, but I'm not sure if someone somewhere relies on this behavior (I hope not). -- Pauli Virtanen

2 3

Indexing changes/deprecations
by Sebastian Berg Sept. 27, 2013

Sept. 27, 2013

Hey, since I am working on the indexing. I was wondering about a few smaller things: * 0-d boolean array, `np.array(0)[True]` (will work now) would give np.array([0]) as a copy, instead of the original array. I guess I could add a FutureWarning or so, but I am not sure and overall the chance of creating bugs seems low. (The boolean index should always add 1 dimension and here, remove 0 dimensions -> 1-d result.) * All index operations return a view; never the object. This means that `v = arr[...]` is slightly slower. But since it does not affect `arr[...] = vals`, I think the speed implications are negligible. * Does anyone have an idea if there is a way to change the subclass logic that view based item setting is implemented as: np.asarray(subclass[index]) = vals I somewhat think the subclass should rather implement `__setitem__` instead of relying on numpy calling its `__getitem__`, but I don't see how it can be changed. * Still thinking a bit about implementing a keepdims keyword or function, to handle matrix type logic mostly in the C-code. And most importantly, is there any behaviour thing in the index machinery that is bugging you, which I may have forgotten until now? - Sebastian

4 5

Re: [Numpy-discussion] Valid algorithm for generating a 3D Wiener Process?
by David Goldsmith Sept. 27, 2013

Sept. 27, 2013

Thanks, guys. Yeah, I realized the problem w/ the uniform-increment-variable-direction approach this morning: physically, it ignores the fact that the particles hitting the particle being tracked are going to have a distribution of momentum, not all the same, just varying in direction. But I don't quite understand Warren's observation: "the 'angles' that describe the position undergo a random walk [actually, it would seem that they don't, since they too fail the varying-as-white-noise test], so the particle tends to move in the same direction over short intervals"--is this just another way of saying that, since I was varying the angles by -1, 0, or 1 unit each time, the simulation is susceptible to "unnaturally" long strings of -1, 0, or 1 increments? Thanks again, DG

4 4

Re: [Numpy-discussion] invalid correlation coefficient from np.ma.corrcoef
by Faraz Mirzaei Sept. 26, 2013

Sept. 26, 2013

Thanks Josef and Nathaniel for your responses. In the application that I have, I don't use the correlation coefficient matrix as a whole (so I don't care if it is PSD or not). I simply read the off-diagonal elements for pair-wise correlation coefficients. I use the pairwise correlation coefficient to test if the data from various sources (i.e., rows of the matrix), agree with each other when present. Right now, I use, ma.corrcoef( x[ i, :] , x[ j, :]) and read the off-diagonal element in a loop over i and j. It is just a bit uglier than calling ma.corrcoef(x). At least for my application, truncation to -1 or +1 (or scaling such that largest values becomes 1 etc) is completely wrong, since it would imply that the two sources completely agree with each other (factoring out a minus sign), which may not the case. For example, consider the first and last rows of the example I provided: >>> print x_ma [[ 7 -4 -1 -- -3 -2] [ 6 -3 0 4 0 5] [-4 -- 7 5 -- --] [-- 5 -- 0 1 4]] >>> np.ma.corrcoef(x_ma)[0,3] -1.6813456149534147 On the other hand, if we supply only the first and third row to the function, we get: >>> np.ma.corrcoef(x_ma[0,:], x_ma[3,:]) masked_array(data = [[1.0 -0.240192230708] [-0.240192230708 1.0]], mask = [[False False] [False False]], fill_value = 1e+20) Interestingly, this is the same as what pandas results as the [3,0] element of the correlation coefficient matrix, and it is equal to pair-wise deletion result: >>> np.corrcoef([-4, -3, -2], [5, 1, 4]) #Note that this is NOT np.ma.corrcoef >>> array([[ 1. , -0.24019223], [-0.24019223, 1. ]]) Also, I don't know why the ma.corrcoef results Josef has mentioned are different than mine. In particular, Josef reports element [2, 0] of the ma.corrcoef result to be -1.19, but I get -- (i.e., missing and masked, probably due to too few samples available). Josef: are you sure that you have entered the example values correctly into python? Along the same lines, the results that Nathaniel has posted from R are different, since the input is not a masked matrix I guess (please note that in the original example, I had masked values less than or equal to -5). In any case, I think the correlation coefficient between two rows of a matrix should not depend on what other rows are supplied. In other words, np.ma.corrcoef(x_ma)[0,3] should be equal to np.ma.corrcoef(x_ma[0,:], x_ma[3,:])[0,1] (which apparently happens to be what pandas reports). This change would need recomputing the mean for every pair-wise coefficient calculation, but since we are computing cross products O(n^2) times, the overall big-O complexity won't change. And please don't remove this functionality. I will volunteer to fix it however we decide :) We can just clarify the behavior in the documentation. Thanks, Faraz

3 2

At.: use less RAM memory and increase execution speed
by Josè Luis Mietta Sept. 26, 2013

Sept. 26, 2013

Hi experts! I wanna use less RAM memory in my Monte Carlo simulations. In my algorithm I use numpy arrays and xrange() function. I hear that I can reduce RAM used in my lagorithm if I do the next: 1) replace xrange() for range(). 2) replace numpya arrays for python lists 3) use reset() function for deleting useless arrays. Is that true? In adition, I wanna increase execution speed of my code (I use numpy and SciPy functions). How can I apply Cython? Will it help? Please help. Thanks a lot!!

3 2

invalid correlation coefficient from np.ma.corrcoef
by Faraz Mirzaei Sept. 26, 2013

Sept. 26, 2013

Hi everyone, I'm using np.ma.corrcoef to compute the correlation coefficients among rows of a masked matrix, where the masked elements are the missing data. I've observed that in some cases, the np.ma.corrcoef gives invalid coefficients that are greater than 1 or less than -1. Here's an example: x = array([[ 7, -4, -1, -7, -3, -2], [ 6, -3, 0, 4, 0, 5], [-4, -5, 7, 5, -7, -7], [-5, 5, -8, 0, 1, 4]]) x_ma = np.ma.masked_less_equal(x , -5) C = np.round(np.ma.corrcoef(x_ma), 2) print C [[1.0 0.73 -- -1.68] [0.73 1.0 -0.86 -0.38] [-- -0.86 1.0 --] [-1.68 -0.38 -- 1.0]] As you can see, the [0,3] element is -1.68 which is not a valid correlation coefficient. (Valid correlation coefficients should be between -1 and 1). I looked at the code for np.ma.corrcoef, and this behavior seems to be due to the way that mean values of the rows of the input matrix are computed and subtracted from them. Apparently, the mean value is individually computed for each row, without masking the elements corresponding to the masked elements of the other row of the matrix, with respect to which the correlation coefficient is being computed. I guess the right way should be to recompute the mean value for each row every time that a correlation coefficient is being computed for two rows after propagating pair-wise masked values. Please let me know what you think. Thanks, Faraz

3 8

Valid algorithm for generating a 3D Wiener Process?
by David Goldsmith Sept. 25, 2013

Sept. 25, 2013

Is this a valid algorithm for generating a 3D Wiener process? (When I graph the results, they certainly look like potential Brownian motion tracks.) def Wiener3D(incr, N): r = incr*(R.randint(3, size=(N,))-1) r[0] = 0 r = r.cumsum() t = 2*np.pi*incr*(R.randint(3, size=(N,))-1) t[0] = 0 t = t.cumsum() p = np.pi*incr*(R.randint(3, size=(N,))-1) p[0] = 0 p = p.cumsum() x = r*np.cos(t)*np.sin(p) y = r*np.sin(t)*np.sin(p) z = r*np.cos(p) return np.array((x,y,z)).T Thanks! DG

3 3

Heads up re: new memory tracking infrastructure in python 3.4
by Nathaniel Smith Sept. 25, 2013

Sept. 25, 2013

Hi all, Possibly of interest to some folks here, it looks like in Python 3.4 there'll be a generic system for hooking and tracking memory allocations: http://www.python.org/dev/peps/pep-0445/ I'm not planning to do anything with this information myself, but if anyone's been thinking about numpy memory profiling, the new numpy malloc hooks, etc., then it might make sense to see how we could integrate with this more general system. -n

1 0

Numpy newbie question: matrix creation
by Edmondo Porcu Sept. 25, 2013

Sept. 25, 2013

Dear all, I am a Newbie with Numpy and I would need some advice on how to create a matrix with certain characteristics : - Every entry should be minimum 0 maximum 1 with a step of 0.1 (legal values are 0,0.1,0.2,0.3 etc) - The number of columns of the matrix is a parameter of this matrix creation algorithm - Only the rows where the sum is 1 must be kept Would great appreciate your advice and suggestions Best Regards Edmondo

4 5

Removal of numarray and oldnumeric packages.
by Charles R Harris Sept. 24, 2013

Sept. 24, 2013

Hi All, I have gotten no feedback on the removal of the numarray and oldnumeric packages. Consequently the removal will take place on 9/28. Scream now or never... Chuck

5 9

Pull request review #3770: Trapezoidal distribution
by Jeremy Hetzel Sept. 24, 2013

Sept. 24, 2013

I've added a trapezoidal distribution to numpy.random for consideration, pull request 3770: https://github.com/numpy/numpy/pull/3770 Similar to the triangular distribution, the trapezoidal distribution may be used where the underlying distribution is not known, but some knowledge of the limits and mode exists. The trapezoidal distribution generalizes the triangular distribution by allowing the modal values to be expressed as a range instead of a point estimate. The trapezoidal distribution implemented, known as the "generalized trapezoidal distribution," has three additional parameters: growth, decay, and boundary ratio. Adjusting these from the default values create trapezoidal-like distributions with non-linear behavior. Examples can be seen in an R vignette ( http://cran.r-project.org/web/packages/trapezoid/vignettes/trapezoid.pdf ), as well as these papers by J.R. van Dorp and colleagues: 1) van Dorp, J. R. and Kotz, S. (2003) Generalized trapezoidal distributions. Metrika. 58(1):85–97. Preprint available: http://www.seas.gwu.edu/~dorpjr/Publications/JournalPapers/Metrika2003VanDo… 2) van Dorp, J. R., Rambaud, S.C., Perez, J. G., and Pleguezuelo, R. H. (2007) An elicitation procedure for the generalized trapezoidal distribution with a uniform central stage. Decision Analysis Journal. 4:156–166. Preprint available: http://www.seas.gwu.edu/~dorpjr/Publications/JournalPapers/DA2007.pdf The docstring for the proposed numpy.random.trapezoidal() is as follows: """ trapezoidal(left, mode1, mode2, right, size=None, m=2, n=2, alpha=1) Draw samples from the generalized trapezoidal distribution. The trapezoidal distribution is defined by minimum (``left``), lower mode (``mode1``), upper mode (``mode1``), and maximum (``right``) parameters. The generalized trapezoidal distribution adds three more parameters: the growth rate (``m``), decay rate (``n``), and boundary ratio (``alpha``) parameters. The generalized trapezoidal distribution simplifies to the trapezoidal distribution when ``m = n = 2`` and ``alpha = 1``. It further simplifies to a triangular distribution when ``mode1 == mode2``. Parameters ---------- left : scalar Lower limit. mode1 : scalar The value where the first peak of the distribution occurs. The value should fulfill the condition ``left <= mode1 <= mode2``. mode2 : scalar The value where the first peak of the distribution occurs. The value should fulfill the condition ``mode1 <= mode2 <= right``. right : scalar Upper limit, should be larger than or equal to `mode2`. size : int or tuple of ints, optional Output shape. Default is None, in which case a single value is returned. m : scalar, optional Growth parameter. n : scalar, optional Decay parameter. alpha : scalar, optional Boundary ratio parameter. Returns ------- samples : ndarray or scalar The returned samples all lie in the interval [left, right]. Notes ----- With ``left``, ``mode1``, ``mode2``, ``right``, ``m``, ``n``, and ``alpha`` parametrized as :math:`a, b, c, d, m, n, \\text{ and } \\alpha`, respectively, the probability density function for the generalized trapezoidal distribution is .. math:: f{\\scriptscriptstyle X}(x\mid\theta) = \\mathcal{C}(\\Theta) \\times \\begin{cases} \\alpha \\left(\\frac{x - \\alpha}{b - \\alpha} \\right)^{m - 1}, & \\text{for } a \\leq x < b \\\\ (1 - \\alpha) \\left(\frac{x - b}{c - b} \\right) + \\alpha, & \\text{for } b \\leq x < c \\\\ \\left(\\frac{d - x}{d - c} \\right)^{n-1}, & \\text{for } c \\leq x \\leq d \\end{cases} with the normalizing constant :math:`\\mathcal{C}(\\Theta)` defined as ..math:: \\mathcal{C}(\\Theta) = \\frac{2mn} {2 \\alpha \\left(b - a\\right) n + \\left(\\alpha + 1 \\right) \\left(c - b \\right)mn + 2 \\left(d - c \\right)m} and where the parameter vector :math:`\\Theta = \\{a, b, c, d, m, n, \\alpha \\}, \\text{ } a \\leq b \\leq c \\leq d, \\text{ and } m, n, \\alpha >0`. Similar to the triangular distribution, the trapezoidal distribution may be used where the underlying distribution is not known, but some knowledge of the limits and mode exists. The trapezoidal distribution generalizes the triangular distribution by allowing the modal values to be expressed as a range instead of a point estimate. The growth, decay, and boundary ratio parameters of the generalized trapezoidal distribution further allow for non-linear behavior to be specified. References ---------- .. [1] van Dorp, J. R. and Kotz, S. (2003) Generalized trapezoidal distributions. Metrika. 58(1):85–97. Preprint available: http://www.seas.gwu.edu/~dorpjr/Publications/JournalPapers/Metrika2003VanDo… .. [2] van Dorp, J. R., Rambaud, S.C., Perez, J. G., and Pleguezuelo, R. H. (2007) An elicitation proce-dure for the generalized trapezoidal distribution with a uniform central stage. Decision AnalysisJournal. 4:156–166. Preprint available: http://www.seas.gwu.edu/~dorpjr/Publications/JournalPapers/DA2007.pdf Examples -------- Draw values from the distribution and plot the histogram: >>> import matplotlib.pyplot as plt >>> h = plt.hist(np.random.triangular(0, 0.25, 0.75, 1, 100000), bins=200, ... normed=True) >>> plt.show() """ I am unsure if NumPy encourages incorporation of new distributions into numpy.random or instead into separate modules, but found the exercise to be helpful regardless. Thanks, Jeremy

3 4

Binary releases
by Charles R Harris Sept. 23, 2013

Sept. 23, 2013

Hi All, Numpy 1.8 is about ready for an rc1, which brings up the question of which binary builds so put up on sourceforge. For Windows maybe 1. 32 bit windows, python 2.6, 2.7, 3.2, 3.3, compiled with MSVC 2. 64 bit windows, python 2.6, 2.7, 3.2, 3.3, compiled with MSVC, linked with MKL I think these should be good for both windows 7 and window 8, correct me if I am wrong. For Mac there is first the question of OS X versions, (10.5?), 10.6, 10.7, 10.8. I don't know if some builds work on more than one OS X version. The 10.5 version is a bit harder to come by than 10.6 and up. It looks like 10.9 is coming up, but it isn't out yet. I have no idea what Python version to match up these, but assuming all of them, then 1. OS X 10.6 python 2.6, 2.7, 3.2, 3.3, compiled with native compiler, linked with Accelerate. 2. OS X 10.7 python 2.6, 2.7, 3.2, 3.3, compiled with native compiler, linked with Accelerate. 3. OS X 10.8 python 2.6, 2.7, 3.2, 3.3, compiled with native compiler, linked with Accelerate. That seems like a lot. It is fairly easy to compile from source on the mac these days, are all those binary packages really needed? I don't know what I am doing with the binary stuff, so any suggestions are welcome. For building it is possible to set up a Macintosh with vmware fusion to handle all of these, but there is time and money involved in that. Any one who is already capable of doing some of these builds is welcome to volunteer. Note, Intel has offered MKL licenses to the open source projects under the NumFocus umbrella, but I don't know of anyone who has taken advantage of that yet or how much time it will take to go through the undoubtedly needed paper work, but I would like to get some of those licenses and move to MSVC so that we can stop rolling around gasping in pain whenever it comes to window builds. Then there is Fortran for windows... Thoughts? Chuck

11 31

Generating a (uniformly distributed) random bit list of length N
by David Goldsmith Sept. 23, 2013

Sept. 23, 2013

Thanks, St?fan, speed: N ~ 1e9. Thanks again. DG ---------------------------------------------------------------------- > > Message: 1 > Date: Sun, 22 Sep 2013 14:04:09 -0700 > From: David Goldsmith <d.l.goldsmith(a)gmail.com> > Subject: [Numpy-discussion] Generating a (uniformly distributed) > random bit list of length N > To: numpy-discussion(a)scipy.org > Message-ID: > < > CAFtPsZqG7Upjy8s04NPwkn8gJVDVGORu3rPjCVLN6ZKRWnPuTQ(a)mail.gmail.com> > Content-Type: text/plain; charset="iso-8859-1" > > Is np.random.randint(2, size=N) the fastest way to do this? Thanks! > > DG >

1 0

openblas?
by Neal Becker Sept. 23, 2013

Sept. 23, 2013

Does numpy/scipy support building with openblas for blas,lapack instead of atlas?

2 1

Runtime Warning by numpy.divide(0.0, 0.0) can't be ignored
by Lorenzo Di Gregorio Sept. 23, 2013

Sept. 23, 2013

Hello, numpy.divide(0.0,0.0) raises a warning which doesn't appear to get caught by numpy.seterr(divide='ignore') This example: import numpy numpy.divide(0.0,0.0) numpy.divide(1.0,0.0) numpy.seterr(divide='ignore') numpy.divide(1.0,0.0) numpy.divide(0.0,0.0) produces this output: Warning (from warnings module): File "C:\Users\ldigrego\Desktop\tst.py", line 2 numpy.divide(0.0,0.0) RuntimeWarning: invalid value encountered in divide Warning (from warnings module): File "C:\Users\ldigrego\Desktop\tst.py", line 3 numpy.divide(1.0,0.0) RuntimeWarning: divide by zero encountered in divide Warning (from warnings module): File "C:\Users\ldigrego\Desktop\tst.py", line 6 numpy.divide(0.0,0.0) RuntimeWarning: invalid value encountered in divide Strangely, on Ipython the numpy.seterr(divide='ignore') seems to catch the warning: import numpy numpy.divide(0.0,0.0) Out[1]: nan C:\eclipse\plugins\org.python.pydev_2.8.1.2013072611\pysrc\pydevconsole.py:1: RuntimeWarning: invalid value encountered in divide from __future__ import nested_scopes #Jython 2.1 support numpy.divide(1.0,0.0) Out[1]: inf C:\eclipse\plugins\org.python.pydev_2.8.1.2013072611\pysrc\pydevconsole.py:1: RuntimeWarning: divide by zero encountered in divide from __future__ import nested_scopes #Jython 2.1 support numpy.seterr(divide='ignore') Out[1]: {'divide': 'warn', 'invalid': 'warn', 'over': 'warn', 'under': 'ignore'} numpy.divide(1.0,0.0) Out[1]: inf numpy.divide(0.0,0.0) Out[1]: nan I could not find information on Google: is this a known problem? Is there a way to suppress this warning? I'm working on a 64b Win7 machine employing numpy-MKL-1.7.0.win-amd64-py2.7.exe. Best Regards, Lorenzo

2 1

Generating a (uniformly distributed) random bit list of length N
by David Goldsmith Sept. 22, 2013

Sept. 22, 2013

Is np.random.randint(2, size=N) the fastest way to do this? Thanks! DG

2 1

Min Error
by David Reed Sept. 22, 2013

Sept. 22, 2013

Hi, I am getting a strange error when finding the minimum of a matrix. The weird thing is I get this while running within iPython shell, and if I do %debug and go to the line where this fails and run the command `a = np.min(D, axis=0)`, I get no error. Here is the trace: draw_lines/main.pyc in haus_distance(A, B) 28 D = D.reshape((len(A), None)) 29 ---> 30 a = np.min(D, axis=0) 31 32 # vp = np.max(np.min(D, axis=0)) anaconda/lib/python2.7/site-packages/numpy/core/fromnumeric.pyc in amin(a, axis, out, keepdims) 1985 else: 1986 return _methods._amin(a, axis=axis, -> 1987 out=out, keepdims=keepdims) 1988 1989 def alen(a): anaconda/lib/python2.7/site-packages/numpy/core/_methods.pyc in _amin(a, axis, out, keepdims) 12 def _amin(a, axis=None, out=None, keepdims=False): 13 return um.minimum.reduce(a, axis=axis, ---> 14 out=out, keepdims=keepdims) 15 16 def _sum(a, axis=None, dtype=None, out=None, keepdims=False): TypeError: an integer is required

2 4