Mailman 3 September 2008 - NumPy-Discussion

ANN: NumPy 1.2.0
by Jarrod Millman Oct. 17, 2008

Oct. 17, 2008

I'm pleased to announce the release of NumPy 1.2.0. NumPy is the fundamental package needed for scientific computing with Python. It contains: * a powerful N-dimensional array object * sophisticated (broadcasting) functions * basic linear algebra functions * basic Fourier transforms * sophisticated random number capabilities * tools for integrating Fortran code. Besides it's obvious scientific uses, NumPy can also be used as an efficient multi-dimensional container of generic data. Arbitrary data-types can be defined. This allows NumPy to seamlessly and speedily integrate with a wide-variety of databases. This minor release comes almost four months after the 1.1.0 release. The major features of this release are a new testing framework and huge amount of documentation work. It also includes a some minor API breakage scheduled in the 1.1 release. Please note that NumPy 1.2.0 requires Python 2.4 or greater. For information, please see the release notes: http://sourceforge.net/project/shownotes.php?group_id=1369&release_id=628858 You can download the release from here: http://sourceforge.net/project/showfiles.php?group_id=1369&package_id=175103 Thank you to everybody who contributed to this release. Enjoy, -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/

11 43

Proposal: scipy.spatial
by Anne Archibald Oct. 12, 2008

Oct. 12, 2008

