Mailman 3 February 2012 - NumPy-Discussion

Re: [Numpy-discussion] Index Array Performance
by Marcel Oliver Feb. 14, 2012

Feb. 14, 2012

Francesc Alted wrote: > On Feb 14, 2012, at 1:50 AM, Wes McKinney wrote: > [clip] > > But: > > > > In [40]: timeit hist[i, j] > > 10000 loops, best of 3: 32 us per loop > > > > So that's roughly 7-8x slower than a simple Cython method, so I > > sincerely hope it could be brought down to the sub 10 microsecond > > level with a little bit of work. > > I vaguely remember this has shown up before. My hunch is that > indexing in NumPy is so powerful, that it has to check for a lot of > different values for indices (integers, tuples, lists, arrays?), > and that it is all these checks what is taking time. Your Cython > wrapper just assumed that the indices where integers, so this is > probably the reason why it is that much faster. Thanks for all the replies. Playing a bit with timeit, it is clear that it cannot just be the overhead of checking the type of the index array, as the overhead grows very roughly propertional to the size of the index array, but remains independent of the size of the indexed array. In [1]: a=arange(1000000) In [2]: i = arange(10) In [3]: timeit a[i] 100000 loops, best of 3: 1.95 us per loop In [4]: i = arange(100) In [5]: timeit a[i] 100000 loops, best of 3: 4.07 us per loop In [6]: i = arange(1000) In [7]: timeit a[i] 10000 loops, best of 3: 23.3 us per loop In [8]: timeit i+i 100000 loops, best of 3: 5.62 us per loop In [9]: i = arange(10000) In [10]: timeit a[i] 1000 loops, best of 3: 220 us per loop In [11]: timeit i+i 10000 loops, best of 3: 28.1 us per loop It would really be nice if this could be made more efficient as it's a very nice and powerful construct. Robert Kern wrote: > Other people have explained that yes, applying index arrays is slow. I > would just like to add the tangential point that this code does not > behave the way that you think it does. You cannot make histograms like > this. The statement "hist[i,j] += 1" gets broken down into three > separate statements by the Python compiler: > > tmp = hist.__getitem__((i,j)) > tmp = tmp.__iadd__(1) > hist.__setitem__((i,j), tmp) > > Note that tmp is a new array with copies of the data in hist at the > (i,j) locations, possibly multiple copies if the i index has > repetitions. Each one of these copies gets incremented by 1, then the > __setitem__() will apply each of those in turn to the appropriate cell > in hist, each one simply overwriting the previous one. I know. The posted code is correct because it updates only one bin in each of the many histograms at a time, which is natural for that particular example. It would indeed be nice if the increment construct would result in multiple increments work when indices repeat - I think that would be generally more useful and expected. But current behavior is documented, so it's fine if only it was efficient... Regards, Marcel

2 1

Index Array Performance
by Marcel Oliver Feb. 14, 2012

Feb. 14, 2012

Hi, I have a short piece of code where the use of an index array "feels right", but incurs a severe performance penalty: It's about an order of magnitude slower than all other operations with arrays of that size. It comes up in a piece of code which is doing a large number of "on the fly" histograms via hist[i,j] += 1 where i is an array with the bin index to be incremented and j is simply enumerating the histograms. I attach a full short sample code below which shows how it's being used in context, and corresponding timeit output from the critical code section. Questions: - Is this a matter of principle, or due to an inefficient implementation? - Is there an equivalent way of doing it which is fast? Regards, Marcel ================================================================= #! /usr/bin/env python # Plot the bifurcation diagram of the logistic map from pylab import * Nx = 800 Ny = 600 I = 50000 rmin = 2.5 rmax = 4.0 ymin = 0.0 ymax = 1.0 rr = linspace (rmin, rmax, Nx) x = 0.5*ones(rr.shape) hist = zeros((Ny+1,Nx), dtype=int) j = arange(Nx) dy = ymax/Ny def f(x): return rr*x*(1.0-x) for n in xrange(1000): x = f(x) for n in xrange(I): x = f(x) i = array(x/dy, dtype=int) hist[i,j] += 1 figure() imshow(hist, cmap='binary', origin='lower', interpolation='nearest', extent=(rmin,rmax,ymin,ymax), norm=matplotlib.colors.LogNorm()) xlabel ('$r$') ylabel ('$x$') title('Attractor of the logistic map $x_{n+1} = r \, x_n (1-x_n)$') show() ==================================================================== In [4]: timeit y=f(x) 10000 loops, best of 3: 19.4 us per loop In [5]: timeit i = array(x/dy, dtype=int) 10000 loops, best of 3: 22 us per loop In [6]: timeit img[i,j] += 1 10000 loops, best of 3: 119 us per loop

