Mailman 3 January 2011 - NumPy-Discussion

Autosummary using numpydoc
by Matthew Turk Jan. 12, 2011

Jan. 12, 2011

Hi there, I've been trying to take the numpy docstring and apply the same methodology to a different project I work on, but there are a couple details that I think I'm unclear on, and I was hoping for some pointers or at least reassurances that it's working as intended, despite Sphinx's protests. The process of applying the numpy docstring method seems to take a few steps: 0) Have sphinx >=1.0 installed. 1) Write a numpydoc-compliant docstring in a class and its methods. 2) Make available to sphinx the numpydoc extension (in this case, I am using the current numpy git tip, which self-reports as 0.4) 3) Add the sphinx.ext.autosummary and numpydoc extensions to conf.py in the appropriate sphinx project 4) Copy the autosummary/class.rst template to the appropriate _template directory in the sphinx project 5) Set up the autosummary_generate variable in conf.py to contain all the files containing autosummary directives. 6) Build docs As a few other quick notes, my reading of the overridden class.rst template is that it comments out the individual method inclusions, using the HACK (comment) directive. The build mostly proceeds correctly when these steps have been taken, but a copious number of warnings are emitted. (The same character of warnings are emitted by the numpy documentation build.) These mostly show up, for instance, as: WARNING: toctree contains reference to nonexisting document 'reference/api/generated/add_phase_object' In this case, add_phase_object is a method hanging off an autosummary'ed class. It looks like a warning gets emitted for every method on every autosummary'ed class -- this can number into the thousands very easily. Additionally, warnings that look like: reference/api/generated/yt.visualization.api.PlotCollection.add_phase_object.rst:: WARNING: document isn't included in any toctree show up, where these reflect the full method name and class name. It looks like the toctree is gaining references to the method names without prefixing them with the class names, but the generated docs all have the full name resolution in their filenames. (There also seems to be a glitch in the output when using the ".. HACK" class.rst template, as within my code it appears to strip the one-line descriptions from the methods.) Is this behaving as expected, or have I perhaps gone through the steps wrong? Having so many warnings can confuse debugging other portions of the build. I appreciate any suggestions you might have -- thanks very much! Best, Matt

1 0

Variable in an array name?
by dstaley Jan. 12, 2011

Jan. 12, 2011

Is it possible to use a variable in an array name? I am looping through a bunch of calculations, and need to have each array as a separate entity. I'm pretty new to python and numpy, so forgive my ignorance. I'm sure there is a simple answer, but I can't seem to find it. let's say i have a variable 'i': i = 5 I would like my array to have the name array5 I know how I could do this manually, but not in a loop where i is redefined several times. any thoughts/comments/suggestions are appreciated. Thanks. -DS -- View this message in context: http://old.nabble.com/Variable-in-an-array-name--tp30645276p30645276.html Sent from the Numpy-discussion mailing list archive at Nabble.com.

4 5

mapping a function to a masked array
by Ben Elliston Jan. 12, 2011

Jan. 12, 2011

I have a masked array of values that I would like to transform through a user-defined function. Naturally, I want to ignore any values that are masked in the initial array. The user-defined function examines other points around the value in question, so I need to use ndenumerate (or similar) to get the array index as I iterate over the array. So, I have two questions: how to make this run without looping in Python, and how to avoid masked values. Here is the clunky solution I have so far: result = ma.copy (data) for i, val in ndenumerate (data): if not data.mask[i]: result[i] = myfunc (data, i, val) Any suggestions? Thanks, Ben

2 1

histogram2d and decreasing bin edges
by Matthias Frank Jan. 11, 2011

Jan. 11, 2011

