Re: [Numpy-discussion] slicing and aliasing
On Thu, May 31, 2012 at 7:30 AM, Neal Becker ndbeck...@gmail.com wrote: That is, will: u[a:b] = u[c:d] always work (assuming the ranges of a:b, d:d are equal, or course) It works most of the time. This thread shows you how to find an example where it does not work: http://mail.scipy.org/pipermail/numpy-discussion/2012-May/062079.html ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] record arrays initialization
On Wed, May 2, 2012 at 4:46 PM, Kevin Jacobs jac...@bioinformed.com bioinfor...@gmail.com wrote: The cKDTree implementation is more than 4 times faster than the brute-force approach: T = scipy.spatial.cKDTree(targets) In [11]: %timeit foo1(element, targets) # Brute force 1000 loops, best of 3: 385 us per loop In [12]: %timeit foo2(element, T) # cKDTree 1 loops, best of 3: 83.5 us per loop In [13]: 385/83.5 Out[13]: 4.610778443113772 Brute force plus cython beats cKDTree for knn with k=1 and Euclidean distance: I[5] targets = np.random.uniform(0, 8, (5000, 7)) I[6] element = np.arange(1, 8, dtype=np.float64) I[7] T = scipy.spatial.cKDTree(targets) I[8] timeit T.query(element) 1 loops, best of 3: 36.1 us per loop I[9] timeit nn(targets, element) 1 loops, best of 3: 28.5 us per loop What about lower dimensions (2 instead of 7) where cKDTree gets faster? I[18] element = np.arange(1,3, dtype=np.float64) I[19] targets = np.random.uniform(0,8,(5000,2)) I[20] T = scipy.spatial.cKDTree(targets) I[21] timeit T.query(element) 1 loops, best of 3: 27.5 us per loop I[22] timeit nn(targets, element) 10 loops, best of 3: 11.6 us per loop I could add nn to the bottleneck package. Is the k=1, Euclidean distance case too specialized? Prototype (messy) code is here: https://github.com/kwgoodman/bottleneck/issues/45 For a smaller speed up (~2x) foo1 could use bottleneck's sum or squares function, bn.ss(). ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] timing results (was: record arrays initialization)
On Thu, May 3, 2012 at 10:38 AM, Moroney, Catherine M (388D) catherine.m.moro...@jpl.nasa.gov wrote: Actually Fortran with correct array ordering - 13 seconds! What horrible python/numpy mistake am I making to cause such a slowdown? For the type of problem you are working on, I'd flip it around and ask what you are doing with Fortran to only get 5x ;) Well, I'd ask that if I knew Fortran. I'm seeing 7.5x with cython/numpy. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] timing results (was: record arrays initialization)
On Thu, May 3, 2012 at 12:46 PM, Paul Anton Letnes paul.anton.let...@gmail.com wrote: Could you show us the code? It's hard to tell otherwise. As Keith Goodman pointed out, if he gets 7.5x with cython, it could be that the Fortran code could be improved as well. Fortran has a reputation of being the gold standard for performance in numerical computation, although one can often be surprised. Picking good algorithms is always more important than the language. Doing array operations with Fortran might slow things down. I got the speed up by calculating the distance one point at a time and not storing the results (only keeping the min value and index). That gets rid of the sorting (argmin) at the end too. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] timing results (was: record arrays initialization)
On Thu, May 3, 2012 at 3:12 PM, Moroney, Catherine M (388D) catherine.m.moro...@jpl.nasa.gov wrote: Here is the python code: def single(element, targets): if (isinstance(element, tuple)): xelement = element[0] elif (isinstance(element, numpy.ndarray)): xelement = element else: return FILL nlen = xelement.shape[0] nvec = targets.data.shape[0] x = xelement.reshape(1, nlen).repeat(nvec, axis=0) repeat is slow. I don't think you need it since broadcasting should take care of things. (But maybe I misunderstand the data.) diffs = ((x - targets.data)**2).sum(axis=1) You could try np.dot instead of sum: one = np.ones(7); diff = np.dot(diff, one). You could even pass one in. diffs = numpy.sqrt(diffs) Since you don't return the distance, no need for the sqrt. If you do need the sqrt only take the sqrt of one element instead of all elements. return int(numpy.argmin(diffs, axis=0)) ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] all elements equal
On Mon, Mar 5, 2012 at 11:14 AM, Neal Becker ndbeck...@gmail.com wrote: What is a simple, efficient way to determine if all elements in an array (in my case, 1D) are equal? How about close? For the exactly equal case, how about: I[1] a = np.array([1,1,1,1]) I[2] np.unique(a).size O[2] 1# All equal I[3] a = np.array([1,1,1,2]) I[4] np.unique(a).size O[4] 2 # All not equal ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] all elements equal
On Mon, Mar 5, 2012 at 11:36 AM, josef.p...@gmail.com wrote: How about numpy.ptp, to follow this line? I would expect it's single pass, but wouldn't short circuit compared to cython of Keith I[1] a = np.ones(10) I[2] timeit (a == a[0]).all() 1000 loops, best of 3: 203 us per loop I[3] timeit a.min() == a.max() 1 loops, best of 3: 106 us per loop I[4] timeit np.ptp(a) 1 loops, best of 3: 106 us per loop I[5] a[1] = 9 I[6] timeit (a == a[0]).all() 1 loops, best of 3: 89.7 us per loop I[7] timeit a.min() == a.max() 1 loops, best of 3: 102 us per loop I[8] timeit np.ptp(a) 1 loops, best of 3: 103 us per loop ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] all elements equal
On Mon, Mar 5, 2012 at 11:52 AM, Benjamin Root ben.r...@ou.edu wrote: Another issue to watch out for is if the array is empty. Technically speaking, that should be True, but some of the solutions offered so far would fail in this case. Good point. For fun, here's the speed of a simple cython allclose: I[2] a = np.ones(10) I[3] timeit a.min() == a.max() 1 loops, best of 3: 106 us per loop I[4] timeit allequal(a) 1 loops, best of 3: 68.9 us per loop I[5] a[1] = 9 I[6] timeit a.min() == a.max() 1 loops, best of 3: 102 us per loop I[7] timeit allequal(a) 100 loops, best of 3: 269 ns per loop where @cython.boundscheck(False) @cython.wraparound(False) def allequal(np.ndarray[np.float64_t, ndim=1] a): cdef: np.float64_t a0 Py_ssize_t i, n=a.size a0 = a[0] for i in range(n): if a[i] != a0: return False return True ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] all elements equal
On Mon, Mar 5, 2012 at 12:06 PM, Neal Becker ndbeck...@gmail.com wrote: But doesn't this one fail on empty array? Yes. I'm optimizing for fun, not for corner cases. This should work for size zero and NaNs: @cython.boundscheck(False) @cython.wraparound(False) def allequal(np.ndarray[np.float64_t, ndim=1] a): cdef: np.float64_t a0 Py_ssize_t i, n=a.size if n == 0: return False # Or would you like True? a0 = a[0] for i in range(n): if a[i] != a0: return False return True ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] all elements equal
On Mon, Mar 5, 2012 at 12:12 PM, Keith Goodman kwgood...@gmail.com wrote: On Mon, Mar 5, 2012 at 12:06 PM, Neal Becker ndbeck...@gmail.com wrote: But doesn't this one fail on empty array? Yes. I'm optimizing for fun, not for corner cases. This should work for size zero and NaNs: @cython.boundscheck(False) @cython.wraparound(False) def allequal(np.ndarray[np.float64_t, ndim=1] a): cdef: np.float64_t a0 Py_ssize_t i, n=a.size if n == 0: return False # Or would you like True? a0 = a[0] for i in range(n): if a[i] != a0: return False return True Sorry for all the posts. I'll go back to being quiet. Seems like np.allclose returns True for empty arrays: I[2] a = np.array([]) I[3] np.allclose(np.array([]), np.array([])) O[3] True The original allequal cython code did the same: I[4] allequal(a) O[4] True ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] speed of array creation from tuples
On Mon, Feb 27, 2012 at 2:55 PM, Skipper Seabold jsseab...@gmail.com wrote: I am surprised by this (though maybe I shouldn't be?) It's always faster to use list comprehension to unpack lists of tuples than np.array/asarray? [~/] [1]: X = [tuple(np.random.randint(10,size=2)) for _ in range(100)] [~/] [2]: timeit np.array([x1 for _,x1 in X]) 1 loops, best of 3: 26.4 us per loop [~/] [3]: timeit np.array(X) 1000 loops, best of 3: 378 us per loop [~/] [4]: timeit x1, x2 = np.array([x for _,x in X]), np.array([y for y,_ in X]) 1 loops, best of 3: 53.4 us per loop [~/] [5]: timeit x1, x2 = np.array(X).T 1000 loops, best of 3: 384 us per loop Here's a similar thread: http://mail.scipy.org/pipermail/numpy-discussion/2010-February/048304.html ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Creating a bool array with Cython
On Sat, Feb 25, 2012 at 7:04 PM, Dag Sverre Seljebotn d.s.seljeb...@astro.uio.no wrote: On 02/25/2012 03:26 PM, Keith Goodman wrote: Is this a reasonable (and fast) way to create a bool array in cython? def makebool(): cdef: int n = 2 np.npy_intp *dims = [n] np.ndarray[np.uint8_t, ndim=1] a a = PyArray_EMPTY(1, dims, NPY_UINT8, 0) a[0] = 1 a[1] = 0 a.dtype = np.bool return a Will a numpy bool array be np.uint8 on all platforms? Someone else will have to answer that for sure, though I expect and always assumed so. How can I do a.dtype=np.bool using the numpy C api? Is there any point (speed) in doing so? Did you try np.ndarray[np.uint8_t, ndim=1, cast=True] a ? You should be able to do that and then pass NPY_BOOL to PyArray_EMPTY, to avoid having to change the dtype later. That looks like the grown up way to write the code. (Plus saves 300 ns in overhead.) Thanks, Dag. (I don't write the code to create the empty output array by hand. I have a templating system that does that. But it does not know how to create bool arrays. Hence the a.dtype = np.bool workaround. Thanks to your advice, I'll try to teach the templating system to create bool arrays so I don't have to remember how to do it.) ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Creating a bool array with Cython
Is this a reasonable (and fast) way to create a bool array in cython? def makebool(): cdef: int n = 2 np.npy_intp *dims = [n] np.ndarray[np.uint8_t, ndim=1] a a = PyArray_EMPTY(1, dims, NPY_UINT8, 0) a[0] = 1 a[1] = 0 a.dtype = np.bool return a Will a numpy bool array be np.uint8 on all platforms? How can I do a.dtype=np.bool using the numpy C api? Is there any point (speed) in doing so? ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] swaxes(0, 1) 10% faster than transpose on 2D matrix?
On Thu, Jan 19, 2012 at 1:37 AM, Mark Bakker mark...@gmail.com wrote: I noticed that swapaxes(0,1) is consistently (on my system) 10% faster than transpose on a 2D matrix. Transpose is faster for me. And a.T is faster than a.transpose() perhaps because a.transpose() checks that the inputs make sense? My guess is that they all do the same thing. It's just a matter of which function has the least overhead. I[10] a = np.random.rand(1000,1000) I[11] timeit a.T 1000 loops, best of 3: 153 ns per loop I[12] timeit a.transpose() 1000 loops, best of 3: 171 ns per loop I[13] timeit a.swapaxes(0,1) 100 loops, best of 3: 227 ns per loop I[14] timeit np.transpose(a) 100 loops, best of 3: 308 ns per loop ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] alterNEP - was: missing data discussion round 2
On Thu, Jun 30, 2011 at 10:51 AM, Nathaniel Smith n...@pobox.com wrote: On Thu, Jun 30, 2011 at 6:31 AM, Matthew Brett matthew.br...@gmail.com wrote: In the interest of making the discussion as concrete as possible, here is my draft of an alternative proposal for NAs and masking, based on Nathaniel's comments. Writing it, it seemed to me that Nathaniel is right, that the ideas become much clearer when the NA idea and the MASK idea are separate. Please do pitch in for things I may have missed or misunderstood: [...] Thanks for writing this up! I stuck it up as a gist so we can edit it more easily: https://gist.github.com/1056379/ This is your initial version: https://gist.github.com/1056379/c809715f4e9765db72908c605468304ea1eb2191 And I made a few changes: https://gist.github.com/1056379/33ba20300e1b72156c8fb655bd1ceef03f8a6583 Specifically, I added a rationale section, changed np.MASKED to np.IGNORE (as per comments in this thread), and added a vowel to propmsk. It might be helpful to make a small toy class in python so that people can play around with NA and IGNORE from the alterNEP. I only had a few minutes, so I only took it this far (1d arrays only): from nary import nary, NA, IGNORE arr = np.array([1,2,3,4,5,6]) nar = nary(arr) nar 1., 2., 3., 4., 5., 6., nar[2] = NA nar 1., 2., NA, 4., 5., 6., nar[4] = IGNORE nar 1., 2., NA, 4., IGNORE, 6., nar[4] IGNORE nar[3] 4 nar[2] NA The gist is here: https://gist.github.com/1057686 It probably just needs an __add__ and a reducing function such as sum, but I'm out of time, or so my family tells me. Implementation? Yes, with masks. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] feedback request: proposal to add masks to the core ndarray
On Thu, Jun 23, 2011 at 3:24 PM, Mark Wiebe mwwi...@gmail.com wrote: On Thu, Jun 23, 2011 at 5:05 PM, Keith Goodman kwgood...@gmail.com wrote: On Thu, Jun 23, 2011 at 1:53 PM, Mark Wiebe mwwi...@gmail.com wrote: Enthought has asked me to look into the missing data problem and how NumPy could treat it better. I've considered the different ideas of adding dtype variants with a special signal value and masked arrays, and concluded that adding masks to the core ndarray appears is the best way to deal with the problem in general. I've written a NEP that proposes a particular design, viewable here: https://github.com/m-paradox/numpy/blob/cmaskedarray/doc/neps/c-masked-array.rst There are some questions at the bottom of the NEP which definitely need discussion to find the best design choices. Please read, and let me know of all the errors and gaps you find in the document. Wow, that is exciting. I wonder about the relative performance of the two possible implementations (mask and NA) in the PEP. I've given that some thought, and I don't think there's a clear way to tell what the performance gap would be without implementations of both to benchmark against each other. I favor the mask primarily because it provides masking for all data types in one go with a single consistent interface to program against. For adding NA signal values, each new data type would need a lot of work to gain the same level of support. If you are, say, doing a calculation along the columns of a 2d array one element at a time, then you will need to grab an element from the array and grab the corresponding element from the mask. I assume the corresponding data and mask elements are not stored together. That would be slow since memory access is usually were time is spent. In this regard NA would be faster. Yes, the masks add more memory traffic and some extra calculation, while the NA signal values just require some additional calculations. I guess a better example would have been summing along rows instead of columns of a large C order array. If one needs to look at both the data and the mask then wouldn't summing along rows in cython be about as slow as it is currently to sum along columns? I currently use NaN as a missing data marker. That adds things like this to my cython code: if a[i] == a[i]: asum += a[i] If NA also had the property NA == NA is False, then it would be easy to use. That's what I believe it should do, and I guess this is a strike against the idea of returning None for a single missing value. If NA == NA is False then I wouldn't need to look at the mask in the example above. Or would ndarray have to look at the mask in order to return NA for a[i]? Which would mean __getitem__ would need to look at the mask? If the missing value is returned as a 0d array (so that NA == NA is False), would that break cython in a fundamental way since it could not always return a same-sized scalar when you index into an array? A mask, on the other hand, would be more difficult for third party packages to support. You have to check if the mask is present and if so do a mask-aware calculation; if is it not present then you have to do a non-mask based calculation. I actually see the mask as being easier for third party packages to support, particularly from C. Having regular C-friendly values with a boolean mask is a lot friendlier than values that require a lot of special casing like the NA signal values would require. So you have two code paths. You also need to check if any of the input arrays have masks and if so apply masks to the other inputs, etc. Most of the time, the masks will transparently propagate or not along with the arrays, with no effort required. In Python, the code you write would be virtually the same between the two approaches. -Mark ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] feedback request: proposal to add masks to the core ndarray
On Fri, Jun 24, 2011 at 7:06 AM, Robert Kern robert.k...@gmail.com wrote: The alternative proposal would be to add a few new dtypes that are NA-aware. E.g. an nafloat64 would reserve a particular NaN value (there are lots of different NaN bit patterns, we'd just reserve one) that would represent NA. An naint32 would probably reserve the most negative int32 value (like R does). Using the NA-aware dtypes signals that you are using NA values; there is no need for an additional flag. I don't understand the numpy design and maintainable issues, but from a user perspective (mine) nafloat64, etc sounds nice. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] feedback request: proposal to add masks to the core ndarray
On Thu, Jun 23, 2011 at 1:53 PM, Mark Wiebe mwwi...@gmail.com wrote: Enthought has asked me to look into the missing data problem and how NumPy could treat it better. I've considered the different ideas of adding dtype variants with a special signal value and masked arrays, and concluded that adding masks to the core ndarray appears is the best way to deal with the problem in general. I've written a NEP that proposes a particular design, viewable here: https://github.com/m-paradox/numpy/blob/cmaskedarray/doc/neps/c-masked-array.rst There are some questions at the bottom of the NEP which definitely need discussion to find the best design choices. Please read, and let me know of all the errors and gaps you find in the document. Wow, that is exciting. I wonder about the relative performance of the two possible implementations (mask and NA) in the PEP. If you are, say, doing a calculation along the columns of a 2d array one element at a time, then you will need to grab an element from the array and grab the corresponding element from the mask. I assume the corresponding data and mask elements are not stored together. That would be slow since memory access is usually were time is spent. In this regard NA would be faster. I currently use NaN as a missing data marker. That adds things like this to my cython code: if a[i] == a[i]: asum += a[i] If NA also had the property NA == NA is False, then it would be easy to use. A mask, on the other hand, would be more difficult for third party packages to support. You have to check if the mask is present and if so do a mask-aware calculation; if is it not present then you have to do a non-mask based calculation. So you have two code paths. You also need to check if any of the input arrays have masks and if so apply masks to the other inputs, etc. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] argmax for top N elements
On Wed, Jun 22, 2011 at 12:08 PM, RadimRehurek radimrehu...@seznam.cz wrote: Date: Wed, 22 Jun 2011 11:30:47 -0400 From: Alex Flint alex.fl...@gmail.com Subject: [Numpy-discussion] argmax for top N elements Is it possible to use argmax or something similar to find the locations of the largest N elements in a matrix? I would also be interested in an O(N) argmax/argmin for indices of top k values in an array. I'm currently using argsort[:k] in a performance sensitive part and it's not ideal. You can try argpartsort from the bottleneck package: a = np.random.rand(10) k = 10 timeit a.argsort()[:k] 100 loops, best of 3: 11.2 ms per loop import bottleneck as bn timeit bn.argpartsort(a, k)[:k] 1000 loops, best of 3: 1.29 ms per loop Output of argpartsort is not ordered: a.argsort()[:k] array([97239, 21091, 25130, 41638, 62323, 57419, 4381, 34905, 94572, 25935]) bn.argpartsort(a, k)[:k] array([21091, 62323, 25130, 97239, 41638, 57419, 4381, 34905, 94572, 25935]) ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] poor performance of sum with sub-machine-word integer types
On Tue, Jun 21, 2011 at 9:46 AM, Zachary Pincus zachary.pin...@yale.edu wrote:
Hello all,
As a result of the fast greyscale conversion thread, I noticed an anomaly
with numpy.ndararray.sum(): summing along certain axes is much slower with
sum() than versus doing it explicitly, but only with integer dtypes and when
the size of the dtype is less than the machine word. I checked in 32-bit and
64-bit modes and in both cases only once the dtype got as large as that did
the speed difference go away. See below...
Is this something to do with numpy or something inexorable about machine /
memory architecture?
Zach
Timings -- 64-bit mode:
--
In [2]: i = numpy.ones((1024,1024,4), numpy.int8)
In [3]: timeit i.sum(axis=-1)
10 loops, best of 3: 131 ms per loop
In [4]: timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
100 loops, best of 3: 2.57 ms per loop
In [5]: i = numpy.ones((1024,1024,4), numpy.int16)
In [6]: timeit i.sum(axis=-1)
10 loops, best of 3: 131 ms per loop
In [7]: timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
100 loops, best of 3: 4.75 ms per loop
In [8]: i = numpy.ones((1024,1024,4), numpy.int32)
In [9]: timeit i.sum(axis=-1)
10 loops, best of 3: 131 ms per loop
In [10]: timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
100 loops, best of 3: 6.37 ms per loop
In [11]: i = numpy.ones((1024,1024,4), numpy.int64)
In [12]: timeit i.sum(axis=-1)
100 loops, best of 3: 16.6 ms per loop
In [13]: timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
100 loops, best of 3: 15.1 ms per loop
Timings -- 32-bit mode:
--
In [2]: i = numpy.ones((1024,1024,4), numpy.int8)
In [3]: timeit i.sum(axis=-1)
10 loops, best of 3: 138 ms per loop
In [4]: timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
100 loops, best of 3: 3.68 ms per loop
In [5]: i = numpy.ones((1024,1024,4), numpy.int16)
In [6]: timeit i.sum(axis=-1)
10 loops, best of 3: 140 ms per loop
In [7]: timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
100 loops, best of 3: 4.17 ms per loop
In [8]: i = numpy.ones((1024,1024,4), numpy.int32)
In [9]: timeit i.sum(axis=-1)
10 loops, best of 3: 22.4 ms per loop
In [10]: timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
100 loops, best of 3: 12.2 ms per loop
In [11]: i = numpy.ones((1024,1024,4), numpy.int64)
In [12]: timeit i.sum(axis=-1)
10 loops, best of 3: 29.2 ms per loop
In [13]: timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
10 loops, best of 3: 23.8 ms per loop
One difference is that i.sum() changes the output dtype of int input
when the int dtype is less than the default int dtype:
i.dtype
dtype('int32')
i.sum(axis=-1).dtype
dtype('int64') # -- dtype changed
(i[...,0]+i[...,1]+i[...,2]+i[...,3]).dtype
dtype('int32')
Here are my timings
i = numpy.ones((1024,1024,4), numpy.int32)
timeit i.sum(axis=-1)
1 loops, best of 3: 278 ms per loop
timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
100 loops, best of 3: 12.1 ms per loop
import bottleneck as bn
timeit bn.func.nansum_3d_int32_axis2(i)
100 loops, best of 3: 8.27 ms per loop
Does making an extra copy of the input explain all of the speed
difference (is this what np.sum does internally?):
timeit i.astype(numpy.int64)
10 loops, best of 3: 29.2 ms per loop
No.
Initializing the output also adds some time:
timeit np.empty((1024,1024,4), dtype=np.int32)
10 loops, best of 3: 2.67 us per loop
timeit np.empty((1024,1024,4), dtype=np.int64)
10 loops, best of 3: 12.8 us per loop
Switching back and forth between the input and output array takes more
memory time too with int64 arrays compared to int32.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] fast SSD
On Tue, Jun 21, 2011 at 5:09 PM, Alex Flint alex.fl...@gmail.com wrote: Is there a fast way to compute an array of sum-of-squared-differences between a (small) K x K array and all K x K sub-arrays of a larger array? (i.e. each element x,y in the output array is the SSD between the small array and the sub-array (x:x+K, y:y+K) My current implementation loops over each sub-array and computes the SSD with something like ((A-B)**2).sum(). I don't know of a clever way. But if ((A-B)**2).sum() is a sizable fraction of the time, then you could use bottleneck: a = np.random.rand(5,5) timeit (a**2).sum() 10 loops, best of 3: 4.63 us per loop import bottleneck as bn timeit bn.ss(a) 100 loops, best of 3: 1.77 us per loop func, b = bn.func.ss_selector(a, axis=None) func built-in function ss_2d_float64_axisNone timeit func(b) 100 loops, best of 3: 830 ns per loop ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Returning the same dtype in Cython as np.argmax
On Wed, Jun 8, 2011 at 9:49 PM, Travis Oliphant oliph...@enthought.com wrote: On Jun 7, 2011, at 3:17 PM, Keith Goodman wrote: What is the rule to determine the dtype returned by numpy functions that return indices such as np.argmax? The return type of indices will be np.intp. Thanks, Travis. I now know a little more about NumPy. I don't know if Cython provides that type or not. If not, then it is either np.int32 or np.int64 depending on the system. This part was answered on the cython list: http://groups.google.com/group/cython-users/browse_thread/thread/f8022ee7ccbf7c5b Short answer is that it is not part of cython 0.14.1 but it is easy to add. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy type mismatch
On Fri, Jun 10, 2011 at 6:35 PM, Charles R Harris
charlesr.har...@gmail.com wrote:
On Fri, Jun 10, 2011 at 5:19 PM, Olivier Delalleau sh...@keba.be wrote:
But isn't it a bug if numpy.dtype('i') != numpy.dtype('l') on a 32 bit
computer where both are int32?
Maybe yes, maybe no ;) They have different descriptors, so from numpy's
perspective they are different, but at the hardware/precision level they are
the same. It's more of a decision as to what != means in this case. Since
numpy started as Numeric with only the c types the current behavior is
consistent, but that doesn't mean it shouldn't change at some point.
Maybe this is the same question, but are you maybe yes, maybe no on this too:
type(np.sum([1, 2, 3], dtype=np.int32)) == np.int32
False
Ben, what happens if you put an axis in there? Like
np.sum([[1, 2, 3], [4,5,6]], axis=0).dtype == np.int32
Just wondering if this is another different-dtype-for-different-axis case.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] [job] Python Job at Hedge Fund
We are looking for help to predict tomorrow's stock returns. The challenge is model selection in the presence of noisy data. The tools are ubuntu, python, cython, c, numpy, scipy, la, bottleneck, git. A quantitative background and experience or interest in model selection, machine learning, and software development are a plus. This is a full time position in Berkeley, California, two blocks from UC Berkeley. If you are interested send a CV or similar (or questions) to '.'.join(['htiek','scitylanayelekreb@namdoog','moc'][::-1])[::-1] ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Returning the same dtype in Cython as np.argmax
What is the rule to determine the dtype returned by numpy functions that return indices such as np.argmax? I assumed that np.argmax() returned the same dtype as np.int_. That works on linux32/64 and win32 but on win-amd64 np.int_ is int32 and np.argmax() returns int64. Someone suggested using intp. But I can't figure out how to do that in Cython since I am unable to cimport NPY_INTP. So this doesn't work: from numpy cimport NPY_INTP cdef np.ndarray[np.intp_t, ndim=1] output = PyArray_EMPTY(1, dims, NPY_INTP, 0) Is there another way? All I can think of is checking whether np.intp is np.int32 or int64 when I template the code. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] k maximal elements
On Sun, Jun 5, 2011 at 11:15 PM, Alex Ter-Sarkissov ater1...@gmail.com wrote: I have a vector of positive integers length n. Is there a simple (i.e. without sorting/ranking) of 'pulling out' k larrgest (or smallest) values. Something like sum(x[sum(x,1)(max(sum(x,1)+min(sum(x,1/2,]) but smarter You could use bottleneck [1]. It has a partial sort written in cython: a = np.arange(10) np.random.shuffle(a) a array([6, 3, 8, 5, 4, 2, 0, 1, 9, 7]) import bottleneck as bn k = 3 b = bn.partsort(a, a.size-k) b[-k:] array([8, 9, 7]) Or you could get the indices: idx = bn.argpartsort(a, a.size-k) idx[-k:] array([2, 8, 9]) [1] The partial sort will be in Bottleneck 0.5.0. There is a beta2 at https://github.com/kwgoodman/bottleneck/downloads. Release version will be posted at http://pypi.python.org/pypi/Bottleneck ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] k maximal elements
On Mon, Jun 6, 2011 at 6:44 AM, Keith Goodman kwgood...@gmail.com wrote: On Sun, Jun 5, 2011 at 11:15 PM, Alex Ter-Sarkissov ater1...@gmail.com wrote: I have a vector of positive integers length n. Is there a simple (i.e. without sorting/ranking) of 'pulling out' k larrgest (or smallest) values. Something like sum(x[sum(x,1)(max(sum(x,1)+min(sum(x,1/2,]) but smarter You could use bottleneck [1]. It has a partial sort written in cython: a = np.arange(10) np.random.shuffle(a) a array([6, 3, 8, 5, 4, 2, 0, 1, 9, 7]) import bottleneck as bn k = 3 b = bn.partsort(a, a.size-k) b[-k:] array([8, 9, 7]) Or you could get the indices: idx = bn.argpartsort(a, a.size-k) idx[-k:] array([2, 8, 9]) [1] The partial sort will be in Bottleneck 0.5.0. There is a beta2 at https://github.com/kwgoodman/bottleneck/downloads. Release version will be posted at http://pypi.python.org/pypi/Bottleneck Oh, and here are some timings: a = np.random.rand(1e6) timeit a.sort() 10 loops, best of 3: 18.3 ms per loop timeit bn.partsort(a, 100) 100 loops, best of 3: 3.71 ms per loop a = np.random.rand(1e4) timeit a.sort() 1000 loops, best of 3: 170 us per loop timeit bn.partsort(a, 10) 1 loops, best of 3: 19.4 us per loop ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] k maximal elements
On Mon, Jun 6, 2011 at 6:57 AM, gary ruben gru...@bigpond.net.au wrote: I learn a lot by watching the numpy and scipy lists (today Olivier taught me about heapq :), but he may not have noticed that Python 2.4 added an nsmallest method) import heapq q = list(x) heapq.heapify(q) k_smallest = heapq.nsmallest(k,q) I learned something new too---heapq. Timings: alist = a.tolist() heapq.heapify(alist) timeit k_smallest = heapq.nsmallest(k, alist) 100 loops, best of 3: 1.84 ms per loop timeit bn.partsort(a, a.size-10) 1 loops, best of 3: 39.3 us per loop ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] New functions.
On Tue, May 31, 2011 at 8:41 PM, Charles R Harris charlesr.har...@gmail.com wrote: On Tue, May 31, 2011 at 8:50 PM, Bruce Southey bsout...@gmail.com wrote: How about including all or some of Keith's Bottleneck package? He has tried to include some of the discussed functions and tried to make them very fast. I don't think they are sufficiently general as they are limited to 2 dimensions. However, I think the moving filters should go into scipy, either in ndimage or maybe signals. Some of the others we can still speed of significantly, for instance nanmedian, by using the new functionality in numpy, i.e., numpy sort has worked with nans for a while now. It looks like call overhead dominates the nanmax times for small arrays and this might improve if the ufunc machinery is cleaned up a bit more, I don't know how far Mark got with that. Currently Bottleneck accelerates 1d, 2d, and 3d input. Anything else falls back to a slower, non-cython version of the function. The same goes for int32, int64, float32, float64. It should not be difficult to extend to higher nd and more dtypes since everything is generated from template. The problem is that there would be a LOT of cython auto-generated C code since there is a separate function for each ndim, dtype, axis combination. Each of the ndim, dtype, axis functions currently has its own copy of the algorithm (such as median). Pulling that out and reusing it should save a lot of trees by reducing the auto-generated C code size. I recently added a partsort and argpartsort. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] random seed replicate 2d randn with 1d loop
On Mon, May 23, 2011 at 11:33 AM, josef.p...@gmail.com wrote:
I have a function in two versions, one vectorized, one with loop
the vectorized function gets all randn variables in one big array
rvs = distr.rvs(args, **{'size':(nobs, nrep)})
the looping version has:
for irep in xrange(nrep):
rvs = distr.rvs(args, **{'size':nobs})
the rest should be identical (except for vectorization
Is there a guarantee that the 2d arrays are filled up in a specific
order so that the loop and vectorized version produce the same result,
given the same seed?
Are you pulling the numbers from rows or columns of the 2d array?
Columns seem to work:
rs = np.random.RandomState([1,2,3])
rs.randn(3,3)
array([[ 0.89858245, 0.25528877, 0.95172625],
[-0.05663392, 0.54721555, 0.11512385],
[ 0.82495129, 0.17252144, 0.74570118]])
which gives the same as
rs = np.random.RandomState([1,2,3])
rs.randn(3)
array([ 0.89858245, 0.25528877, 0.95172625])
rs.randn(3)
array([-0.05663392, 0.54721555, 0.11512385])
rs.randn(3)
array([ 0.82495129, 0.17252144, 0.74570118])
I get similar results with np.random.seed
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] random seed replicate 2d randn with 1d loop
On Mon, May 23, 2011 at 11:42 AM, Keith Goodman kwgood...@gmail.com wrote:
On Mon, May 23, 2011 at 11:33 AM, josef.p...@gmail.com wrote:
I have a function in two versions, one vectorized, one with loop
the vectorized function gets all randn variables in one big array
rvs = distr.rvs(args, **{'size':(nobs, nrep)})
the looping version has:
for irep in xrange(nrep):
rvs = distr.rvs(args, **{'size':nobs})
the rest should be identical (except for vectorization
Is there a guarantee that the 2d arrays are filled up in a specific
order so that the loop and vectorized version produce the same result,
given the same seed?
Are you pulling the numbers from rows or columns of the 2d array?
Columns seem to work:
Sorry, I meant rows.
rs = np.random.RandomState([1,2,3])
rs.randn(3,3)
array([[ 0.89858245, 0.25528877, 0.95172625],
[-0.05663392, 0.54721555, 0.11512385],
[ 0.82495129, 0.17252144, 0.74570118]])
which gives the same as
rs = np.random.RandomState([1,2,3])
rs.randn(3)
array([ 0.89858245, 0.25528877, 0.95172625])
rs.randn(3)
array([-0.05663392, 0.54721555, 0.11512385])
rs.randn(3)
array([ 0.82495129, 0.17252144, 0.74570118])
I get similar results with np.random.seed
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] random seed replicate 2d randn with 1d loop
On Mon, May 23, 2011 at 12:34 PM, josef.p...@gmail.com wrote: Obviously I was working by columns, using a transpose worked, but rewriting to axis=1 instead of axis=0 which should be more efficient since I had almost all calculations by columns, I needed params = map(lambda x: np.expand_dims(x, 1), params) to get around broadcast errors that work automatically with axis=0 Here's another way to work by columns (note the input parameter oder='F'): rs = np.random.RandomState([1,2,3]) rs.randn(9).reshape(3,3,order='F') array([[ 0.89858245, -0.05663392, 0.82495129], [ 0.25528877, 0.54721555, 0.17252144], [ 0.95172625, 0.11512385, 0.74570118]]) rs = np.random.RandomState([1,2,3]) rs.randn(3) array([ 0.89858245, 0.25528877, 0.95172625]) rs.randn(3) array([-0.05663392, 0.54721555, 0.11512385]) rs.randn(3) array([ 0.82495129, 0.17252144, 0.74570118]) ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Mapping of dtype to C types
On Mon, May 9, 2011 at 1:46 AM, Pauli Virtanen p...@iki.fi wrote: Sun, 08 May 2011 14:45:45 -0700, Keith Goodman wrote: I'm writing a function that accepts four possible dtypes: int32, int64, float32, float64. The function will call a C extension (wrapped in Cython). What are the equivalent C types? int, long, float, double, respectively? Will that work on all systems? Long can be 32-bit or 64-bit, depending on the platform. The types available in Numpy are listed here, including the information which of them are compatible with which C types: http://docs.scipy.org/doc/numpy/reference/arrays.scalars.html The C long seems not to be listed -- but it's the same as Python int, i.e., np.int_ will work. IIRC, the Numpy type codes, dtype.kind, map directly to C types. Does this mapping look right? np.float32 -- float np.float64 -- double if np.int_ == np.int32: # np.int32 -- long # np.int64 -- long long elif np.int_ == np.int64: # np.int32 -- int # np.int64 -- long ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Mapping of dtype to C types
I'm writing a function that accepts four possible dtypes: int32, int64, float32, float64. The function will call a C extension (wrapped in Cython). What are the equivalent C types? int, long, float, double, respectively? Will that work on all systems? ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] ANN: Numpy 1.6.0 release candidate 2
On Tue, May 3, 2011 at 11:18 AM, Ralf Gommers
ralf.gomm...@googlemail.com wrote:
I am pleased to announce the availability of the second release
candidate of NumPy 1.6.0.
I get one failure when I run from ipython but not python. In ipython I
import a few packages at startup. One of those packages must be
changing the numpy error settings (would that explain it?). Should
that be a failure?
The failing test:
class TestSeterr(TestCase):
def test_default(self):
err = geterr()
self.assertEqual(err, dict(
divide='warn',
invalid='warn',
over='warn',
under='ignore',
))
==
FAIL: test_default (test_numeric.TestSeterr)
--
Traceback (most recent call last):
File
/usr/local/lib/python2.6/dist-packages/numpy/core/tests/test_numeric.py,
line 232, in test_default
under='ignore',
AssertionError: {'over': 'warn', 'divide': 'ignore', 'invalid':
'ignore', 'under': 'ignore'} != {'over': 'warn', 'divide': 'warn',
'invalid': 'warn', 'under': 'ignore'}
--
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] ANN: Numpy 1.6.0 release candidate 2
On Tue, May 3, 2011 at 11:18 AM, Ralf Gommers ralf.gomm...@googlemail.com wrote: I am pleased to announce the availability of the second release candidate of NumPy 1.6.0. nanmin and nanmax are much faster in Numpy 1.6. Plus they now return an object that has dtype, etc attributes when the input is all NaN. Broke Bottleneck unit tests and makes the benchmarks less impressive but: Nice. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Array as Variable using from cdms2 import MV2 as MV
On Mon, Apr 25, 2011 at 12:17 AM, josef.p...@gmail.com wrote: On Mon, Apr 25, 2011 at 2:50 AM, dileep kunjaai dileepkunj...@gmail.com wrote: Dear sir, I am have 2 mxn numpy array say obs fcst. I have to calculate sum of squre of (obs[i, j]-fcst[i, j]) using from cdms2 import MV2 as MV in CDAT without using for loop. For example: obs= [0.6 1.1 0.02 0.2 0.2 0.8 0. 0. 0.4 0.8 0.5 5.5 1.5 0.5 1.5 3.5 0.5 1.5 5.0 2.6 5.1 4.1 3.2 2.3 1.5 4.4 0.9 1.5 2. 2.3 1.1 1.1 1.5 12.6 1.3 2.2 12 1.7 1.6 15 1.9 1.5 0.9 2.5 5.5 ] fcst= [0.7 0.1 0.2 0.2 0.2 0.3 0.8 0. 0. 0. 0.5 0.5 0.5 0.5 0.5 0.7 1. 1.5 2. 2.6 5.1 4.1 3.2 2.3 1.5 0.7 1. 1.5 2. 2.3 1.1 1.1 1.1 12.7 1.3 2.2 2. 1.7 1.6 1.5 1.9 1.5 0.9 0.5 7.5] here obs and fcst are numpy array I give obs=MV.array(obs) fcst=MV.array(fcst) Then it become sbst=obs-fcst subst= [[ -0.1 1. -0.18 0. 0. ] [ 0.5 -0.8 0. 0.4 0.8 ] [ 0. 5. 1. 0. 1. ] [ 2.8 -0.5 0. 3. 0. ] [ 0. 0. 0. 0. 0. ] [ 3.7 -0.1 0. 0. 0. ] [ 0. 0. 0.4 -0.1 0. ] [ 0. 10. 0. 0. 13.5 ] [ 0. 0. 0. 2. -2. ]] But i dont know how to find sum of squre of each term(Actually my aim is to finding MEAN SQUARED ERROR) (sbst**2).sum() or with sum along columns (sbst**2).sum(0) explanation is in the documentation If speed is an issue then the development version of bottleneck (0.5.dev) has a fast sum of squares function: a = np.random.rand(1000, 10) timeit (a * a).sum(1) 1 loops, best of 3: 143 us per loop timeit bn.ss(a, 1) 10 loops, best of 3: 14.3 us per loop It is a faster replacement for scipy.stats.ss: from scipy import stats timeit stats.ss(a, 1) 1 loops, best of 3: 159 us per loop Speed up factor depends on the shape of the input array and the axis: timeit stats.ss(a, 0) 1 loops, best of 3: 42.8 us per loop timeit bn.ss(a, 0) 10 loops, best of 3: 17.3 us per loop ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] argmin and argmax without nan
On Thu, Mar 24, 2011 at 6:19 AM, Ralf Gommers ralf.gomm...@googlemail.com wrote: 2011/3/24 Dmitrey tm...@ukr.net: hi, is there any way to get argmin and argmax of an array w/o nans? Currently I have from numpy import * argmax([10,nan,100]) 1 argmin([10,nan,100]) 1 But it's not the values I would like to get. The walkaround I use: get all indeces of nans, replace them by -inf, get argmax, replace them by inf, get argmin. Is there any better way? (BTW, I invoke argmin/argmax along of a chosen axis) D. In [3]: np.nanargmax([10, np.nan, 100]) Out[3]: 2 In [4]: np.nanargmin([10, np.nan, 100]) Out[4]: 0 And if speed is an issue (it usually isn't) you can use the nanargmax from Bottleneck: a = np.random.rand(1) a[a 0.5] = np.nan timeit np.nanargmax(a) 1 loops, best of 3: 127 us per loop import bottleneck as bn timeit bn.nanargmax(a) 10 loops, best of 3: 12.4 us per loop ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] moving window product
On Mon, Mar 21, 2011 at 10:10 AM, Brent Pedersen bpede...@gmail.com wrote: hi, is there a way to take the product along a 1-d array in a moving window? -- similar to convolve, with product in place of sum? currently, i'm column_stacking the array with offsets of itself into window_size columns and then taking the product at axis 1. like:: w = np.column_stack(a[i:-window_size+i] for i in range(0, window_size)) window_product = np.product(w, axis=1) but then there are the edge effects/array size issues--like those handled in np.convolve. it there something in numpy/scipy that addresses this. or that does the column_stacking with an offset? The Bottleneck package has a fast moving window sum (bn.move_sum and bn.move_nansum). You could use that along with a = np.random.rand(5) a.prod() 0.015877866878931741 np.exp(np.log(a).sum()) 0.015877866878931751 Or you could use strides or scipy.ndimage as in https://github.com/kwgoodman/bottleneck/blob/master/bottleneck/slow/move.py ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] moving window product
On Mon, Mar 21, 2011 at 10:34 AM, Brent Pedersen bpede...@gmail.com wrote: On Mon, Mar 21, 2011 at 11:19 AM, Keith Goodman kwgood...@gmail.com wrote: On Mon, Mar 21, 2011 at 10:10 AM, Brent Pedersen bpede...@gmail.com wrote: hi, is there a way to take the product along a 1-d array in a moving window? -- similar to convolve, with product in place of sum? currently, i'm column_stacking the array with offsets of itself into window_size columns and then taking the product at axis 1. like:: w = np.column_stack(a[i:-window_size+i] for i in range(0, window_size)) window_product = np.product(w, axis=1) but then there are the edge effects/array size issues--like those handled in np.convolve. it there something in numpy/scipy that addresses this. or that does the column_stacking with an offset? The Bottleneck package has a fast moving window sum (bn.move_sum and bn.move_nansum). You could use that along with a = np.random.rand(5) a.prod() 0.015877866878931741 np.exp(np.log(a).sum()) 0.015877866878931751 Or you could use strides or scipy.ndimage as in https://github.com/kwgoodman/bottleneck/blob/master/bottleneck/slow/move.py ah yes, of course. thank you. def moving_product(a, window_size, mode=same): return np.exp(np.convolve(np.log(a), np.ones(window_size), mode)) i'll have a closer look at your strided version in bottleneck as well. I don't know what size problem you are working on or if speed is an issue, but here are some timings: a = np.random.rand(100) window_size = 1000 timeit np.exp(np.convolve(np.log(a), np.ones(window_size), 'same')) 1 loops, best of 3: 889 ms per loop timeit np.exp(bn.move_sum(np.log(a), window_size)) 10 loops, best of 3: 82.5 ms per loop Most all that time is spent in np.exp(np.log(a)): timeit bn.move_sum(a, window_size) 100 loops, best of 3: 3.72 ms per loop So I assume if I made a bn.move_prod the time would be around 200x compared to convolve. BTW, you could do the exp inplace: timeit b = bn.move_sum(np.log(a), window_size); np.exp(b, b) 10 loops, best of 3: 76.3 ms per loop ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] moving window product
On Mon, Mar 21, 2011 at 11:27 AM, Brent Pedersen bpede...@gmail.com wrote:
my current use-case is to do this 24 times on arrays of about 200K elements.
file IO is the major bottleneck.
Would using h5py or pytables help? I get about 3 ms for a write-read
cycle for a 200K array. That's much faster than the convolve speed.
import h5py
import time
f = h5py.File('speed.h5py')
a = np.random.rand(20)
t1 = time.time(); f['a'] = a; b = f['a'][:]; t2 = time.time()
print %0.2f ms % (1000*(t2 - t1))
3.30 ms
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] assert_almost_equal bug?
On Sat, Mar 12, 2011 at 4:16 AM, Ralf Gommers ralf.gomm...@googlemail.com wrote: On Sat, Mar 12, 2011 at 12:10 AM, Keith Goodman kwgood...@gmail.com wrote: assert_almost_equal() and assert_array_almost_equal() raise a ValueError instead of an AssertionError when the array contains np.inf: That's a bug, is fixed in 45269ee1. assert_array_compare was checking for nans but not for infs. Thanks, Ralf! ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] assert_almost_equal bug?
assert_almost_equal() and assert_array_almost_equal() raise a ValueError instead of an AssertionError when the array contains np.inf: a = np.array([[1., 2.], [3., 4.]]) b = a.copy() np.testing.assert_almost_equal(a, b) b[0,0] = np.inf np.testing.assert_almost_equal(a, b) snip ValueError: Arrays are not almost equal x: array([[ 1., 2.], [ 3., 4.]]) y: array([[ inf, 2.], [ 3., 4.]]) np.testing.assert_array_almost_equal(a, b) snip ValueError: Arrays are not almost equal x: array([[ 1., 2.], [ 3., 4.]]) y: array([[ inf, 2.], [ 3., 4.]]) ticket: http://projects.scipy.org/numpy/ticket/1769 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] How to get the prices of Moving Averages Crosses?
On Tue, Mar 1, 2011 at 8:07 AM, Andre Lopes lopes80an...@gmail.com wrote: Hi, I'm new to Numpy. I'm doing some tests with some Stock Market Quotes My struggle right now is how to get the values of the moving averages crosses, I send an image in attach to illustrate what I'm trying to get. I'm using the this computation to get when the moving averages crosses, but when I look at the graph, the values doesn't seem ok. [quote] # Get when the ma20 cross ma50 equal = np.round(ma20,2)==np.round(ma50,2) dates_cross = (dates[equal]) prices_cross = (prices[equal]) [/quote] The full code is this: [quote] # Modules import datetime import numpy as np import matplotlib.finance as finance import matplotlib.mlab as mlab import matplotlib.pyplot as plot # Define quote startdate = datetime.date(2008,10,1) today = enddate = datetime.date.today() ticker = 'uso' # Catch CSV fh = finance.fetch_historical_yahoo(ticker, startdate, enddate) # From CSV to REACARRAY r = mlab.csv2rec(fh); fh.close() # Order by Desc r.sort() ### Methods Begin def moving_average(x, n, type='simple'): compute an n period moving average. type is 'simple' | 'exponential' x = np.asarray(x) if type=='simple': weights = np.ones(n) else: weights = np.exp(np.linspace(-1., 0., n)) weights /= weights.sum() a = np.convolve(x, weights, mode='full')[:len(x)] a[:n] = a[n] return a ### Methods End prices = r.adj_close dates = r.date ma20 = moving_average(prices, 20, type='simple') ma50 = moving_average(prices, 50, type='simple') # Get when the ma20 cross ma50 equal = np.round(ma20,2)==np.round(ma50,2) dates_cross = (dates[equal]) prices_cross = (prices[equal]) # Ver se a ma20 ma50 # ma20_greater_than_ma50 = np.round(ma20,2) np.round(ma50,2) # dates_ma20_greater_than_ma50 = (dates[ma20_greater_than_ma50]) # prices_ma20_greater_than_ma50 = (prices[ma20_greater_than_ma50]) print dates_cross print prices_cross #print dates_ma20_greater_than_ma50 #print prices_ma20_greater_than_ma50 plot.plot(prices) plot.plot(ma20) plot.plot(ma50) plot.show() [/quote] Someone can give me some clues? I gave it a try using a labeled array: from la.data.yahoo import quotes data = quotes(['uso'], date1=(2008,10,1)) close = data.lix[:,['close']] ma20 = close.mov_mean(20) ma50 = close.mov_mean(50) lo = ma20 = ma50 hi = ma20 = ma50 cross_up = lo.lag(1) hi cross_dn = hi.lag(1) lo Dates when ma20 goes above ma50: cross_up[cross_up.x] label_0 2009-03-24 2009-08-14 2009-10-20 2010-01-15 2010-03-09 2010-07-16 2010-10-08 x array([ True, True, True, True, True, True, True], dtype=bool) Dates when it goes below: cross_dn[cross_dn.x] label_0 2009-07-17 2009-09-23 2009-12-09 2010-02-04 2010-05-10 2010-08-27 2011-02-04 x array([ True, True, True, True, True, True, True], dtype=bool) I didn't double check my logic, but hopefully it will give you some ideas to play with. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] When memory access is a bottleneck
A topic that often comes up on the list is that arr.sum(axis=-1) is faster than arr.sum(axis=0). For C ordered arrays, moving along the last axis moves the smallest amount in memory. And moving small amounts in memory keeps the data in cache longer. Can I use that fact to speed up calculations along axis=0? As a test I used the nanmean function from the Bottleneck package. I created a second version of the function that calculates the nanmean of two columns at a time (code below). I figured that would be a more efficient use of cache. It is faster for many array sizes but not for all. Is there anything I can do about that? Is the speed up really due to cache? The smallest arrays I tried should fit entirely in cache, so perhaps unrolling the loop is helping the compiler be more efficient? Yeah, so obviously, I don't understand what is going on. Faster: a = np.random.rand(10,10) timeit nanmean_2d_float64_axis0(a) 100 loops, best of 3: 1.24 us per loop timeit nanmean_2d_float64_axis0_double(a) 100 loops, best of 3: 1.17 us per loop a = np.random.rand(100,100) timeit nanmean_2d_float64_axis0(a) 10 loops, best of 3: 16 us per loop timeit nanmean_2d_float64_axis0_double(a) 10 loops, best of 3: 15.7 us per loop a = np.random.rand(1000,1000) timeit nanmean_2d_float64_axis0(a) 100 loops, best of 3: 8.57 ms per loop timeit nanmean_2d_float64_axis0_double(a) 100 loops, best of 3: 5.02 ms per loop a = np.random.rand(1,100) timeit nanmean_2d_float64_axis0(a) 100 loops, best of 3: 3.52 ms per loop timeit nanmean_2d_float64_axis0_double(a) 100 loops, best of 3: 3.35 ms per loop Slower: a = np.random.rand(100,1) timeit nanmean_2d_float64_axis0(a) 100 loops, best of 3: 2.22 ms per loop timeit nanmean_2d_float64_axis0_double(a) 100 loops, best of 3: 2.57 ms per loop Code (Bottleneck simplified BSD license): @cython.boundscheck(False) @cython.wraparound(False) def nanmean_2d_float64_axis0(np.ndarray[np.float64_t, ndim=2] a): Mean of 2d array with dtype=float64 along axis=0 ignoring NaNs. cdef int count = 0 cdef np.float64_t asum = 0, ai cdef Py_ssize_t i0, i1 cdef np.npy_intp *dim dim = PyArray_DIMS(a) cdef int n0 = dim[0] cdef int n1 = dim[1] cdef np.npy_intp *dims = [n1] cdef np.ndarray[np.float64_t, ndim=1] y = PyArray_EMPTY(1, dims, NPY_float64, 0) for i1 in range(n1): asum = 0 count = 0 for i0 in range(n0): ai = a[i0, i1] if ai == ai: asum += ai count += 1 if count 0: y[i1] = asum / count else: y[i1] = NAN return y @cython.boundscheck(False) @cython.wraparound(False) def nanmean_2d_float64_axis0_double(np.ndarray[np.float64_t, ndim=2] a): Mean of 2d array with dtype=float64 along axis=0 ignoring NaNs. cdef int count = 0, count2 = 0 cdef np.float64_t asum = 0, asum2 = 0, ai cdef Py_ssize_t i0, i1, i11 cdef np.npy_intp *dim dim = PyArray_DIMS(a) cdef int n0 = dim[0] cdef int n1 = dim[1] cdef np.npy_intp *dims = [n1] cdef np.ndarray[np.float64_t, ndim=1] y = PyArray_EMPTY(1, dims, NPY_float64, 0) for i1 in range(0,n1,2): asum = 0 count = 0 asum2 = 0 count2 = 0 i11 = i1 + 1 for i0 in range(n0): ai = a[i0, i1] if ai == ai: asum += ai count += 1 ai = a[i0, i11] if ai == ai: asum2 += ai count2 += 1 if count 0: y[i1] = asum / count else: y[i1] = NAN if count2 0: y[i11] = asum / count else: y[i11] = NAN return y ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] NaN value processing in weave.inline code
On Fri, Jan 14, 2011 at 12:03 PM, Joon Ro joonp...@gmail.com wrote:
Hi,
I was wondering if it is possible to process (in if statement - check if the
given value is NaN) numpy NaN value inside the weave.inline c code.
testcode = '''
if (test(0)) {
return_val = test(0);
}
'''
err = weave.inline(testcode,
['test'],
type_converters = converters.blitz, force = 0, verbose = 1)
with test(0) = nan returns err = nan correctly, but I don't know how to
check the nan value inside the c inline c code. Is there any way I can get
similar functionality as isnan?
To check if a scalar, x, is NaN:
if x == x:
# No, it is not a NaN
else:
# Yes, it is a NaN
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Output dtype
On Wed, Jan 12, 2011 at 8:20 AM, Bruce Southey bsout...@gmail.com wrote:
On 12/13/2010 04:53 PM, Keith Goodman wrote:
On Mon, Dec 13, 2010 at 12:20 PM, Bruce Southeybsout...@gmail.com wrote:
Unless something has changed since the docstring was written, this is
probably an inherited 'bug' from np.mean() as the author expected that
the docstring of mean was correct. For my 'old' 2.0 dev version:
np.mean( np.array([[0,1,2,3,4,5]], dtype='float32'), axis=1).dtype
dtype('float32')
np.mean( np.array([[0,1,2,3,4,5]], dtype='float32')).dtype
dtype('float64')
Are you saying the bug is in the doc string, the output, or both? I
think it is both; I expect the second result above to be float32.
Sorry as I filed a bug for this as 1710
http://projects.scipy.org/numpy/ticket/1710
but this is the same as ticket 518 that is listed as won't fix:
http://projects.scipy.org/numpy/ticket/518
I fixed ticket 518 in bottleneck:
a = np.array([1,2,3], dtype='float32')
bn.median(a).dtype
dtype('float32')
np.median(a).dtype
dtype('float64')
Not sure I would have done that if I knew that numpy has a won't fix on it.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Rolling window (moving average, moving std, and more)
On Tue, Jan 4, 2011 at 8:14 AM, Keith Goodman kwgood...@gmail.com wrote: On Tue, Jan 4, 2011 at 8:06 AM, Sebastian Haase seb.ha...@gmail.com wrote: On Mon, Jan 3, 2011 at 5:32 PM, Erik Rigtorp e...@rigtorp.com wrote: On Mon, Jan 3, 2011 at 11:26, Eric Firing efir...@hawaii.edu wrote: Instead of calculating statistics independently each time the window is advanced one data point, the statistics are updated. I have not done any benchmarking, but I expect this approach to be quick. This might accumulate numerical errors. But could be fine for many applications. The code is old; I have not tried to update it to take advantage of cython's advances over pyrex. If I were writing it now, I might not bother with the C level at all; it could all be done in cython, probably with no speed penalty, and maybe even with reduced overhead. No doubt this would be faster, I just wanted to offer a general way to this in NumPy. ___ BTW, some of these operations can be done using scipy's ndimage - right ? Any comments ? How does the performance compare ? ndimage might have more options regarding edge handling, or ? Take a look at the moving window function in the development version of the la package: https://github.com/kwgoodman/la/blob/master/la/farray/mov.py Many of the moving window functions offer three calculation methods: filter (ndimage), strides (the strides trick discussed in this thread), and loop (a simple python loop). For example: a = np.random.rand(500,2000) timeit la.farray.mov_max(a, window=252, axis=-1, method='filter') 1 loops, best of 3: 336 ms per loop timeit la.farray.mov_max(a, window=252, axis=-1, method='strides') 1 loops, best of 3: 609 ms per loop timeit la.farray.mov_max(a, window=252, axis=-1, method='loop') 1 loops, best of 3: 638 ms per loop No one method is best for all situations. That is one of the reasons I started the Bottleneck package. I figured Cython could beat them all. I added four new function to Bottleneck: move_min, move_max, move_nanmin, move_nanmax. They are much faster than using SciPy's ndimage.maximum_filter1d or the strides trick: a = np.random.rand(500,2000) timeit la.farray.mov_max(a, window=252, axis=-1, method='filter') # ndimage 1 loops, best of 3: 336 ms per loop timeit bn.move_max(a, window=252, axis=-1) # bottleneck 100 loops, best of 3: 14.1 ms per loop That looks too good to be true. Are the outputs the same? a1 = la.farray.mov_max(a, window=252, axis=-1, method='filter') a2 = bn.move_max(a, window=252, axis=-1) np.testing.assert_array_almost_equal(a1, a2) Yes. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] aa.astype(int) truncates and doesn't round
On Thu, Jan 6, 2011 at 2:14 AM, josef.p...@gmail.com wrote: just something I bumped into and wasn't aware of aa array([ 1., 1., 1., 1., 1.]) aa.astype(int) array([0, 1, 0, 0, 0]) aa - 1 array([ -2.22044605e-16, 2.22044605e-16, -2.22044605e-16, -3.33066907e-16, -3.33066907e-16]) np.round(aa).astype(int) array([1, 1, 1, 1, 1]) a = np.ones(100) a.astype(int) array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) My default numpy int is 64 bits. Try 32 bits: a = np.ones(100, np.int32) a.astype(int) array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Rolling window (moving average, moving std, and more)
On Tue, Jan 4, 2011 at 8:06 AM, Sebastian Haase seb.ha...@gmail.com wrote: On Mon, Jan 3, 2011 at 5:32 PM, Erik Rigtorp e...@rigtorp.com wrote: On Mon, Jan 3, 2011 at 11:26, Eric Firing efir...@hawaii.edu wrote: Instead of calculating statistics independently each time the window is advanced one data point, the statistics are updated. I have not done any benchmarking, but I expect this approach to be quick. This might accumulate numerical errors. But could be fine for many applications. The code is old; I have not tried to update it to take advantage of cython's advances over pyrex. If I were writing it now, I might not bother with the C level at all; it could all be done in cython, probably with no speed penalty, and maybe even with reduced overhead. No doubt this would be faster, I just wanted to offer a general way to this in NumPy. ___ BTW, some of these operations can be done using scipy's ndimage - right ? Any comments ? How does the performance compare ? ndimage might have more options regarding edge handling, or ? Take a look at the moving window function in the development version of the la package: https://github.com/kwgoodman/la/blob/master/la/farray/mov.py Many of the moving window functions offer three calculation methods: filter (ndimage), strides (the strides trick discussed in this thread), and loop (a simple python loop). For example: a = np.random.rand(500,2000) timeit la.farray.mov_max(a, window=252, axis=-1, method='filter') 1 loops, best of 3: 336 ms per loop timeit la.farray.mov_max(a, window=252, axis=-1, method='strides') 1 loops, best of 3: 609 ms per loop timeit la.farray.mov_max(a, window=252, axis=-1, method='loop') 1 loops, best of 3: 638 ms per loop No one method is best for all situations. That is one of the reasons I started the Bottleneck package. I figured Cython could beat them all. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Rolling window (moving average, moving std, and more)
On Fri, Dec 31, 2010 at 8:29 PM, Erik Rigtorp e...@rigtorp.com wrote: Implementing moving average, moving std and other functions working over rolling windows using python for loops are slow. This is a effective stride trick I learned from Keith Goodman's kwgood...@gmail.com Bottleneck code but generalized into arrays of any dimension. This trick allows the loop to be performed in C code and in the future hopefully using multiple cores. I like using strides for moving window functions. The one downside I found is that it is slow when window * (arr.shape[axis] - window) is large: a = np.random.rand(100) b = rolling_window(a, 5000) import bottleneck as bn timeit bn.nanmean(b, axis=1) 1 loops, best of 3: 7.1 s per loop timeit bn.move_nanmean(a, window=5000, axis=0) 100 loops, best of 3: 7.99 ms per loop ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Rolling window (moving average, moving std, and more)
On Mon, Jan 3, 2011 at 5:37 AM, Erik Rigtorp e...@rigtorp.com wrote: It's only a view of the array, no copying is done. Though some operations like np.std() will copy the array, but that's more of a bug. In general It's hard to imagine any speedup gains by copying a 10GB array. I don't think that np.std makes a copy of the input data if the input is an array. If the input is, for example, a list, then an array is created. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Rolling window (moving average, moving std, and more)
On Mon, Jan 3, 2011 at 7:41 AM, Erik Rigtorp e...@rigtorp.com wrote: On Mon, Jan 3, 2011 at 10:36, Keith Goodman kwgood...@gmail.com wrote: On Mon, Jan 3, 2011 at 5:37 AM, Erik Rigtorp e...@rigtorp.com wrote: It's only a view of the array, no copying is done. Though some operations like np.std() will copy the array, but that's more of a bug. In general It's hard to imagine any speedup gains by copying a 10GB array. I don't think that np.std makes a copy of the input data if the input is an array. If the input is, for example, a list, then an array is created. When I tried it on a big array, it tried to allocate a huge amount of memory. As I said it's probably a bug. Yes, that would be a big bug. np.std does have to initialize the output array. If the window size is small compared to arr.shape[axis] then the memory taken by the output array is of the same order as that of the input array. Could that be what you are seeing? a = np.arange(10) Small window, output array shape (8,): rolling_window(a, 2) array([[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9]]) Big window, output array shape (2,): rolling_window(a, 9) array([[0, 1, 2, 3, 4, 5, 6, 7, 8], [1, 2, 3, 4, 5, 6, 7, 8, 9]]) ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Faster NaN functions
On Fri, Dec 31, 2010 at 8:21 AM, Lev Givon l...@columbia.edu wrote: Received from Erik Rigtorp on Fri, Dec 31, 2010 at 08:52:53AM EST: Hi, I just send a pull request for some faster NaN functions, https://github.com/rigtorp/numpy. I implemented the following generalized ufuncs: nansum(), nancumsum(), nanmean(), nanstd() and for fun mean() and std(). It turns out that the generalized ufunc mean() and std() is faster than the current numpy functions. I'm also going to add nanprod(), nancumprod(), nanmax(), nanmin(), nanargmax(), nanargmin(). The current implementation is not optimized in any way and there are probably some speedups possible. I hope we can get this into numpy 2.0, me and people around me seems to have a need for these functions. Erik ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion How does this compare to Bottleneck? http://pypi.python.org/pypi/Bottleneck/ I had all sorts of problems with ABI differences (this is the first time I've tried numpy 2.0). So I couldn't get ipython, etc to work with Erik's new nan functions. That's why my speed comparison below might be hard to follow and only tests one example. For timing I used bottleneck's autotimeit function: from bottleneck.benchmark.autotimeit import autotimeit First Erik's new nanmean: stmt = nanmean2(a.flat) setup = import numpy as np; from numpy.core.umath_tests import nanmean as nanmean2; rs=np.random.RandomState([1,2,3]); a = rs.rand(100,100) autotimeit(stmt, setup) 5.1356482505798338e-05 Bottleneck's low level nanmean: stmt = nanmean(a) setup = import numpy as np; from bottleneck.func import nanmean_2d_float64_axisNone as nanmean; rs=np.random.RandomState([1,2,3]); a = rs.rand(100,100) autotimeit(stmt, setup) 1.5422070026397704e-05 Bottleneck's high level nanmean: setup = import numpy as np; from bottleneck.func import nanmean; rs=np.random.RandomState([1,2,3]); a = rs.rand(100,100) autotimeit(stmt, setup) 1.7850480079650879e-05 Numpy's mean: setup = import numpy as np; from numpy import mean; rs=np.random.RandomState([1,2,3]); a = rs.rand(100,100) stmt = mean(a) autotimeit(stmt, setup) 1.6718170642852782e-05 Scipy's nanmean: setup = import numpy as np; from scipy.stats import nanmean; rs=np.random.RandomState([1,2,3]); a = rs.rand(100,100) stmt = nanmean(a) autotimeit(stmt, setup) 0.00024667191505432128 The tests above should be repeated for arrays that contain NaNs, and for different array sizes and different axes. Bottleneck's benchmark suite can be modified to do all that but I can't import Erik's new numpy and bottleneck at the same time at the moment. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] NumPy C-API equivalent of np.float64()
On Tue, Dec 28, 2010 at 11:22 PM, Robert Bradshaw
rober...@math.washington.edu wrote:
On Tue, Dec 28, 2010 at 8:10 PM, John Salvatier
jsalv...@u.washington.edu wrote:
Wouldn't that be a cast? You do casts in Cython with double(expression)
and that should be the equivalent of float64 I think.
Or even numpy.float64_t (expression) if you've cimported numpy
(though as mentioned this is the same as double on every platform I
know of). Even easier is just to use the expression in a the right
context and it will convert it for you.
That will give me a float object but it will not have dtype, shape,
ndim, etc methods.
m = np.mean([1,2,3])
m
2.0
m.dtype
dtype('float64')
m.ndim
0
using np.float64_t gives:
AttributeError: 'float' object has no attribute 'dtype'
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] NumPy C-API equivalent of np.float64()
On Wed, Dec 29, 2010 at 9:37 AM, Robert Bradshaw
rober...@math.washington.edu wrote:
On Wed, Dec 29, 2010 at 9:05 AM, Keith Goodman kwgood...@gmail.com wrote:
On Tue, Dec 28, 2010 at 11:22 PM, Robert Bradshaw
rober...@math.washington.edu wrote:
On Tue, Dec 28, 2010 at 8:10 PM, John Salvatier
jsalv...@u.washington.edu wrote:
Wouldn't that be a cast? You do casts in Cython with double(expression)
and that should be the equivalent of float64 I think.
Or even numpy.float64_t (expression) if you've cimported numpy
(though as mentioned this is the same as double on every platform I
know of). Even easier is just to use the expression in a the right
context and it will convert it for you.
That will give me a float object but it will not have dtype, shape,
ndim, etc methods.
m = np.mean([1,2,3])
m
2.0
m.dtype
dtype('float64')
m.ndim
0
using np.float64_t gives:
AttributeError: 'float' object has no attribute 'dtype'
Well, in this case I doubt your'e going to be able to do much better
than np.float64(expr), as the bulk or the time is probably spent in
object allocation (and you're really asking for an object here). If
you knew the right C calls, you might be able to get a 2x speedup.
A factor of 2 would be great! A tenth of a micro second is a lot of
overhead for small input arrays.
I'm guessing it is one of these functions but I don't understand the
signatures (nor ref counting):
PyObject* PyArray_Scalar(void* data, PyArray_Descr* dtype, PyObject* itemsize)
Return an array scalar object of the given enumerated typenum and
itemsize by copying from memory pointed to by data . If swap is
nonzero then this function will byteswap the data if appropriate to
the data-type because array scalars are always in correct machine-byte
order.
PyObject* PyArray_ToScalar(void* data, PyArrayObject* arr)
Return an array scalar object of the type and itemsize indicated
by the array object arr copied from the memory pointed to by data and
swapping if the data in arr is not in machine byte-order.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] NumPy C-API equivalent of np.float64()
On Wed, Dec 29, 2010 at 9:48 AM, Matthew Brett matthew.br...@gmail.com wrote:
Hi,
On Wed, Dec 29, 2010 at 5:37 PM, Robert Bradshaw
rober...@math.washington.edu wrote:
On Wed, Dec 29, 2010 at 9:05 AM, Keith Goodman kwgood...@gmail.com wrote:
On Tue, Dec 28, 2010 at 11:22 PM, Robert Bradshaw
rober...@math.washington.edu wrote:
On Tue, Dec 28, 2010 at 8:10 PM, John Salvatier
jsalv...@u.washington.edu wrote:
Wouldn't that be a cast? You do casts in Cython with double(expression)
and that should be the equivalent of float64 I think.
Or even numpy.float64_t (expression) if you've cimported numpy
(though as mentioned this is the same as double on every platform I
know of). Even easier is just to use the expression in a the right
context and it will convert it for you.
That will give me a float object but it will not have dtype, shape,
ndim, etc methods.
m = np.mean([1,2,3])
m
2.0
m.dtype
dtype('float64')
m.ndim
0
using np.float64_t gives:
AttributeError: 'float' object has no attribute 'dtype'
Forgive me if I haven't understood your question, but can you use
PyArray_DescrFromType with e.g NPY_FLOAT64 ?
I'm pretty hopeless here. I don't know how to put all that together in
a function.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] NumPy C-API equivalent of np.float64()
On Wed, Dec 29, 2010 at 10:13 AM, Matthew Brett matthew.br...@gmail.com wrote: Forgive me if I haven't understood your question, but can you use PyArray_DescrFromType with e.g NPY_FLOAT64 ? I'm pretty hopeless here. I don't know how to put all that together in a function. That might be because I'm not understanding you very well, but I was thinking that: cdef dtype descr = PyArray_DescrFromType(NPY_FLOAT64) would give you the float64 dtype that I thought you wanted? I'm shooting from the hip here, in between nieces competing for the computer and my attention. I think I need a function. One that does this: n = 10.0 hasattr(n, 'ndim') False m = np.float64(n) hasattr(m, 'ndim') True np.float64 is fast, just hoping someone had a C-API inline version of np.float64() that is faster. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] NumPy C-API equivalent of np.float64()
On Wed, Dec 29, 2010 at 11:43 AM, Matthew Brett matthew.br...@gmail.com wrote: Hi, That might be because I'm not understanding you very well, but I was thinking that: cdef dtype descr = PyArray_DescrFromType(NPY_FLOAT64) would give you the float64 dtype that I thought you wanted? I'm shooting from the hip here, in between nieces competing for the computer and my attention. I think I need a function. One that does this: n = 10.0 hasattr(n, 'ndim') False m = np.float64(n) hasattr(m, 'ndim') True Now the nieces have gone, I see that I did completely misunderstand. I think you want the C-API calls to be able to create a 0-dim ndarray object from a python float. There was a thread on C-API array creation on the cython list a little while ago: http://www.mail-archive.com/cython-dev@codespeak.net/msg07703.html Code in scipy here: https://github.com/scipy/scipy-svn/blob/master/scipy/io/matlab/mio5_utils.pyx See around line 36 there, and 432, and the header file I copied from Dag Sverre: https://github.com/scipy/scipy-svn/blob/master/scipy/io/matlab/numpy_rephrasing.h As you can see, it's a little horrible, in that you have to take care to get the references right to the dtype and to the data. I actually did not investigate in detail whether this lower-level array creation was speeding my code up much. I hope that's more useful... Wow! That's a mouthful of code. Yes, very handy to have an example to work from. Thank you. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] NumPy C-API equivalent of np.float64()
On Wed, Dec 29, 2010 at 11:54 AM, Pauli Virtanen p...@iki.fi wrote: Keith Goodman wrote: np.float64 is fast, just hoping someone had a C-API inline version of np.float64() that is faster. You're looking for PyArrayScalar_New and _ASSIGN. See https://github.com/numpy/numpy/blob/master/numpy/core/include/numpy/arrayscalars.h Undocumented (bad), but AFAIK public. Those look nice. I'm stuck since I can't cimport them. I'll have to read up on how to tell cython about those functions. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] NumPy C-API equivalent of np.float64()
I'm looking for the C-API equivalent of the np.float64 function, something that I could use inline in a Cython function. I don't know how to write the function. Anyone have one sitting around? I'd like to use it, if it is faster than np.float64 (np.int32, np.float32, ...) in the Bottleneck package when the output is a scalar, for example bn.median(arr, axis=None). ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] How construct custom slice
On Mon, Dec 27, 2010 at 10:36 AM, Mario Moura moura.ma...@gmail.com wrote: Hi Folks a = np.zeros((4,3,5,55,5),dtype='|S8') myLen = 4 # here I use myLen = len(something) li = [3,2,4] # li from a list.append(something) sl = slice(0,myLen) tmpIndex = tuple(li) + sl + 4 # == Here my problem a[tmpIndex] # So What I want is: fillMe = np.array(['foo','bar','hello','world']) # But I cant contruct by hand like this a[3,2,4,:4,4] = fillMe a Again. I need construct custom slice from here tmpIndex = tuple(li) + sl + 4 a[tmpIndex] First let's do it by hand: a = np.zeros((4,3,5,55,5),dtype='|S8') fillMe = np.array(['foo','bar','hello','world']) a[3,2,4,:4,4] = fillMe Now let's try using an index: b = np.zeros((4,3,5,55,5),dtype='|S8') myLen = 4 li = [3,2,4] sl = slice(0,myLen) Make index: idx = range(a.ndim) idx[:3] = li idx[3] = sl idx[4] = 4 idx = tuple(idx) Compare results: b[idx] = fillMe (a == b).all() True ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Output dtype
From the np.median doc string: If the input contains integers, or
floats of smaller precision than 64, then the output data-type is
float64.
arr = np.array([[0,1,2,3,4,5]], dtype='float32')
np.median(arr, axis=0).dtype
dtype('float32')
np.median(arr, axis=1).dtype
dtype('float32')
np.median(arr, axis=None).dtype
dtype('float64')
So the output doesn't agree with the doc string.
What is the desired dtype of the accumulator and the output for when
the input dtype is less than float64? Should it depend on axis?
I'm trying to duplicate the behavior of np.median (and other
numpy/scipy functions) in the Bottleneck package and am running into a
few corner cases while unit testing.
Here's another one:
np.sum([np.nan]).dtype
dtype('float64')
np.nansum([1,np.nan]).dtype
dtype('float64')
np.nansum([np.nan]).dtype
snip
AttributeError: 'float' object has no attribute 'dtype'
I just duplicated the numpy behavior for that one since it was easy to do.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Output dtype
On Mon, Dec 13, 2010 at 12:20 PM, Bruce Southey bsout...@gmail.com wrote:
On 12/13/2010 11:59 AM, Keith Goodman wrote:
From the np.median doc string: If the input contains integers, or
floats of smaller precision than 64, then the output data-type is
float64.
arr = np.array([[0,1,2,3,4,5]], dtype='float32')
np.median(arr, axis=0).dtype
dtype('float32')
np.median(arr, axis=1).dtype
dtype('float32')
np.median(arr, axis=None).dtype
dtype('float64')
So the output doesn't agree with the doc string.
What is the desired dtype of the accumulator and the output for when
the input dtype is less than float64? Should it depend on axis?
I'm trying to duplicate the behavior of np.median (and other
numpy/scipy functions) in the Bottleneck package and am running into a
few corner cases while unit testing.
Here's another one:
np.sum([np.nan]).dtype
dtype('float64')
np.nansum([1,np.nan]).dtype
dtype('float64')
np.nansum([np.nan]).dtype
snip
AttributeError: 'float' object has no attribute 'dtype'
I just duplicated the numpy behavior for that one since it was easy to do.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Unless something has changed since the docstring was written, this is
probably an inherited 'bug' from np.mean() as the author expected that
the docstring of mean was correct. For my 'old' 2.0 dev version:
np.mean( np.array([[0,1,2,3,4,5]], dtype='float32'), axis=1).dtype
dtype('float32')
np.mean( np.array([[0,1,2,3,4,5]], dtype='float32')).dtype
dtype('float64')
Same issue with np.std and np.var.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Output dtype
On Mon, Dec 13, 2010 at 12:20 PM, Bruce Southey bsout...@gmail.com wrote:
Unless something has changed since the docstring was written, this is
probably an inherited 'bug' from np.mean() as the author expected that
the docstring of mean was correct. For my 'old' 2.0 dev version:
np.mean( np.array([[0,1,2,3,4,5]], dtype='float32'), axis=1).dtype
dtype('float32')
np.mean( np.array([[0,1,2,3,4,5]], dtype='float32')).dtype
dtype('float64')
Are you saying the bug is in the doc string, the output, or both? I
think it is both; I expect the second result above to be float32.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] A Cython apply_along_axis function
On Wed, Dec 1, 2010 at 6:07 PM, Keith Goodman kwgood...@gmail.com wrote: On Wed, Dec 1, 2010 at 5:53 PM, David da...@silveregg.co.jp wrote: On 12/02/2010 04:47 AM, Keith Goodman wrote: It's hard to write Cython code that can handle all dtypes and arbitrary number of dimensions. The former is typically dealt with using templates, but what do people do about the latter? The only way that I know to do that systematically is iterator. There is a relatively simple example in scipy/signal (lfilter.c.src). I wonder if it would be possible to add better support for numpy iterators in cython... Thanks for the tip. I'm starting to think that for now I should just template both dtype and ndim. I ended up templating both dtype and axis. For the axis templating I used two functions: looper and loop_cdef. LOOPER Make a 3d loop template: loop = ''' for iINDEX0 in range(nINDEX0): for iINDEX1 in range(nINDEX1): amin = MAXDTYPE for iINDEX2 in range(nINDEX2): ai = a[INDEXALL] if ai = amin: amin = ai y[INDEXPOP] = amin ''' Import the looper function: from bottleneck.src.template.template import looper Make a loop over axis=0: print looper(loop, ndim=3, axis=0) for i1 in range(n1): for i2 in range(n2): amin = MAXDTYPE for i0 in range(n0): ai = a[i0, i1, i2] if ai = amin: amin = ai y[i1, i2] = amin Make a loop over axis=1: print looper(loop, ndim=3, axis=1) for i0 in range(n0): for i2 in range(n2): amin = MAXDTYPE for i1 in range(n1): ai = a[i0, i1, i2] if ai = amin: amin = ai y[i0, i2] = amin LOOP_CDEF Define parameters: ndim = 3 dtype = 'float64' axis = 1 is_reducing_function = True Import loop_cdef: from bottleneck.src.template.template import loop_cdef Make loop initialization code: print loop_cdef(ndim, dtype, axis, is_reducing_function) cdef Py_ssize_t i0, i1, i2 cdef int n0 = a.shape[0] cdef int n1 = a.shape[1] cdef int n2 = a.shape[2] cdef np.npy_intp *dims = [n0, n2] cdef np.ndarray[np.float64_t, ndim=2] y = PyArray_EMPTY(2, dims, NPY_float64, 0) ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] np.var() and ddof
Why does ddof=2 and ddof=3 give the same result? np.var([1, 2, 3], ddof=0) 0.3 np.var([1, 2, 3], ddof=1) 1.0 np.var([1, 2, 3], ddof=2) 2.0 np.var([1, 2, 3], ddof=3) 2.0 np.var([1, 2, 3], ddof=4) -2.0 I expected NaN for ddof=3. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] np.var() and ddof
On Fri, Dec 10, 2010 at 2:26 PM, josef.p...@gmail.com wrote: On Fri, Dec 10, 2010 at 4:42 PM, Keith Goodman kwgood...@gmail.com wrote: Why does ddof=2 and ddof=3 give the same result? np.var([1, 2, 3], ddof=0) 0.3 np.var([1, 2, 3], ddof=1) 1.0 np.var([1, 2, 3], ddof=2) 2.0 np.var([1, 2, 3], ddof=3) 2.0 np.var([1, 2, 3], ddof=4) -2.0 I expected NaN for ddof=3. It's a floating point calculation, so I would expect np.inf Right, NAFN (F=Finite). Unless, of course, the numerator is zero too. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] A Cython apply_along_axis function
It's hard to write Cython code that can handle all dtypes and arbitrary number of dimensions. The former is typically dealt with using templates, but what do people do about the latter? I'm trying to take baby steps towards writing an apply_along_axis function that takes as input a cython function, numpy array, and axis. I'm using the following numpy ticket as a guide but I'm really just copying and pasting without understanding: http://projects.scipy.org/numpy/attachment/ticket/1213/_selectmodule.pyx Can anyone spot why I get a segfault on the call to nanmean_1d in apply_along_axis? import numpy as np cimport numpy as np import cython cdef double NAN = double np.nan ctypedef np.float64_t (*func_t)(void *buf, np.npy_intp size, np.npy_intp s) def apply_along_axis(np.ndarray[np.float64_t, ndim=1] a, int axis): cdef func_t nanmean_1d cdef np.npy_intp stride, itemsize cdef int ndim = a.ndim cdef np.float64_t out itemsize = np.npy_intp a.itemsize if ndim == 1: stride = a.strides[0] // itemsize # convert stride bytes -- items out = nanmean_1d(a.data, a.shape[0], stride) else: raise ValueError(Not yet coded) return out cdef np.float64_t nanmean_1d(void *buf, np.npy_intp n, np.npy_intp s): nanmean of buffer. cdef np.float64_t *a = np.float64_t * buf # cdef np.npy_intp i, count = 0 cdef np.float64_t asum, ai if s == 1: for i in range(n): ai = a[i] if ai == ai: asum += ai count += 1 else: for i in range(n): ai = a[i*s] if ai == ai: asum += ai count += 1 if count 0: return asum / count else: return NAN ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] A Cython apply_along_axis function
On Wed, Dec 1, 2010 at 5:53 PM, David da...@silveregg.co.jp wrote: On 12/02/2010 04:47 AM, Keith Goodman wrote: It's hard to write Cython code that can handle all dtypes and arbitrary number of dimensions. The former is typically dealt with using templates, but what do people do about the latter? The only way that I know to do that systematically is iterator. There is a relatively simple example in scipy/signal (lfilter.c.src). I wonder if it would be possible to add better support for numpy iterators in cython... Thanks for the tip. I'm starting to think that for now I should just template both dtype and ndim. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Warning: invalid value encountered in subtract
After upgrading from numpy 1.4.1 to 1.5.1 I get warnings like Warning: invalid value encountered in subtract when I run unit tests (or timeit) using python -c 'blah' but not from an interactive session. How can I tell the warnings to go away? ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] A faster median (Wirth's method)
On Tue, Sep 1, 2009 at 2:37 PM, Sturla Molden stu...@molden.no wrote: Dag Sverre Seljebotn skrev: Nitpick: This will fail on large arrays. I guess numpy.npy_intp is the right type to use in this case? By the way, here is a more polished version, does it look ok? http://projects.scipy.org/numpy/attachment/ticket/1213/generate_qselect.py http://projects.scipy.org/numpy/attachment/ticket/1213/quickselect.pyx This is my favorite numpy/scipy ticket. So I am happy that I can contribute in a small way by pointing out a bug. The search for the k-th smallest element is only done over the first k elements (that's the bug) instead of over the entire array. Specifically while l k should be while l r. I added a median function to the Bottleneck package: https://github.com/kwgoodman/bottleneck Timings: import bottleneck as bn arr = np.random.rand(100, 100) timeit np.median(arr) 1000 loops, best of 3: 762 us per loop timeit bn.median(arr) 1 loops, best of 3: 198 us per loop What other functions could be built from a selection algorithm? nanmedian scoreatpercentile quantile knn select others? But before I add more functions to the package I need to figure out how to make a cython apply_along_axis function. For the first release I am hand coding the 1d, 2d, and 3d cases. Boring to write, hard to maintain, and doesn't solve the nd case. Does anyone have a cython apply_along_axis that takes a cython reducing function as input? The ticket has an example but I couldn't get it to run. If no one has one (the horror!) I'll begin to work on one sometime after the first release. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] A faster median (Wirth's method)
On Tue, Nov 30, 2010 at 11:25 AM, John Salvatier jsalv...@u.washington.edu wrote: I am very interested in this result. I have wanted to know how to do an My first thought was to write the reducing function like this cdef np.float64_t namean(np.ndarray[np.float64_t, ndim=1] a): but cython doesn't allow np.ndarray in a cdef. That's why the ticket (URL earlier in the thread) uses a (the?) buffer interface to the array. The particular example in the ticket is not a reducing function, it works on the data in place. We can change that. I can set up a sandbox in the Bottleneck project if anyone (John!) is interested in working on this. I plan to get a first (preview) release out soon and then take a break. The first project for the second release is this; second project is templating. Then the third release is just turning the crank and adding more functions. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] A faster median (Wirth's method)
On Tue, Nov 30, 2010 at 11:58 AM, Matthew Brett matthew.br...@gmail.com wrote: Hi, On Tue, Nov 30, 2010 at 11:35 AM, Keith Goodman kwgood...@gmail.com wrote: On Tue, Nov 30, 2010 at 11:25 AM, John Salvatier jsalv...@u.washington.edu wrote: I am very interested in this result. I have wanted to know how to do an My first thought was to write the reducing function like this cdef np.float64_t namean(np.ndarray[np.float64_t, ndim=1] a): but cython doesn't allow np.ndarray in a cdef. Sorry for the ill-considered hasty reply, but do you mean that this: import numpy as np cimport numpy as cnp cdef cnp.float64_t namean(cnp.ndarray[cnp.float64_t, ndim=1] a): return np.nanmean(a) # just a placeholder is not allowed? It works for me. Is it a cython version thing? (I've got 0.13), Oh, that's nice! I'm using 0.11.2. OK, time to upgrade. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Warning: invalid value encountered in subtract
On Tue, Nov 30, 2010 at 2:25 PM, Robert Kern robert.k...@gmail.com wrote:
On Tue, Nov 30, 2010 at 16:22, Keith Goodman kwgood...@gmail.com wrote:
On Tue, Nov 30, 2010 at 1:41 PM, Skipper Seabold jsseab...@gmail.com wrote:
On Tue, Nov 30, 2010 at 1:34 PM, Keith Goodman kwgood...@gmail.com wrote:
After upgrading from numpy 1.4.1 to 1.5.1 I get warnings like
Warning: invalid value encountered in subtract when I run unit tests
(or timeit) using python -c 'blah' but not from an interactive
session. How can I tell the warnings to go away?
If it's this type of floating point related stuff, you can use np.seterr
In [1]: import numpy as np
In [2]: np.log(1./np.array(0))
Warning: divide by zero encountered in divide
Out[2]: inf
In [3]: orig_settings = np.seterr()
In [4]: np.seterr(all=ignore)
Out[4]: {'divide': 'print', 'invalid': 'print', 'over': 'print', 'under':
'ignor
e'}
In [5]: np.log(1./np.array(0))
Out[5]: inf
In [6]: np.seterr(**orig_settings)
Out[6]: {'divide': 'ignore', 'invalid': 'ignore', 'over': 'ignore',
'under': 'ig
nore'}
In [7]: np.log(1./np.array(0))
Warning: divide by zero encountered in divide
Out[7]: inf
I have been using the orig_settings so that I can take over the
control of this from the user and then set it back to how it was.
Thank, Skipper. That works. Do you wrap it in a try...except? And then
raise whatever brought you to the exception? Sounds like a pain.
Is it considered OK for a package to change the state of np.seterr if
there is an error? Silly question. I'm just looking for an easy fix.
with np.errstate(invalid='ignore'):
Ah! Thank you, Robert!
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Does np.std() make two passes through the data?
On Sun, Nov 21, 2010 at 5:56 PM, Robert Kern robert.k...@gmail.com wrote: On Sun, Nov 21, 2010 at 19:49, Keith Goodman kwgood...@gmail.com wrote: But this sample gives a difference: a = np.random.rand(100) a.var() 0.080232196646619805 var(a) 0.080232196646619791 As you know, I'm trying to make a drop-in replacement for scipy.stats.nanstd. Maybe I'll have to add an asterisk to the drop-in part. Either that, or suck it up and store the damn mean. The difference is less than eps. Quite possibly, the one-pass version is even closer to the true value than the two-pass version. I wrote 3 cython prototype implementations of nanstd for 1d float64 arrays: a = np.random.rand(100) # numpy; doesn't take care of NaNs a.std() 0.28852169850186793 # cython of np.sqrt(((arr - arr.mean())**2).mean()) nanstd_twopass(a, ddof=0) 0.28852169850186798 # cython of np.sqrt((arr*arr).mean() - arr.mean()**2) nanstd_simple(a, ddof=0) 0.28852169850187437 # http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#On-line_algorithm nanstd_online(a, ddof=0) 0.28852169850187243 # My target, scipy version from scipy.stats import nanstd nanstd(a, bias=True) 0.28852169850186798 Timing: timeit nanstd(a, bias=True) 10 loops, best of 3: 27.8 ms per loop timeit a.std() 100 loops, best of 3: 11.5 ms per loop timeit nanstd_twopass(a, ddof=0) 100 loops, best of 3: 3.24 ms per loop timeit nanstd_simple(a, ddof=0) 1000 loops, best of 3: 1.6 ms per loop timeit nanstd_online(a, ddof=0) 100 loops, best of 3: 10.8 ms per loop nanstd_simple is the fastest but I assume the algorithm is no good for general use? I think I'll go with nanstd_twopass. It will most closely match numpy/scipy, is more robust than nanstd_simple, and is the second fastest. Here's the code. Improvements welcomed. @cython.boundscheck(False) @cython.wraparound(False) def nanstd_simple(np.ndarray[np.float64_t, ndim=1] a, int ddof): nanstd of 1d numpy array with dtype=np.float64 along axis=0. cdef Py_ssize_t i cdef int a0 = a.shape[0], count = 0 cdef np.float64_t asum = 0, a2sum=0, ai for i in range(a0): ai = a[i] if ai == ai: asum += ai a2sum += ai * ai count += 1 if count 0: asum = asum * asum return sqrt((a2sum - asum / count) / (count - ddof)) else: return np.float64(NAN) @cython.boundscheck(False) @cython.wraparound(False) def nanstd_online(np.ndarray[np.float64_t, ndim=1] a, int ddof): nanstd of 1d numpy array with dtype=np.float64 along axis=0. cdef Py_ssize_t i cdef int a0 = a.shape[0], n = 0 cdef np.float64_t mean = 0, M2 = 0, delta, x for i in range(a0): x = a[i] if x == x: n += 1 delta = x - mean mean = mean + delta / n M2 = M2 + delta * (x - mean) if n 0: return np.float64(sqrt(M2 / (n - ddof))) else: return np.float64(NAN) @cython.boundscheck(False) @cython.wraparound(False) def nanstd_twopass(np.ndarray[np.float64_t, ndim=1] a, int ddof): nanstd of 1d numpy array with dtype=np.float64 along axis=0. cdef Py_ssize_t i cdef int a0 = a.shape[0], count = 0 cdef np.float64_t asum = 0, a2sum=0, amean, ai, da for i in range(a0): ai = a[i] if ai == ai: asum += ai count += 1 if count 0: amean = asum / count asum = 0 for i in range(a0): ai = a[i] if ai == ai: da = ai - amean asum += da a2sum += (da * da) asum = asum * asum return sqrt((a2sum - asum / count) / (count - ddof)) else: return np.float64(NAN) ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Does np.std() make two passes through the data?
On Mon, Nov 22, 2010 at 9:13 AM, josef.p...@gmail.com wrote: Two pass would provide precision that we would expect in numpy, but I don't know if anyone ever tested the NIST problems for basic statistics. Here are the results for their most difficult dataset. But I guess running one test doesn't mean anything. http://www.itl.nist.gov/div898/strd/univ/addinfo/numacc4.html np.absolute(a.std(ddof=1) - 0.1) 5.5884095961911129e-10 np.absolute(nanstd_online(a, ddof=1) - 0.1) 5.5890501948763216e-10 np.absolute(nanstd_simple(a, ddof=1) - 0.1) nan # Ha! np.absolute(nanstd_twopass(a, ddof=1) - 0.1) 5.5879308125117433e-10 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Does np.std() make two passes through the data?
On Mon, Nov 22, 2010 at 10:51 AM, josef.p...@gmail.com wrote: On Mon, Nov 22, 2010 at 1:39 PM, Keith Goodman kwgood...@gmail.com wrote: On Mon, Nov 22, 2010 at 10:32 AM, josef.p...@gmail.com wrote: On Mon, Nov 22, 2010 at 1:26 PM, Keith Goodman kwgood...@gmail.com wrote: On Mon, Nov 22, 2010 at 9:03 AM, Keith Goodman kwgood...@gmail.com wrote: @cython.boundscheck(False) @cython.wraparound(False) def nanstd_twopass(np.ndarray[np.float64_t, ndim=1] a, int ddof): nanstd of 1d numpy array with dtype=np.float64 along axis=0. cdef Py_ssize_t i cdef int a0 = a.shape[0], count = 0 cdef np.float64_t asum = 0, a2sum=0, amean, ai, da for i in range(a0): ai = a[i] if ai == ai: asum += ai count += 1 if count 0: amean = asum / count asum = 0 for i in range(a0): ai = a[i] if ai == ai: da = ai - amean asum += da a2sum += (da * da) asum = asum * asum return sqrt((a2sum - asum / count) / (count - ddof)) else: return np.float64(NAN) This is 5% faster: @cython.boundscheck(False) @cython.wraparound(False) def nanstd_1d_float64_axis0_2(np.ndarray[np.float64_t, ndim=1] a, int ddof): nanstd of 1d numpy array with dtype=np.float64 along axis=0. cdef Py_ssize_t i cdef int a0 = a.shape[0], count = 0 cdef np.float64_t asum = 0, amean, ai for i in range(a0): ai = a[i] if ai == ai: asum += ai count += 1 if count 0: amean = asum / count asum = 0 for i in range(a0): ai = a[i] if ai == ai: ai -= amean asum += (ai * ai) return sqrt(asum / (count - ddof)) else: return np.float64(NAN) I think it would be better to write nanvar instead of nanstd and take the square root only in a delegating nanstd, instead of the other way around. (Also a change that should be made in scipy.stats) Yeah, I noticed that numpy does that. I was planning to have separate var and std functions. Here's why (from the readme file, but maybe I should template it, the sqrt automatically converts large ddof to NaN): I'm not sure what you are saying, dropping the squareroot in the function doesn't require nan handling in the inner loop. If you want to return nan when count-ddof=0, then you could just replace if count 0: ... by if count -ddof 0: ... Or am I missing the point? Yes, sorry. Ignore the sqrt/nan comment. The point is that I want to be able to return the underlying, low-level cython function (see below). So I need separate var and std versions to do that (unless I can modify default kwargs on the fly). Under the hood Nanny uses a separate Cython function for each combination of ndim, dtype, and axis. A lot of the overhead in ny.nanmax, for example, is in checking that your axis is within range, converting non-array data to an array, and selecting the function to use to calculate nanmax. You can get rid of the overhead by doing all this before you, say, enter an inner loop: arr = np.random.rand(10,10) axis = 0 func, a = ny.func.nanmax_selector(arr, axis) func.__name__ 'nanmax_2d_float64_axis0' Let's see how much faster than runs: timeit np.nanmax(arr, axis=0) 1 loops, best of 3: 25.7 us per loop timeit ny.nanmax(arr, axis=0) 10 loops, best of 3: 5.25 us per loop timeit func(a) 10 loops, best of 3: 2.5 us per loop Note that func is faster than the Numpy's non-nan version of max: timeit arr.max(axis=0) 10 loops, best of 3: 3.28 us per loop So adding NaN protection to your inner loops has a negative cost! ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Does np.std() make two passes through the data?
On Mon, Nov 22, 2010 at 11:00 AM, josef.p...@gmail.com wrote: I don't think that works for complex numbers. (statsmodels has now a preference that calculations work also for complex numbers) I'm only supporting int32, int64, float64 for now. Getting the other ints and floats should be easy. I don't have plans for complex numbers. If it's not a big change then someone can add support later. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] [ANN] Nanny, faster NaN functions
On Sun, Nov 21, 2010 at 10:25 AM, Wes McKinney wesmck...@gmail.com wrote: On Sat, Nov 20, 2010 at 7:24 PM, Keith Goodman kwgood...@gmail.com wrote: On Sat, Nov 20, 2010 at 3:54 PM, Wes McKinney wesmck...@gmail.com wrote: Keith (and others), What would you think about creating a library of mostly Cython-based domain specific functions? So stuff like rolling statistical moments, nan* functions like you have here, and all that-- NumPy-array only functions that don't necessarily belong in NumPy or SciPy (but could be included on down the road). You were already talking about this on the statsmodels mailing list for larry. I spent a lot of time writing a bunch of these for pandas over the last couple of years, and I would have relatively few qualms about moving these outside of pandas and introducing a dependency. You could do the same for larry-- then we'd all be relying on the same well-vetted and tested codebase. I've started working on moving window statistics cython functions. I plan to make it into a package called Roly (for rolling). The signatures are: mov_sum(arr, window, axis=-1) and mov_nansum(arr, window, axis=-1), etc. I think of Nanny and Roly as two separate packages. A narrow focus is good for a new package. But maybe each package could be a subpackage in a super package? Would the function signatures in Nanny (exact duplicates of the corresponding functions in Numpy and Scipy) work for pandas? I plan to use Nanny in larry. I'll try to get the structure of the Nanny package in place. But if it doesn't attract any interest after that then I may fold it into larry. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion Why make multiple packages? It seems like all these functions are somewhat related: practical tools for real-world data analysis (where observations are often missing). I suspect having everything under one hood would create more interest than chopping things up-- would be very useful to folks in many different disciplines (finance, economics, statistics, etc.). In R, for example, NA-handling is just a part of every day life. Of course in R there is a special NA value which is distinct from NaN-- many folks object to the use of NaN for missing values. The alternative is masked arrays, but in my case I wasn't willing to sacrifice so much performance for purity's sake. I could certainly use the nan* functions to replace code in pandas where I've handled things in a somewhat ad hoc way. A package focused on NaN-aware functions sounds like a good idea. I think a good plan would be to start by making faster, drop-in replacements for the NaN functions that are already in numpy and scipy. That is already a lot of work. After that, one possibility is to add stuff like nancumsum, nanprod, etc. After that moving window stuff? ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] [ANN] Nanny, faster NaN functions
On Sun, Nov 21, 2010 at 12:30 PM, josef.p...@gmail.com wrote: On Sun, Nov 21, 2010 at 2:48 PM, Keith Goodman kwgood...@gmail.com wrote: On Sun, Nov 21, 2010 at 10:25 AM, Wes McKinney wesmck...@gmail.com wrote: On Sat, Nov 20, 2010 at 7:24 PM, Keith Goodman kwgood...@gmail.com wrote: On Sat, Nov 20, 2010 at 3:54 PM, Wes McKinney wesmck...@gmail.com wrote: Keith (and others), What would you think about creating a library of mostly Cython-based domain specific functions? So stuff like rolling statistical moments, nan* functions like you have here, and all that-- NumPy-array only functions that don't necessarily belong in NumPy or SciPy (but could be included on down the road). You were already talking about this on the statsmodels mailing list for larry. I spent a lot of time writing a bunch of these for pandas over the last couple of years, and I would have relatively few qualms about moving these outside of pandas and introducing a dependency. You could do the same for larry-- then we'd all be relying on the same well-vetted and tested codebase. I've started working on moving window statistics cython functions. I plan to make it into a package called Roly (for rolling). The signatures are: mov_sum(arr, window, axis=-1) and mov_nansum(arr, window, axis=-1), etc. I think of Nanny and Roly as two separate packages. A narrow focus is good for a new package. But maybe each package could be a subpackage in a super package? Would the function signatures in Nanny (exact duplicates of the corresponding functions in Numpy and Scipy) work for pandas? I plan to use Nanny in larry. I'll try to get the structure of the Nanny package in place. But if it doesn't attract any interest after that then I may fold it into larry. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion Why make multiple packages? It seems like all these functions are somewhat related: practical tools for real-world data analysis (where observations are often missing). I suspect having everything under one hood would create more interest than chopping things up-- would be very useful to folks in many different disciplines (finance, economics, statistics, etc.). In R, for example, NA-handling is just a part of every day life. Of course in R there is a special NA value which is distinct from NaN-- many folks object to the use of NaN for missing values. The alternative is masked arrays, but in my case I wasn't willing to sacrifice so much performance for purity's sake. I could certainly use the nan* functions to replace code in pandas where I've handled things in a somewhat ad hoc way. A package focused on NaN-aware functions sounds like a good idea. I think a good plan would be to start by making faster, drop-in replacements for the NaN functions that are already in numpy and scipy. That is already a lot of work. After that, one possibility is to add stuff like nancumsum, nanprod, etc. After that moving window stuff? and maybe group functions after that? Yes, group functions are on my list. If there is a lot of repetition, you could use templating. Even simple string substitution, if it is only replacing the dtype, works pretty well. It would at least reduce some copy-paste. Unit test coverage should be good enough to mess around with trying templating. What's a good way to go? Write my own script that creates the .pyx file and call it from the make file? Or are there packages for doing the templating? I added nanmean (the first scipy function to enter nanny) and nanmin. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Does np.std() make two passes through the data?
Does np.std() make two passes through the data? Numpy: arr = np.random.rand(10) arr.std() 0.3008736260967052 Looks like an algorithm that makes one pass through the data (one for loop) wouldn't match arr.std(): np.sqrt((arr*arr).mean() - arr.mean()**2) 0.30087362609670526 But a slower two-pass algorithm would match arr.std(): np.sqrt(((arr - arr.mean())**2).mean()) 0.3008736260967052 Is there a way to get the same result as arr.std() in one pass (cython for loop) of the data? ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Does np.std() make two passes through the data?
On Sun, Nov 21, 2010 at 4:18 PM, josef.p...@gmail.com wrote: On Sun, Nov 21, 2010 at 6:43 PM, Keith Goodman kwgood...@gmail.com wrote: Does np.std() make two passes through the data? Numpy: arr = np.random.rand(10) arr.std() 0.3008736260967052 Looks like an algorithm that makes one pass through the data (one for loop) wouldn't match arr.std(): np.sqrt((arr*arr).mean() - arr.mean()**2) 0.30087362609670526 But a slower two-pass algorithm would match arr.std(): np.sqrt(((arr - arr.mean())**2).mean()) 0.3008736260967052 Is there a way to get the same result as arr.std() in one pass (cython for loop) of the data? reference several times pointed to on the list is the wikipedia page, e.g. http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#On-line_algorithm Unfortunately it doesn't give the same result as numpy's two pass algorithm. From wikipedia: def var(arr): n = 0 mean = 0 M2 = 0 for x in arr: n = n + 1 delta = x - mean mean = mean + delta / n M2 = M2 + delta * (x - mean) return M2 / n This set of random samples gives matching variance: a = np.random.rand(100) a.var() 0.07867478939716277 var(a) 0.07867478939716277 But this sample gives a difference: a = np.random.rand(100) a.var() 0.080232196646619805 var(a) 0.080232196646619791 As you know, I'm trying to make a drop-in replacement for scipy.stats.nanstd. Maybe I'll have to add an asterisk to the drop-in part. Either that, or suck it up and store the damn mean. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Does np.std() make two passes through the data?
On Sun, Nov 21, 2010 at 5:56 PM, Robert Kern robert.k...@gmail.com wrote: On Sun, Nov 21, 2010 at 19:49, Keith Goodman kwgood...@gmail.com wrote: But this sample gives a difference: a = np.random.rand(100) a.var() 0.080232196646619805 var(a) 0.080232196646619791 As you know, I'm trying to make a drop-in replacement for scipy.stats.nanstd. Maybe I'll have to add an asterisk to the drop-in part. Either that, or suck it up and store the damn mean. The difference is less than eps. Quite possibly, the one-pass version is even closer to the true value than the two-pass version. Good, it passes the Kern test. Here's an even more robust estimate: var(a - a.mean()) 0.080232196646619819 Which is better, numpy's two pass or the one pass on-line method? test() NumPy error: 9.31135e-18 Nanny error: 6.5745e-18 -- One pass wins! def test(n=10): numpy = 0 nanny = 0 for i in range(n): a = np.random.rand(10) truth = var(a - a.mean()) numpy += np.absolute(truth - a.var()) nanny += np.absolute(truth - var(a)) print 'NumPy error: %g' % (numpy / n) print 'Nanny error: %g' % (nanny / n) print ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] [ANN] Nanny, faster NaN functions
On Fri, Nov 19, 2010 at 7:42 PM, Keith Goodman kwgood...@gmail.com wrote: I should make a benchmark suite. ny.benchit(verbose=False) Nanny performance benchmark Nanny 0.0.1dev Numpy 1.4.1 Speed is numpy time divided by nanny time NaN means all NaNs Speed TestShapedtypeNaN? 6.6770 nansum(a, axis=-1) (500,500)int64 4.6612 nansum(a, axis=-1) (1,) float64 9.0351 nansum(a, axis=-1) (500,500)int32 3.0746 nansum(a, axis=-1) (500,500)float64 11.5740 nansum(a, axis=-1) (1,) int32 6.4484 nansum(a, axis=-1) (1,) int64 51.3917 nansum(a, axis=-1) (500,500)float64 NaN 13.8692 nansum(a, axis=-1) (1,) float64 NaN 6.5327 nanmax(a, axis=-1) (500,500)int64 8.8222 nanmax(a, axis=-1) (1,) float64 0.2059 nanmax(a, axis=-1) (500,500)int32 6.9262 nanmax(a, axis=-1) (500,500)float64 5.0688 nanmax(a, axis=-1) (1,) int32 6.5605 nanmax(a, axis=-1) (1,) int64 48.4850 nanmax(a, axis=-1) (500,500)float64 NaN 14.6289 nanmax(a, axis=-1) (1,) float64 NaN You can also use the makefile to run the benchmark: make bench ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] [ANN] Nanny, faster NaN functions
On Sat, Nov 20, 2010 at 3:54 PM, Wes McKinney wesmck...@gmail.com wrote: Keith (and others), What would you think about creating a library of mostly Cython-based domain specific functions? So stuff like rolling statistical moments, nan* functions like you have here, and all that-- NumPy-array only functions that don't necessarily belong in NumPy or SciPy (but could be included on down the road). You were already talking about this on the statsmodels mailing list for larry. I spent a lot of time writing a bunch of these for pandas over the last couple of years, and I would have relatively few qualms about moving these outside of pandas and introducing a dependency. You could do the same for larry-- then we'd all be relying on the same well-vetted and tested codebase. I've started working on moving window statistics cython functions. I plan to make it into a package called Roly (for rolling). The signatures are: mov_sum(arr, window, axis=-1) and mov_nansum(arr, window, axis=-1), etc. I think of Nanny and Roly as two separate packages. A narrow focus is good for a new package. But maybe each package could be a subpackage in a super package? Would the function signatures in Nanny (exact duplicates of the corresponding functions in Numpy and Scipy) work for pandas? I plan to use Nanny in larry. I'll try to get the structure of the Nanny package in place. But if it doesn't attract any interest after that then I may fold it into larry. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] [ANN] Nanny, faster NaN functions
= Nanny = Nanny uses the magic of Cython to give you a faster, drop-in replacement for the NaN functions in NumPy and SciPy. For example:: import nanny as ny import numpy as np arr = np.random.rand(100, 100) timeit np.nansum(arr) 1 loops, best of 3: 67.5 us per loop timeit ny.nansum(arr) 10 loops, best of 3: 18.2 us per loop Let's not forget to add some NaNs:: arr[arr 0.5] = np.nan timeit np.nansum(arr) 1000 loops, best of 3: 411 us per loop timeit ny.nansum(arr) 1 loops, best of 3: 65 us per loop Nanny uses a separate Cython function for each combination of ndim, dtype, and axis. You can get rid of a lot of overhead (useful in an inner loop, e.g.) by directly importing the function that matches your problem:: arr = np.random.rand(10, 10) from nansum import nansum_2d_float64_axis1 timeit np.nansum(arr, axis=1) 1 loops, best of 3: 25.5 us per loop timeit ny.nansum(arr, axis=1) 10 loops, best of 3: 5.15 us per loop timeit nansum_2d_float64_axis1(arr) 100 loops, best of 3: 1.75 us per loop I put together Nanny as a way to learn Cython. It currently only supports: - functions: nansum - Operating systems: 64-bit (accumulator for int32 is hard coded to int64) - dtype: int32, int64, float64 - ndim: 1, 2, and 3 If there is interest in the project, I could continue adding the remaining NaN functions from NumPy and SciPy: nanmin, nanmax, nanmean, nanmedian (using a partial sort), nanstd. But why stop there? How about nancumsum or nanprod? Or anynan, which could short-circuit once a NaN is found? Feedback on the code or the direction of the project are welcomed. So is coding help---without that I doubt the package will ever be completed. Once nansum is complete, many of the remaining functions will be copy, paste, touch up operations. Remember, Nanny quickly protects your precious data from the corrupting influence of Mr. Nan. License === Nanny is distributed under a Simplified BSD license. Parts of NumPy, which has a BSD licenses, are included in Nanny. See the LICENSE file, which is distributed with Nanny, for details. Installation You can grab Nanny at http://github.com/kwgoodman/nanny. nansum of ints is only supported by 64-bit operating systems at the moment. **GNU/Linux, Mac OS X, et al.** To install Nanny:: $ python setup.py build $ sudo python setup.py install Or, if you wish to specify where Nanny is installed, for example inside ``/usr/local``:: $ python setup.py build $ sudo python setup.py install --prefix=/usr/local **Windows** In order to compile the C code in Nanny you need a Windows version of the gcc compiler. MinGW (Minimalist GNU for Windows) contains gcc and has been used to successfully compile Nanny on Windows. Install MinGW and add it to your system path. Then install Nanny with the commands:: python setup.py build --compiler=mingw32 python setup.py install **Post install** After you have installed Nanny, run the suite of unit tests:: import nanny nanny.test() snip Ran 1 tests in 0.008s OK nose.result.TextTestResult run=1 errors=0 failures=0 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] [ANN] Nanny, faster NaN functions
On Fri, Nov 19, 2010 at 10:55 AM, Nathaniel Smith n...@pobox.com wrote: On Fri, Nov 19, 2010 at 10:33 AM, Keith Goodman kwgood...@gmail.com wrote: Nanny uses the magic of Cython to give you a faster, drop-in replacement for the NaN functions in NumPy and SciPy. Neat! Why not make this a patch to numpy/scipy instead? My guess is that having separate underlying functions for each dtype, ndim, and axis would be a nightmare for a large project like Numpy. But manageable for a focused project like nanny. Nanny uses a separate Cython function for each combination of ndim, dtype, and axis. You can get rid of a lot of overhead (useful in an inner loop, e.g.) by directly importing the function that matches your problem:: arr = np.random.rand(10, 10) from nansum import nansum_2d_float64_axis1 If this is really useful, then better to provide a function that finds the correct function for you? best_nansum = ny.get_best_nansum(ary[0, :, :], axis=1) for i in xrange(ary.shape[0]): best_nansum(ary[i, :, :], axis=1) That would be useful. It is what nanny.nansum does but it returns the sum instead of the function. - functions: nansum - Operating systems: 64-bit (accumulator for int32 is hard coded to int64) - dtype: int32, int64, float64 - ndim: 1, 2, and 3 What does it even mean to do NaN operations on integers? (I'd sometimes find it *really convenient* if there were a NaN value for standard computer integers... but there isn't?) Well, sometimes you write functions without knowing the dtype of the input. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] [ANN] Nanny, faster NaN functions
On Fri, Nov 19, 2010 at 12:10 PM, josef.p...@gmail.com wrote: What's the speed advantage of nanny compared to np.nansum that you have if the arrays are larger, say (1000,10) or (1,100) axis=0 ? Good point. In the small examples I showed so far maybe the speed up was all in overhead. Fortunately, that's not the case: arr = np.random.rand(1000, 1000) timeit np.nansum(arr) 100 loops, best of 3: 4.79 ms per loop timeit ny.nansum(arr) 1000 loops, best of 3: 1.53 ms per loop arr[arr 0.5] = np.nan timeit np.nansum(arr) 10 loops, best of 3: 44.5 ms per loop timeit ny.nansum(arr) 100 loops, best of 3: 6.18 ms per loop timeit np.nansum(arr, axis=0) 10 loops, best of 3: 52.3 ms per loop timeit ny.nansum(arr, axis=0) 100 loops, best of 3: 12.2 ms per loop np.nansum makes a copy of the input array and makes a mask (another copy) and then uses the mask to set the NaNs to zero in the copy. So not only is nanny faster, but it uses less memory. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] [ANN] Nanny, faster NaN functions
On Fri, Nov 19, 2010 at 12:19 PM, Pauli Virtanen p...@iki.fi wrote: Fri, 19 Nov 2010 11:19:57 -0800, Keith Goodman wrote: [clip] My guess is that having separate underlying functions for each dtype, ndim, and axis would be a nightmare for a large project like Numpy. But manageable for a focused project like nanny. Might be easier to migrate the nan* functions to using Ufuncs. Unless I'm missing something, np.nanmax - np.fmax.reduce np.nanmin - np.fmin.reduce For `nansum`, we'd need to add an ufunc `nanadd`, and for `nanargmax/min`, we'd need `argfmin/fmax'. How about that! I wasn't aware of fmax/fmin. Yes, I'd like a nanadd, please. arr = np.random.rand(1000, 1000) arr[arr 0.5] = np.nan np.nanmax(arr) 0.4625409581072 np.fmax.reduce(arr, axis=None) snip TypeError: an integer is required np.fmax.reduce(np.fmax.reduce(arr, axis=0), axis=0) 0.4625409581072 timeit np.fmax.reduce(np.fmax.reduce(arr, axis=0), axis=0) 100 loops, best of 3: 12.7 ms per loop timeit np.nanmax(arr) 10 loops, best of 3: 39.6 ms per loop timeit np.nanmax(arr, axis=0) 10 loops, best of 3: 46.5 ms per loop timeit np.fmax.reduce(arr, axis=0) 100 loops, best of 3: 12.7 ms per loop ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] [ANN] Nanny, faster NaN functions
On Fri, Nov 19, 2010 at 12:29 PM, Keith Goodman kwgood...@gmail.com wrote: On Fri, Nov 19, 2010 at 12:19 PM, Pauli Virtanen p...@iki.fi wrote: Fri, 19 Nov 2010 11:19:57 -0800, Keith Goodman wrote: [clip] My guess is that having separate underlying functions for each dtype, ndim, and axis would be a nightmare for a large project like Numpy. But manageable for a focused project like nanny. Might be easier to migrate the nan* functions to using Ufuncs. Unless I'm missing something, np.nanmax - np.fmax.reduce np.nanmin - np.fmin.reduce For `nansum`, we'd need to add an ufunc `nanadd`, and for `nanargmax/min`, we'd need `argfmin/fmax'. How about that! I wasn't aware of fmax/fmin. Yes, I'd like a nanadd, please. arr = np.random.rand(1000, 1000) arr[arr 0.5] = np.nan np.nanmax(arr) 0.4625409581072 np.fmax.reduce(arr, axis=None) snip TypeError: an integer is required np.fmax.reduce(np.fmax.reduce(arr, axis=0), axis=0) 0.4625409581072 timeit np.fmax.reduce(np.fmax.reduce(arr, axis=0), axis=0) 100 loops, best of 3: 12.7 ms per loop timeit np.nanmax(arr) 10 loops, best of 3: 39.6 ms per loop timeit np.nanmax(arr, axis=0) 10 loops, best of 3: 46.5 ms per loop timeit np.fmax.reduce(arr, axis=0) 100 loops, best of 3: 12.7 ms per loop Cython is faster than np.fmax.reduce. I wrote a cython version of np.nanmax, called nanmax below. (It only handles the 2d, float64, axis=None case, but since the array is large I don't think that explains the time difference). Note that fmax.reduce is slower than np.nanmax when there are no NaNs: arr = np.random.rand(1000, 1000) timeit np.nanmax(arr) 100 loops, best of 3: 5.82 ms per loop timeit np.fmax.reduce(np.fmax.reduce(arr)) 100 loops, best of 3: 9.14 ms per loop timeit nanmax(arr) 1000 loops, best of 3: 1.17 ms per loop arr[arr 0.5] = np.nan timeit np.nanmax(arr) 10 loops, best of 3: 45.5 ms per loop timeit np.fmax.reduce(np.fmax.reduce(arr)) 100 loops, best of 3: 12.7 ms per loop timeit nanmax(arr) 1000 loops, best of 3: 1.17 ms per loop ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] [ANN] Nanny, faster NaN functions
On Fri, Nov 19, 2010 at 7:19 PM, Charles R Harris charlesr.har...@gmail.com wrote: On Fri, Nov 19, 2010 at 1:50 PM, Keith Goodman kwgood...@gmail.com wrote: On Fri, Nov 19, 2010 at 12:29 PM, Keith Goodman kwgood...@gmail.com wrote: On Fri, Nov 19, 2010 at 12:19 PM, Pauli Virtanen p...@iki.fi wrote: Fri, 19 Nov 2010 11:19:57 -0800, Keith Goodman wrote: [clip] My guess is that having separate underlying functions for each dtype, ndim, and axis would be a nightmare for a large project like Numpy. But manageable for a focused project like nanny. Might be easier to migrate the nan* functions to using Ufuncs. Unless I'm missing something, np.nanmax - np.fmax.reduce np.nanmin - np.fmin.reduce For `nansum`, we'd need to add an ufunc `nanadd`, and for `nanargmax/min`, we'd need `argfmin/fmax'. How about that! I wasn't aware of fmax/fmin. Yes, I'd like a nanadd, please. arr = np.random.rand(1000, 1000) arr[arr 0.5] = np.nan np.nanmax(arr) 0.4625409581072 np.fmax.reduce(arr, axis=None) snip TypeError: an integer is required np.fmax.reduce(np.fmax.reduce(arr, axis=0), axis=0) 0.4625409581072 timeit np.fmax.reduce(np.fmax.reduce(arr, axis=0), axis=0) 100 loops, best of 3: 12.7 ms per loop timeit np.nanmax(arr) 10 loops, best of 3: 39.6 ms per loop timeit np.nanmax(arr, axis=0) 10 loops, best of 3: 46.5 ms per loop timeit np.fmax.reduce(arr, axis=0) 100 loops, best of 3: 12.7 ms per loop Cython is faster than np.fmax.reduce. I wrote a cython version of np.nanmax, called nanmax below. (It only handles the 2d, float64, axis=None case, but since the array is large I don't think that explains the time difference). Note that fmax.reduce is slower than np.nanmax when there are no NaNs: arr = np.random.rand(1000, 1000) timeit np.nanmax(arr) 100 loops, best of 3: 5.82 ms per loop timeit np.fmax.reduce(np.fmax.reduce(arr)) 100 loops, best of 3: 9.14 ms per loop timeit nanmax(arr) 1000 loops, best of 3: 1.17 ms per loop arr[arr 0.5] = np.nan timeit np.nanmax(arr) 10 loops, best of 3: 45.5 ms per loop timeit np.fmax.reduce(np.fmax.reduce(arr)) 100 loops, best of 3: 12.7 ms per loop timeit nanmax(arr) 1000 loops, best of 3: 1.17 ms per loop There seem to be some odd hardware/compiler dependencies. I get quite a different pattern of times: In [1]: arr = np.random.rand(1000, 1000) In [2]: timeit np.nanmax(arr) 100 loops, best of 3: 10.4 ms per loop In [3]: timeit np.fmax.reduce(arr.flat) 100 loops, best of 3: 2.09 ms per loop In [4]: arr[arr 0.5] = np.nan In [5]: timeit np.nanmax(arr) 100 loops, best of 3: 12.9 ms per loop In [6]: timeit np.fmax.reduce(arr.flat) 100 loops, best of 3: 7.09 ms per loop I've tweaked fmax with the reduce loop option but the nanmax times don't look like yours at all. I'm also a bit surprised that you don't see any difference in times when the array contains a lot of nans. I'm running on AMD Phenom, gcc 4.4.5. Ubuntu 10.04 64 bit, numpy 1.4.1. Difference in which times? nanny.nanmax with and wintout NaNs? The code doesn't explictily check for NaNs (it does check for all NaNs). It basically loops through the data and does: allnan = 1 ai = ai[i,k] if ai amax: amax = ai allnan = 0 I should make a benchmark suite. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] [ANN] Nanny, faster NaN functions
On Fri, Nov 19, 2010 at 7:51 PM, josef.p...@gmail.com wrote: does this give you the correct answer? 1np.nan False What's the starting value for amax? -inf? Because 1 np.nan is False, the current running max does not get updated, which is what we want. import nanny as ny np.nanmax([1, np.nan]) 1.0 np.nanmax([np.nan, 1]) 1.0 np.nanmax([np.nan, 1, np.nan]) 1.0 Starting value is -np.inf for floats and stuff like this for ints: cdef np.int32_t MININT32 = np.iinfo(np.int32).min cdef np.int64_t MININT64 = np.iinfo(np.int64).min Numpy does this: np.nanmax([]) snip ValueError: zero-size array to ufunc.reduce without identity Nanny does this: ny.nanmax([]) nan So I haven't taken care of that corner case yet. I'll commit nanmax to github in case anyone wants to give it a try. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] [ANN] Nanny, faster NaN functions
On Fri, Nov 19, 2010 at 8:05 PM, Charles R Harris charlesr.har...@gmail.com wrote: This doesn't look right: @cython.boundscheck(False) @cython.wraparound(False) def nanmax_2d_float64_axisNone(np.ndarray[np.float64_t, ndim=2] a): nanmax of 2d numpy array with dtype=np.float64 along axis=None. cdef Py_ssize_t i, j cdef int arow = a.shape[0], acol = a.shape[1], allnan = 1 cdef np.float64_t amax = 0, aij for i in range(arow): for j in range(acol): aij = a[i,j] if aij == aij: amax += aij allnan = 0 if allnan == 0: return np.float64(amax) else: return NAN It's doing a sum, not a comparison. That was a placeholder. Looks at the latest commit. Sorry for the confusion. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] [ANN] Nanny, faster NaN functions
On Fri, Nov 19, 2010 at 8:33 PM, josef.p...@gmail.com wrote: -np.inf-np.inf False If the only value is -np.inf, you will return nan, I guess. np.nanmax([-np.inf, np.nan]) -inf That's a great corner case. Thanks, Josef. This looks like it would fix it: change if ai amax: amax = ai to if ai = amax: amax = ai ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Returning numpy scalars in cython functions
The cython function below returns a long int: @cython.boundscheck(False) def mysum(np.ndarray[np.int64_t, ndim=1] a): sum of 1d numpy array with dtype=np.int64. cdef Py_ssize_t i cdef int asize = a.shape[0] cdef np.int64_t asum = 0 for i in range(asize): asum += a[i] return asum What's the best way to make it return a numpy long int, or whatever it is called, that has dtype, ndim, size, etc. class methods? The only thing I could come up with is changing the last line to return np.array(asum)[()] It works. And adds some overhead: a = np.arange(10) timeit mysum(a) 1000 loops, best of 3: 167 ns per loop timeit mysum2(a) 100 loops, best of 3: 984 ns per loop And for scale: timeit np.sum(a) 10 loops, best of 3: 3.3 us per loop I'm new to cython. Did I miss any optimizations in the mysum function above? ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Returning numpy scalars in cython functions
On Thu, Nov 18, 2010 at 10:08 AM, Francesc Alted fal...@pytables.org wrote: A Thursday 18 November 2010 18:51:04 Keith Goodman escrigué: What's the best way to make it return a numpy long int, or whatever it is called, that has dtype, ndim, size, etc. class methods? The only thing I could come up with is changing the last line to return np.array(asum)[()] Perhaps the scalar constructor is your best bet: type(np.array(2)[()]) type 'numpy.int64' type(np.int_(2)) type 'numpy.int64' timeit np.array(2)[()] 100 loops, best of 3: 791 ns per loop timeit np.int_(2) 100 loops, best of 3: 234 ns per loop Perfect! Thank you. a = np.arange(10) timeit mysum2(a) 100 loops, best of 3: 1.16 us per loop timeit mysum2_francesc(a) 100 loops, best of 3: 451 ns per loop I also added @cython.wraparound(False). ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] getitem and slice
There is more than one way to form the same slice. For example, a[:2] and a[0:2] and a[0:2:] pass slice objects to getitem with the same start, stop and step. Is there any way to get a hold of the exact character sequence the user used to form the slice? That is, I'd like to know if the user entered :2, 0:2 or 0:2:. Context: Say I have a 2d array-like object with axis 0 named 'space' and axis 1 named 'time', I'd like to be able to make the following distinctions: a['time':2] -- a[:,2] a['time':2:] -- a[:,2:] a['space':2] -- a[2] a['space':2:] -- a[2:] ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] speed of array vs matrix
On Mon, Oct 25, 2010 at 6:48 AM, Citi, Luca lc...@essex.ac.uk wrote: Hello, I have noticed a significant speed difference between the array and the matrix implementation of the dot product, especially for not-so-big matrices. For example: In [1]: import numpy as np In [2]: b = np.random.rand(104,1) In [3]: bm = np.mat(b) In [4]: a = np.random.rand(8, 104) In [5]: am = np.mat(a) In [6]: %timeit np.dot(a, b) 100 loops, best of 3: 1.74 us per loop In [7]: %timeit am * bm 10 loops, best of 3: 6.38 us per loop The results for two different PCs (PC1 with windows/EPD6.2-2 and PC2 with ubuntu/numpy-1.3.0) and two different sizes are below: array matrix 8x104 * 104x1 PC1 1.74us 6.38us PC2 1.23us 5.85us 8x10 * 10x5 PC1 2.38us 7.55us PC2 1.56us 6.01us For bigger matrices the timings seem to asymptotically approach. Is it something worth trying to fix or should I just accept this as a fact and, when working with small matrices, stick to array? I think the fixed overhead comes from the subclassing of arrays. The subclassing is done in Python and if an operation creates a matrix then __array_finalize__ is called. All that adds up to overhead. http://github.com/numpy/numpy/blob/master/numpy/matrixlib/defmatrix.py I wrote a mean-variance optimizer with matrices. Switching to arrays gave me a big speed up. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] DataArray fixes
On Wed, Oct 13, 2010 at 6:02 AM, Lluís xscr...@gmx.net wrote: Fernando Perez writes: By the way Lluis, I replied to one of your requests (Keith merged the other one, thanks!) here: How can I see which one got merged without checking one by one? It was this one: http://github.com/fperez/datarray/commit/a66b871af0ac7d50ca7a4f4094875b588074a3e9 I added this unit test: http://github.com/fperez/datarray/commit/1015c1906f1c20de51f9f68c9ce9d605f230aa46 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