Hi, Once again there has been a thread on the numpy/scipy mailing lists requesting (essentially) some form of spatial data structure. Pointers have been posted to ANN (sadly LGPLed and in C++) as well as a handful of pure-python implementations of kd-trees. I suggest the creation of a new submodule of scipy, scipy.spatial, to contain spatial data structures and algorithms. Specifically, I propose it contain a kd-tree implementation providing nearest-neighbor, approximate nearest-neighbor, and all-points-near queries. I have a few other suggestions for things it might contain, but kd-trees seem like a good first step. 2008/9/27 Nathan Bell <wnbell(a)gmail.com>: > On Sat, Sep 27, 2008 at 11:18 PM, Anne Archibald > <peridot.faceted(a)gmail.com> wrote: >> >> I think a kd-tree implementation would be a valuable addition to >> scipy, perhaps in a submodule scipy.spatial that might eventually >> contain other spatial data structures and algorithms. What do you >> think? Should we have one? Should it be based on Sturla Molden's code, >> if the license permits? I am willing to contribute one, if not. > > +1 Judging that your vote and mine are enough in the absence of dissenting voices, I have written an implementation based on yours, Sturla Molden's, and the algorithms described by the authors of the ANN library. Before integrating it into scipy, though, I'd like to send it around for comments. Particular issues: * It's pure python right now, but with some effort it could be partially or completely cythonized. This is probably a good idea in the long run. In the meantime I have crudely vectorized it so that users can at least avoid looping in their own code. * It is able to return the r nearest neighbors, optionally subject to a maximum distance limit, as well as approximate nearest neighbors. * It is not currently set up to conveniently return all points within some fixed distance of the target point, but this could easily be added. * It returns distances and indices of nearest neighbors in the original array. * The testing code is, frankly, a mess. I need to look into using nose in a sensible fashion. * The license is the scipy license. I am particularly concerned about providing a convenient return format. The natural return from a query is a list of neighbors, since it may have variable length (there may not be r elements in the tree, or you might have supplied a maximum distance which doesn't contain r points). For a single query, it's simple to return a python list (should it be sorted? currently it's a heap). But if you want to vectorize the process, storing the results in an array becomes cumbersome. One option is an object array full of lists; another, which, I currently use, is an array with nonexistent points marked with an infinite distance and an invalid index. A third would be to return masked arrays. How do you recommend handling this variable (or potentially-variable) sized output? > If you're implementing one, I would highly recommend the > "left-balanced" kd-tree. > http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/2535/pdf/imm2535.pdf Research suggests that at least in high dimension a more geometric balancing criterion can produce vastly better results. I used the "sliding midpoint" rule, which doesn't allow such a numerical balancing but does ensure that you don't have too many long skinny cells (since a sphere tends to cut very many of these). I've also been thinking about what else would go in scipy.spatial. I think it would be valuable to have a reasonably efficient distance matrix function (we seem to get that question a lot, and the answer's not trivial) as well as a sparse distance matrix function based on the kd-trees. The module could potentially also swallow the (currently sandboxed?) delaunay code. Anne

11 22

Generating random samples without repeats
by Paul Moore Oct. 3, 2008

Oct. 3, 2008

I want to generate a series of random samples, to do simulations based on them. Essentially, I want to be able to produce a SAMPLESIZE * N matrix, where each row of N values consists of either 1. Integers between 1 and M (simulating M rolls of an N-sided die), or 2. A sample of N numbers between 1 and M without repeats (simulating deals of N cards from an M-card deck). Example (1) is easy, numpy.random.random_integers(1, M, (SAMPLESIZE, N)) But I can't find an obvious equivalent for (2). Am I missing something glaringly obvious? I'm using numpy - is there maybe something in scipy I should be looking at? Also, in evaluating samples, I'm likely to want to calculate combinatorial functions, such as the list of all pairs of items from a sample (imagine looking at how many pairs add up to 15 in a cribbage hand). Clearly, I can write a normal Python function which does this for one row, and use apply_along_axis - but that's *slow*. I'm looking for a function that, given an N*M array and a sample size S, gives a C(N,S)*S*M array of all the combinations, which runs at array-processing speeds (preferably without having to code it in C myself!!) Is there anywhere with this type of function available? This type of combinatorial simulation seems to me to be a fairly good fit for numpy's capabilities, and yet I can't seem to find things that seem relevant. Is it simly not something that people use numpy for? Or am I looking in the wrong places in the documentation? Thanks for any help, Paul.

16 40

1.2.0rc2 tagged! --PLEASE TEST--
by Jarrod Millman Oct. 3, 2008

Oct. 3, 2008

Hello, The 1.2.0rc2 is now available: http://svn.scipy.org/svn/numpy/tags/1.2.0rc2 The source tarball is here: https://cirl.berkeley.edu/numpy/numpy-1.2.0rc2.tar.gz Here is the universal Mac binary: https://cirl.berkeley.edu/numpy/numpy-1.2.0rc2-py2.5-macosx10.5.dmg Here are the Window's binaries: http://www.ar.media.kyoto-u.ac.jp/members/david/archives/numpy/numpy-1.2.0r… http://www.ar.media.kyoto-u.ac.jp/members/david/archives/numpy/numpy-1.2.0r… Here are the release notes: http://scipy.org/scipy/numpy/milestone/1.2.0 Please test this release ASAP and let us know if there are any problems. If there are no show stoppers, this will become the 1.2.0 release. Thanks, -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/

13 27

numpy.random.logseries - incorrect convergence for k=1, k=2
by joep Oct. 3, 2008

Oct. 3, 2008

random numbers generated by numpy.random.logseries do not converge to theoretical distribution: for probability paramater pr = 0.8, the random number generator converges to a frequency for k=1 at 39.8 %, while the theoretical probability mass is 49.71 k=2 is oversampled, other k's look ok check frequency of k=1 and k=2 at N = 1000000 0.398406 0.296465 pmf at k = 1 and k=2 with formula [ 0.4971 0.1988] for probability paramater pr = 0.3, the results are not as bad, but still off: frequency for k=1 at 82.6 %, while the theoretical probability mass is 84.11 check frequency of k=1 and k=2 at N = 1000000 0.826006 0.141244 pmf at k = 1 and k=2 with formula [ 0.8411 0.1262] below is a quick script for checking this Josef {{{ import numpy as np from scipy import stats pr = 0.8 np.set_printoptions(precision=2, suppress=True) # calculation for N=1million takes some time for N in [1000, 10000, 10000, 1000000]: rvsn=np.random.logseries(pr,size=N) fr=stats.itemfreq(rvsn) pmfs=stats.logser.pmf(fr[:,0],pr)*100 print 'log series sample frequency and pmf (in %) with N = ', N print np.column_stack((fr[:,0],fr[:,1]*100.0/N,pmfs)) np.set_printoptions(precision=4, suppress=True) print 'check frequency of k=1 and k=2 at N = ', N print np.sum(rvsn==1)/float(N), print np.sum(rvsn==2)/float(N) k = np.array([1,2]) print 'pmf at k = 1 and k=2 with formula' print -pr**k * 1.0 / k / np.log(1-pr) }}}

1 1

numpy.random.hypergeometric - strange results
by joep Oct. 3, 2008

Oct. 3, 2008

In my fuzz testing of scipy stats, I get sometimes a test failure. I think there is something wrong with numpy.random.hypergeometric for some cases: Josef >>> import numpy.random as mtrand >>> mtrand.hypergeometric(3,17,12,size=10) # there are only 3 good balls in urn array([16, 17, 16, 16, 15, 16, 17, 16, 17, 16]) >>> mtrand.hypergeometric(17,3,12,size=10) #negative result array([-3, -4, -3, -4, -3, -3, -4, -4, -5, -4]) >>> mtrand.hypergeometric(4,3,12,size=10) >>> np.version.version '1.2.0rc2' I did not find any clear pattern when trying out different parameter values: >>> mtrand.hypergeometric(10,10,12,size=10) array([5, 6, 4, 4, 8, 5, 4, 6, 7, 4]) >>> mtrand.hypergeometric(10,10,20,size=10) array([10, 10, 10, 10, 10, 10, 10, 10, 10, 10]) >>> mtrand.hypergeometric(10,10,19,size=10) array([10, 9, 9, 9, 9, 9, 10, 9, 9, 9]) >>> mtrand.hypergeometric(10,10,5,size=10) array([3, 5, 2, 2, 1, 2, 2, 4, 3, 1]) >>> mtrand.hypergeometric(10,2,5,size=10) array([4, 5, 4, 5, 5, 5, 4, 3, 4, 4]) >>> mtrand.hypergeometric(2,10,5,size=10) array([0, 2, 1, 0, 2, 2, 1, 1, 1, 1]) >>> mtrand.hypergeometric(17,3,12,size=10) array([-5, -3, -4, -4, -4, -3, -4, -4, -3, -3]) >>> mtrand.hypergeometric(3,17,12,size=10) array([15, 16, 17, 16, 15, 16, 15, 15, 17, 17]) >>> mtrand.hypergeometric(18,3,12,size=10) array([-5, -6, -6, -4, -4, -4, -5, -3, -5, -5]) >>> mtrand.hypergeometric(18,3,5,size=10) array([4, 5, 5, 5, 5, 5, 4, 5, 4, 3]) >>> mtrand.hypergeometric(18,3,19,size=10) array([1, 1, 2, 1, 1, 1, 1, 3, 1, 1])

1 1

PyArray_Resize reference counting
by Ravi Oct. 1, 2008

Oct. 1, 2008

Hi, Consider the following code: PyArrayObject *array = get_me_an_array_along_with_a_reference(); /* 1 */ PyArray_Dims *dims = get_me_some_dims(); /* 2 */ array = (PyArrayObject *)PyArray_Resize( array, dims, 1 ); /* 3 */ Is the code above valid? In particular, is the assigment on line 3 valid? Travis' book says that a reference to the new array is returned. Does that mean that I need to store a pointer to the original array somewhere else and decrement its reference count after the resize call? Alternatively, if I don't need to store the pointer because the returned pointer is the same as the one passed into the resize function, do I need to decrement the reference count since I already have a reference to the array? In other words, if the reference count is, say, 16 before line 1 and 17 after line 1, will it be 18 after line 3? In that case, I need to decrement the reference once to get it back to 17 because I should not need a reference as I already have one. Regards, Ravi

2 4

xml-rpc with numpy arrays
by Brian Blais Oct. 1, 2008

Oct. 1, 2008

Hello, I am trying to use xml-rpc to be able to run some simulations remotely. I am running into a problem with the transfer of numpy arrays. my server code looks like: #!/usr/bin/env python def again(x): # test out the sending of data return [x,x] from SimpleXMLRPCServer import SimpleXMLRPCServer SimpleXMLRPCServer.allow_reuse_address = 1 server = SimpleXMLRPCServer(("", 8000)) server.register_function(again) try: print "Serving..." server.serve_forever() # Start the server finally: print "done." server.server_close() my client code looks like: import numpy from xmlrpclib import ServerProxy server=ServerProxy('http://localhost:8000') server.again(5) # this works b=numpy.random.rand(5,5) server.again(b) # this gives an error this gives the error: <type 'exceptions.TypeError'>: cannot marshal <type 'numpy.ndarray'> objects which seems to be a deficiency of the marshal library, or perhaps I am doing something wrong. Is there a way to fix this? Is there another approach that I should be using? thanks, Brian Blais -- Brian Blais bblais(a)bryant.edu http://web.bryant.edu/~bblais

6 10

Upper and lower envelopes
by bevan Oct. 1, 2008

Oct. 1, 2008

Hello, I have some XY data. I would like to generate the equations for an upper and lower envelope that excludes a percentage of the data points. I would like to define the slope of the envelope line (say 3) and then have my code find the intercept that fits my requirements (say 5% of data below the lower envelope). This would then give me the equation and I could plot the upper and lower envelopes. I hope this makes sense. Thanks for any help.

4 4

Re: [Numpy-discussion] Help needed with the priority of 0d-arrays and np.scalars
by Pierre GM Sept. 30, 2008

Sept. 30, 2008

Sorry to bump my own post, I know it's rude... However, I'm in dire need for some pointers: what are the priority rules when multiplying numpy scalars and 0d ndarrays ? Thanks a lot in advance

2 2