Hi all, I've noticed a change in numpy.histogram2d between (possibly very much) older versions and the current one: The function can no longer handle the situation where bin edges decrease instead of increasing monotonically. The reason for this seems to be the handling of outliers histogramdd, see the output of minimal example below. If I understand correctly, this is the only place where histogramdd implicitly assumes monotonically increasing bin edges. If so, this could be fixed to work with increasing and decreasing bin edges by taking abs(dedges[i]).min() when calculating the rounding precision. If not, it might be more consistent, and produce a more meaningful error message, if histogram2d asserted that bin edges increase monotonically and otherwise raised an AttributeError as the 1-d histogram() function does in that case (see below) Matthias In [1]: import numpy In [2]: numpy.__version__ Out[2]: '1.5.1' In [3]: ascending=numpy.array([0,1]) In [4]: descending=numpy.array([1,0]) In [5]: numpy.histogram2d([0.5],[0.5],bins=(ascending,ascending)) Out[5]: (array([[ 1.]]), array([ 0., 1.]), array([ 0., 1.])) In [6]: numpy.histogram2d([0.5],[0.5],bins=(descending,descending)) Warning: invalid value encountered in log10 --------------------------------------------------------------------------- ValueError Traceback (most recent call last) /lib/python2.6/site-packages/numpy/lib/twodim_base.pyc in histogram2d(x, y, bins, range, normed, weights) 613 xedges = yedges = asarray(bins, float) 614 bins = [xedges, yedges] --> 615 hist, edges = histogramdd([x,y], bins, range, normed, weights) 616 return hist, edges[0], edges[1] 617 /lib/python2.6/site-packages/numpy/lib/function_base.pyc in histogramdd(sample, bins, range, normed, weights) 312 for i in arange(D): 313 # Rounding precision --> 314 decimal = int(-log10(dedges[i].min())) +6 315 # Find which points are on the rightmost edge. 316 on_edge = where(around(sample[:,i], decimal) == around(edges[i][-1], ValueError: cannot convert float NaN to integer Behavior of the 1-d histogram() In [8]: numpy.histogram([0.5],bins=ascending) Out[8]: (array([1]), array([0, 1])) In [9]: numpy.histogram([0.5],bins=descending) --------------------------------------------------------------------------- AttributeError Traceback (most recent call last) /lib/python2.6/site-packages/numpy/lib/function_base.pyc in histogram(a, bins, range, normed, weights) 160 if (np.diff(bins) < 0).any(): 161 raise AttributeError( --> 162 'bins must increase monotonically.') 163 164 # Histogram is an integer or a float array depending on the weights. AttributeError: bins must increase monotonically.

1 0

Rolling window (moving average, moving std, and more)
by Erik Rigtorp Jan. 11, 2011

Jan. 11, 2011

Hi, 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 <kwgoodman(a)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. import numpy as np def rolling_window(a, window): """ Make an ndarray with a rolling window of the last dimension Parameters ---------- a : array_like Array to add rolling window to window : int Size of rolling window Returns ------- Array that is a view of the original array with a added dimension of size w. Examples -------- >>> x=np.arange(10).reshape((2,5)) >>> rolling_window(x, 3) array([[[0, 1, 2], [1, 2, 3], [2, 3, 4]], [[5, 6, 7], [6, 7, 8], [7, 8, 9]]]) Calculate rolling mean of last dimension: >>> np.mean(rolling_window(x, 3), -1) array([[ 1., 2., 3.], [ 6., 7., 8.]]) """ if window < 1: raise ValueError, "`window` must be at least 1." if window > a.shape[-1]: raise ValueError, "`window` is too long." shape = a.shape[:-1] + (a.shape[-1] - window + 1, window) strides = a.strides + (a.strides[-1],) return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides) Using np.swapaxes(-1, axis) rolling aggregations over any axis can be computed. I submitted a pull request to add this to the stride_tricks module. Erik

4 12

Drawing circles in a numpy array
by totonixsame＠gmail.com Jan. 11, 2011

Jan. 11, 2011