5 8

Typecasting changes from 1.5.1 to 1.6.1
by Travis Oliphant Feb. 14, 2012

Feb. 14, 2012

Here is the code I used to determine the coercion table of types. I first used *all* of the numeric_ops, narrowed it down to those with 2 inputs and 1 output, and then determined the run-time coercion table. Then, I removed ops that had the same tables until I was left with binary ops that had different coercion tables. Some operations were NotImplemented and I used 'X' in the table for those combinations. The table for each op is a dictionary with keys given by (type1, type2) and values given by a length-4 list of the types of the result between: [scalar-scalar, scalar-array, array-scalar, array-array] where the first term is type1 and the second term is type2. This resulting dictionary of tables for each op is then saved to a file. I ran this code for NumPy 1.5.1 64-bit and then again for NumPy 1.6.1 64-bit. I also ran this code for NumPy 1.4.1 64-bit and NumPy 1.3.1.dev 64-bit. The code to compare them is also attached. I'm attaching also the changes that have occurred between 1.3.1.dev and 1.4.1, 1.4.1 to 1.5.1, and finally 1.5.1 to 1.6.1 As you can see there were changes in each release. Most of these were minor prior to the change from 1.5.1 to 1.6.1. I am still reviewing the changes from 1.5.1 to 1.6.1. At first blush, it looks like there are a lot of changes to swallow that are not necessarily minor. I really would like to just say all is well, and it's no big deal. I hope that users really don't care and nobody's code is really relying on array-scalar combination conversions. -Travis

4 5

Change in scalar upcasting rules for 1.6.x?
by Matthew Brett Feb. 14, 2012

Feb. 14, 2012

Hi, I recently noticed a change in the upcasting rules in numpy 1.6.0 / 1.6.1 and I just wanted to check it was intentional. For all versions of numpy I've tested, we have: >>> import numpy as np >>> Adata = np.array([127], dtype=np.int8) >>> Bdata = np.int16(127) >>> (Adata + Bdata).dtype dtype('int8') That is - adding an integer scalar of a larger dtype does not result in upcasting of the output dtype, if the data in the scalar type fits in the smaller. For numpy < 1.6.0 we have this: >>> Bdata = np.int16(128) >>> (Adata + Bdata).dtype dtype('int8') That is - even if the data in the scalar does not fit in the dtype of the array to which it is being added, there is no upcasting. For numpy >= 1.6.0 we have this: >>> Bdata = np.int16(128) >>> (Adata + Bdata).dtype dtype('int16') There is upcasting... I can see why the numpy 1.6.0 way might be preferable but it is an API change I suppose. Best, Matthew

10 30

addition to numpy discussion list
by pulkit yadav Feb. 14, 2012

Feb. 14, 2012

Hello, I am a Python enthusiast and developer. Please add me to numpy mailing list so that I can contribute to the FLOSS community. -- Pulkit

2 1

Re: [Numpy-discussion] @Dag re numpy.pxd
by Gordon L. Kindlmann Feb. 14, 2012

Feb. 14, 2012

Hello, This (below) caught my eye and I'm wondering what further information is available? I very much value the ability to wrap underlying array data from numpy for processing in non-python libraries, as well as the ability to wrap numpy arrays around array data allocated by non-python libraries. Is this capability going to be removed? Gordon > Message: 6 > Date: Sat, 11 Feb 2012 13:31:51 -0700 > From: Charles R Harris <charlesr.harris(a)gmail.com> > Subject: [Numpy-discussion] @Dag re numpy.pxd > To: numpy-discussion <numpy-discussion(a)scipy.org> > Message-ID: > <CAB6mnxJVERqK64UpG3GOmaWDQHUn-z-=WJMQju0iU4KO9sefTg(a)mail.gmail.com> > Content-Type: text/plain; charset="iso-8859-1" > > Hi Dag, > > This probably needs to be on the cython mailing list at some point, but I > thought I'd start the discussion here. Numpy is going to begin deprecating > direct access to ndarray/dtype internals, ala arr->data etc. There are > currently macros/functions for many of these operations in the numpy > development branch and I expect more to go in over the coming year. Also, > some of the macros have been renamed. I don't know the best way for Cython > to support this, but the current version (0.15 here) generates code that > will fail if the deprecated things are excluded. Ideally, numpy.pxd would > have numpy version dependent parts but I don't know if that is possible. In > any case, I'd like your thoughts on the best way to coordinate this > migration with Cython. > > Chuck