Hi all, I have this problem: Given some point draw a circle centered in this point with radius r. I'm doing that using numpy this way (Snippet code from here [1]): >>> # Create the initial black and white image >>> import numpy as np >>> from scipy import ndimage >>> a = np.zeros((512, 512)).astype(uint8) #unsigned integer type needed by watershed >>> y, x = np.ogrid[0:512, 0:512] >>> m1 = ((y-200)**2 + (x-100)**2 < 30**2) >>> m2 = ((y-350)**2 + (x-400)**2 < 20**2) >>> m3 = ((y-260)**2 + (x-200)**2 < 20**2) >>> a[m1+m2+m3]=1 >>> imshow(a, cmap = cm.gray)# left plot in the image above The problem is that it have to evaluate all values from 0 to image size (in snippet, from 0 to 512 in X and Y dimensions). There is a faster way of doing that? Without evaluate all that values? For example: only evaluate from 0 to 30, in a circle centered in (0, 0) with radius 30. Thanks! Thiago Franco de Moraes [1] - http://www.scipy.org/Cookbook/Watershed

1 1

numexpr with the new iterator
by Mark Wiebe Jan. 11, 2011

Jan. 11, 2011

As a benchmark of C-based iterator usage and to make it work properly in a multi-threaded context, I've updated numexpr to use the new iterator. In addition to some performance improvements, this also made it easy to add optional out= and order= parameters to the evaluate function. The numexpr repository with this update is available here: https://github.com/m-paradox/numexpr To use it, you need the new_iterator branch of NumPy from here: https://github.com/m-paradox/numpy In all cases tested, the iterator version of numexpr's evaluate function matches or beats the standard version. The timing results are below, with some explanatory comments placed inline: -Mark In [1]: import numexpr as ne # numexpr front page example In [2]: a = np.arange(1e6) In [3]: b = np.arange(1e6) In [4]: timeit a**2 + b**2 + 2*a*b 1 loops, best of 3: 121 ms per loop In [5]: ne.set_num_threads(1) # iterator version performance matches standard version In [6]: timeit ne.evaluate("a**2 + b**2 + 2*a*b") 10 loops, best of 3: 24.8 ms per loop In [7]: timeit ne.evaluate_iter("a**2 + b**2 + 2*a*b") 10 loops, best of 3: 24.3 ms per loop In [8]: ne.set_num_threads(2) # iterator version performance matches standard version In [9]: timeit ne.evaluate("a**2 + b**2 + 2*a*b") 10 loops, best of 3: 21 ms per loop In [10]: timeit ne.evaluate_iter("a**2 + b**2 + 2*a*b") 10 loops, best of 3: 20.5 ms per loop # numexpr front page example with a 10x bigger array In [11]: a = np.arange(1e7) In [12]: b = np.arange(1e7) In [13]: ne.set_num_threads(2) # the iterator version performance improvement is due to # a small task scheduler tweak In [14]: timeit ne.evaluate("a**2 + b**2 + 2*a*b") 1 loops, best of 3: 282 ms per loop In [15]: timeit ne.evaluate_iter("a**2 + b**2 + 2*a*b") 1 loops, best of 3: 255 ms per loop # numexpr front page example with a Fortran contiguous array In [16]: a = np.arange(1e7).reshape(10,100,100,100).T In [17]: b = np.arange(1e7).reshape(10,100,100,100).T In [18]: timeit a**2 + b**2 + 2*a*b 1 loops, best of 3: 3.22 s per loop In [19]: ne.set_num_threads(1) # even with a C-ordered output, the iterator version performs better In [20]: timeit ne.evaluate("a**2 + b**2 + 2*a*b") 1 loops, best of 3: 3.74 s per loop In [21]: timeit ne.evaluate_iter("a**2 + b**2 + 2*a*b") 1 loops, best of 3: 379 ms per loop In [22]: timeit ne.evaluate_iter("a**2 + b**2 + 2*a*b", order='C') 1 loops, best of 3: 2.03 s per loop In [23]: ne.set_num_threads(2) # the standard version just uses 1 thread here, I believe # the iterator version performs the same as for the flat 1e7-sized array In [24]: timeit ne.evaluate("a**2 + b**2 + 2*a*b") 1 loops, best of 3: 3.92 s per loop In [25]: timeit ne.evaluate_iter("a**2 + b**2 + 2*a*b") 1 loops, best of 3: 254 ms per loop In [26]: timeit ne.evaluate_iter("a**2 + b**2 + 2*a*b", order='C') 1 loops, best of 3: 1.74 s per loop