3 2

_import_array()
by Mads Ipsen Feb. 14, 2012

Feb. 14, 2012

Hi, I have C++ module (OpenGL) that extracts data from numpy arrays. The interface is pure read-only: It never returns any Python objects but only extracts data from numpy arrays. Eg: #include "numpy/arrayobject.h" void PrimitiveManager::deleteAtoms(PyObject * numpy_indices) { // Extract number of indices int const n = static_cast<int>(PyArray_DIMS(numpy_indices)[0]); long * const indices = (long *) PyArray_DATA(numpy_indices); // Delete atoms in buffer for (int i = 0; i < n; ++i) { // Do stuff } } Now, when I compile the code with g++, I get the following warning: numpy/core/include/numpy/__multiarray_api.h:1532: warning: 'int _import_array()' defined but not used Do I need to call '_import_array()' somewhere? Am I doing something potentially nasty? Best regards, Mads -- +-----------------------------------------------------+ | Mads Ipsen | +----------------------+------------------------------+ | Gåsebæksvej 7, 4. tv | | | DK-2500 Valby | phone: +45-29716388 | | Denmark | email: mads.ipsen(a)gmail.com | +----------------------+------------------------------+

4 5

Fwd: Re: Creating parallel curves
by Andrea Gavana Feb. 14, 2012

Feb. 14, 2012

---------- Forwarded message ---------- From: "Andrea Gavana" <andrea.gavana(a)gmail.com> Date: Feb 13, 2012 11:31 PM Subject: Re: [Numpy-discussion] Creating parallel curves To: "Jonathan Hilmer" <jkhilmer(a)gmail.com> Thank you Jonathan for this, it's exactly what I was looking for. I' ll try it tomorrow on the 768 well trajectories I have and I'll let you know if I stumble upon any issue. If someone could shed some light on my problem number 2 (how to adjust the scaling/distance) so that the curves look parallel on a matplotlib graph even though the axes scales are different, I'd be more than grateful. Thank you in advance. Andrea. On Feb 13, 2012 4:32 AM, "Jonathan Hilmer" <jkhilmer(a)gmail.com> wrote: > Andrea, > > This is playing some tricks with 2D array expansion to make a tradeoff > in memory for speed. Given two sets of matching vectors (one > reference, given first, and a newly-expanded one, given second), it > removes all points from the expanded vectors that aren't needed to > describe the new contour. > > def filter_expansion(x, y, x_expan, y_expan, distance_target, tol=1e-6): > > target_xx, expansion_xx = scipy.meshgrid(x, x_expan) > target_yy, expansion_yy = scipy.meshgrid(y, y_expan) > > distance = scipy.sqrt((expansion_yy - target_yy)**2 + (expansion_xx > - > target_xx)**2) > > valid = distance.min(axis=1) > distance_target*(1.-tol) > > return x_expan.compress(valid), y_expan.compress(valid) > # > > > Jonathan > > > On Sun, Feb 12, 2012 at 2:31 PM, Robert Kern <robert.kern(a)gmail.com> > wrote: > > On Sun, Feb 12, 2012 at 20:26, Andrea Gavana <andrea.gavana(a)gmail.com> > wrote: > > > >> I know, my definition of "parallel" was probably not orthodox enough. > >> What I am looking for is to generate 2 curves that look "graphically > >> parallel enough" to the original one, and not "parallel" in the true > >> mathematical sense. > > > > There is a rigorous way to define the curve that you are looking for, > > and fortunately it gives some hints for implementation. For each point > > (x,y) in space, associate with it the nearest distance D from that > > point to the reference curve. The "parallel" curves are just two sides > > of the level set where D(x,y) is equal to the specified distance > > (possibly removing the circular caps that surround the ends of the > > reference curve). > > > > If performance is not a constraint, then you could just evaluate that > > D(x,y) function on a fine-enough grid and do marching squares to find > > the level set. matplotlib's contour plotting routines can help here. > > > > There is a hint in the PyX page that you linked to that you should > > consider. Angles in the reference curve become circular arcs in the > > "parallel" curves. So if your reference curve is just a bunch of line > > segments, then what you can do is take each line segment, and make > > parallel copies the same length to either side. Now you just need to > > connect up these parallel segments with each other. You do this by > > using circular arcs centered on the vertices of the reference curve. > > Do this on both sides. On the "outer" side, the arcs will go "forward" > > while on the "inner" side, the arcs will go "backwards" just like the > > cusps that you saw in your attempt. Now let's take care of that. You > > will have two self-intersecting curves consisting of alternating line > > segments and circular arcs. Parts of these curves will be too close to > > the reference curve. You will have to go through these curves to find > > the locations of self-intersection and remove the parts of the > > segments and arcs that are too close to the reference curve. This is > > tricky to do, but the formulae for segment-segment, segment-arc, and > > arc-arc intersection can be found online. > > > > -- > > Robert Kern > > > > "I have come to believe that the whole world is an enigma, a harmless > > enigma that is made terrible by our own mad attempt to interpret it as > > though it had an underlying truth." > > -- Umberto Eco > > _______________________________________________ > > NumPy-Discussion mailing list > > NumPy-Discussion(a)scipy.org > > http://mail.scipy.org/mailman/listinfo/numpy-discussion >

2 1

Re: [Numpy-discussion] numpy.arange() error?
by Matteo Malosio Feb. 14, 2012

Feb. 14, 2012

I think the problem is quite easy to solve, without changing the "documentation" behaviour. The doc says: Help on built-in function arange in module numpy.core.multiarray: / arange(...) arange([start,] stop[, step,], dtype=None) Return evenly spaced values within a given interval. Values are generated within the half-open interval ``[start, stop)`` (in other words, the interval including `start` but excluding `stop`). For integer arguments the function is equivalent to the Python built-in `range <http://docs.python.org/lib/built-in-funcs.html>`_ function, but returns a ndarray rather than a list. / stop is exclusive "by definition". So substracting a very small value to stop when processing "stop" I think is the best way. Matteo Il 10/02/2012 02:22, Drew Frank ha scritto: > On Thu, Feb 9, 2012 at 3:40 PM, Benjamin Root <ben.root(a)ou.edu > <mailto:[email protected]>> wrote: > > > > On Thursday, February 9, 2012, Sturla Molden <sturla(a)molden.no > <mailto:[email protected]>> wrote: > > > > > > Den 9. feb. 2012 kl. 22:44 skrev eat <e.antero.tammi(a)gmail.com > <mailto:[email protected]>>: > > > >> > > Maybe this issue is raised also earlier, but wouldn't it be more > consistent to let arange operate only with integers (like Python's > range) and let linspace handle the floats as well? > > > > > > Perhaps. Another possibility would be to let arange take decimal > arguments, possibly entered as text strings. > > Sturla > > > Personally, I treat arange() to mean, "give me a sequence of > values from x to y, exclusive, with a specific step size". > Nowhere in that statement does it guarantee a particular number > of elements. Whereas linspace() means, "give me a sequence of > evenly spaced numbers from x to y, optionally inclusive, such that > there are exactly N elements". They complement each other well. > > > I agree -- both functions are useful and I think about them the same > way. The unfortunate part is that tiny precision errors in y can make > arange appear to be "sometimes-exclusive" rather than always > exclusive. I've always imagined there to be a sort of duality between > the two functions, where arange(low, high, step) == linspace(low, > high-step, round((high-low)/step)) in cases where (high - low)/step is > integral, but it turns out this is not the case. > > > There are times when I intentionally will specify a range where > the step size will not nicely fit. i.e.- np.arange(1, 7, 3.5). I > wouldn't want this to change. > > > Nor would I. What I meant to express earlier is that I like how > Matlab addresses this particular class of floating point precision > errors, not that I think arange output should somehow include both > endpoints. > > > My vote is that if users want matlab-colon-like behavior, we could > make a new function - maybe erange() for "exact range"? > > > Ben Root > > > That could work; it would completely replace arange for me in every > circumstance I can think of, but I understand we can't just go > changing the behavior of core functions. > > Drew > > > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion(a)scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion -- ------------------------------------------------------- Matteo Malosio, Eng. Researcher ITIA-CNR (www.itia.cnr.it) Institute of Industrial Technologies and Automation National Research Council via Bassini 15, 20133 MILANO, ITALY Ph: +39 0223699625 Fax: +39 0223699925 e-mail:[email protected] -------------------------------------------------------

1 0

Initializing an array to a constant value
by Pierre Haessig Feb. 13, 2012

Feb. 13, 2012

I have a pretty silly question about initializing an array a to a given scalar value, say A. Most of the time I use a=np.ones(shape)*A which seems the most widespread idiom, but I got recently interested in getting some performance improvement. I tried a=np.zeros(shape)+A, based on broadcasting but it seems to be equivalent in terms of speed. Now, the fastest : a = np.empty(shape) a.fill(A) but it is a two-steps instruction to do one thing, which I feel doesn't look very nice. Did I miss an all-in-one function like numpy.fill(shape, A) ? Best, Pierre

2 2