3 11

Re: [Numpy-discussion] indexing of rank-0 structured arrays: why not?
by Nils Becker Jan. 10, 2011

Jan. 10, 2011

Robert, your answer does work: after indexing with () I can then further index into the datatype. In [115]: a_rank_0[()][0] Out[115]: 0.0 I guess I just found the fact confusing that a_rank_1[0] and a_rank_0 compare and print equal but behave differently under indexing. More precisely if I do In [117]: b = a_rank_1[0] then In [118]: b.shape Out[118]: () and In [120]: a_rank_0 == b Out[120]: True but In [119]: b[0] Out[119]: 0.0 works but a_rank_0[0] doesn't. I thought b is a rank-0 array which it apparently is not since it can be indexed. So maybe b[0] should fail for consistency? N.

2 1

indexing of rank-0 structured arrays: why not?
by Nils Becker Jan. 10, 2011

Jan. 10, 2011

Hi, I noticed that I can index into a dtype when I take an element of a rank-1 array but not if I make a rank-0 array directly. This seems inconsistent. A bug? Nils In [76]: np.version.version Out[76]: '1.5.1' In [78]: dt = np.dtype([('x', '<f8'), ('y', '<f8')]) In [80]: a_rank_1 = np.zeros((1,), dtype=dt) In [81]: a_rank_0 = np.zeros((), dtype=dt) In [83]: a_rank_1[0] Out[83]: (0.0, 0.0) In [84]: a_rank_1[0] == a_rank_0 Out[84]: True In [85]: a_rank_1[0][0] Out[85]: 0.0 In [86]: a_rank_0[0] --------------------------------------------------------------------------- IndexError Traceback (most recent call last) /home/pabra/ramp/testramp/<ipython console> in <module>() IndexError: 0-d arrays can't be indexed In [87]: a_rank_0['x'] Out[87]: array(0.0)

2 1

Re: [Numpy-discussion] speed of numpy.ndarray compared to Numeric.array
by EMMEL Thomas Jan. 10, 2011

Jan. 10, 2011

To John: > Did you try larger arrays/tuples? I would guess that makes a significant > difference. No I didn't, due to the fact that these values are coordinates in 3D (x,y,z). In fact I work with a list/array/tuple of arrays with 100000 to 1M of elements or more. What I need to do is to calculate the distance of each of these elements (coordinates) to a given coordinate and filter for the nearest. The brute force method would look like this: #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ def bruteForceSearch(points, point): minpt = min([(vec2Norm(pt, point), pt, i) for i, pt in enumerate(points)], key=itemgetter(0)) return sqrt(minpt[0]), minpt[1], minpt[2] #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ def vec2Norm(pt1,pt2): xDis = pt1[0]-pt2[0] yDis = pt1[1]-pt2[1] zDis = pt1[2]-pt2[2] return xDis*xDis+yDis*yDis+zDis*zDis I have a more clever method but it still takes a lot of time in the vec2norm-function. If you like I can attach a running example. To Ben: > Don't know how much of an impact it would have, but those timeit statements > for array creation include the import process, which are going to be > different for each module and are probably not indicative of the speed of > array creation. No, the timeit statements counts the time for the statement in the first argument only, the import-thing isn't included in the time. Thomas This email and any attachments are intended solely for the use of the individual or entity to whom it is addressed and may be confidential and/or privileged. If you are not one of the named recipients or have received this email in error, (i) you should not read, disclose, or copy it, (ii) please notify sender of your receipt by reply email and delete this email and all attachments, (iii) Dassault Systemes does not accept or assume any liability or responsibility for any use of or reliance on this email.For other languages, go to http://www.3ds.com/terms/email-disclaimer.

5 8