Mailman 3 May 2011 - NumPy-Discussion

convert integer into bit array
by Nikolas Tautenhahn May 17, 2011

May 17, 2011

Hi, for some research, I need to convert lots of integers into their bit representation - but as a bit array, not a string like numpy.binary_repr() returns it. So instead of In [22]: numpy.binary_repr(23) Out[22]: '10111 I'd need: numpy.binary_magic(23) Out: array([ True, False, True, True, True], dtype=bool) is there any way to do this efficiently? best regards, Nik

3 5

PyArray_ContiguousFromObject & PyDECREF
by Michiel Jan Laurens de Hoon May 17, 2011

May 17, 2011

Dear pythoneers, In Chapter 12 ("Writing a C extension") in the Numerical Python manual, in the first example PyArray_ContiguousFromObject is used together with PyDECREF, while in the last example PyArray_ContiguousFromObject is used but the arrays are not PyDECREF'ed before the function returns. From Chapter 13 ("C API Reference") it seems that a PyDECREF is needed in both cases. Or am I missing something here? Is PyDECREF not needed for some reason in the last example? Thanks in advance, Michiel de Hoon U Tokyo First example: ============== static PyObject * trace(PyObject *self, PyObject *args) { PyObject *input; PyArrayObject *array; double sum; int i, n; if (!PyArg_ParseTuple(args, "O", &input)) return NULL; array = (PyArrayObject *) PyArray_ContiguousFromObject(input, PyArray_DOUBLE, 2, 2); if (array == NULL) return NULL; n = array->dimensions[0]; if (n > array->dimensions[1]) n = array->dimensions[1]; sum = 0.; for (i = 0; i < n; i++) sum += *(double *)(array->data + i*array->strides[0] + i*array->strides[1]); Py_DECREF(array); return PyFloat_FromDouble(sum); } Last example: ============= static PyObject * matrix_vector(PyObject *self, PyObject *args) { PyObject *input1, *input2; PyArrayObject *matrix, *vector, *result; int dimensions[1]; double factor[1]; double real_zero[1] = {0.}; long int_one[1] = {1}; long dim0[1], dim1[1]; extern dgemv_(char *trans, long *m, long *n, double *alpha, double *a, long *lda, double *x, long *incx, double *beta, double *Y, long *incy); if (!PyArg_ParseTuple(args, "dOO", factor, &input1, &input2)) return NULL; matrix = (PyArrayObject *) PyArray_ContiguousFromObject(input1, PyArray_DOUBLE, 2, 2); if (matrix == NULL) return NULL; vector = (PyArrayObject *) PyArray_ContiguousFromObject(input2, PyArray_DOUBLE, 1, 1); if (vector == NULL) return NULL; if (matrix->dimensions[1] != vector->dimensions[0]) { PyErr_SetString(PyExc_ValueError, "array dimensions are not compatible"); return NULL; } dimensions[0] = matrix->dimensions[0]; result = (PyArrayObject *)PyArray_FromDims(1, dimensions, PyArray_DOUBLE); if (result == NULL) return NULL; dim0[0] = (long)matrix->dimensions[0]; dim1[0] = (long)matrix->dimensions[1]; dgemv_("T", dim1, dim0, factor, (double *)matrix->data, dim1, (double *)vector->data, int_one, real_zero, (double *)result->data, int_one); return PyArray_Return(result); } -- Michiel de Hoon, Assistant Professor University of Tokyo, Institute of Medical Science Human Genome Center 4-6-1 Shirokane-dai, Minato-ku Tokyo 108-8639 Japan http://bonsai.ims.u-tokyo.ac.jp/~mdehoon

3 2

Need to eliminate a nested loop
by Elfnor May 16, 2011

May 16, 2011

Hi The following code produces the desired result but has a slow triple loop iterating over the matrix multiplication. I'm sure it can be eliminated with a neat indexing trick but I can't figure out how. Any suggestions please? ----------------------------- import numpy #define domain of function x = numpy.linspace(-5,5,64) y = numpy.linspace(-5,5,64) z = numpy.linspace(-5,5,64) #calculate f at each point in domain a = 5.0 b = 3.0 c = 2.0 #ellipsoid E = numpy.array([[1/a**2, 0 , 0 , 0 ], [ 0 ,1/b**2 , 0 , 0 ], [ 0 , 0 ,1/c**2, 0 ], [ 0 , 0 , 0 , -1 ]]) f = numpy.zeros((x.size,y.size,z.size)) for i,xi in enumerate(x): for j,yj in enumerate(y): for k,zk in enumerate(z): X = numpy.array([xi,yj,zk,1]) f[i,j,k] = numpy.dot(numpy.dot(X.T,E),X) ----------------------------------- Thanks Eleanor -- View this message in context: http://old.nabble.com/Need-to-eliminate-a-nested-loop-tp31591457p31591457.h… Sent from the Numpy-discussion mailing list archive at Nabble.com.

2 2

Performance issues with numpy cross tabulation
by Justin Thomas May 16, 2011

May 16, 2011

Hello, I am having trouble with performance when trying to create a cross tabulation using numpy. Ideally, I would calculate each cell in the cross tabulation separately because this gives me the greatest amount of flexibility. I have included some sample code as a reference and am really looking for better approaches to the simpleLoop method. So far the histogram2d and histogramdd methods seem to outperform any code I write by a factor of about 100, at least. I chalk this up to I just don't understand enough about numpy, yet. Any help would be appreciated. Here is the test code: import numpy as np import time import random # Create a simple loop and count up the number of matching cases # Basic cross tabulation or histogram of the data # This approach is prefered because of the need to customize the calculation potentially for each cell. def simpleLoop(c): #number of items per inner loop a_cnt = len(np.unique(c[:,0])) b_cnt = len(np.unique(c[:,1])) idx = 0 result = np.zeros(b_cnt * a_cnt) for i in np.unique(c[:,0]): for j in np.unique(c[:,1]): result[idx] = np.sum(1*(c[:,0] == i) & (c[:,1] == j)) idx += 1 result.resize(len(result)/b_cnt,b_cnt) return result # Use numpys histogram method to calculate the matrix of combinations and the number of cases in each one. def simpleHistogram(c): #number of items per inner loop return np.histogramdd((c[:,0],c[:,1]), bins=[np.unique(c[:, 0]),range(1,11)]) # Variation1 of simple histogram def simpleHistogram1(c): #number of items per inner loop results = [] for i in np.unique(c[:,1]): results.append(np.histogramdd((c[:,0][c[:,1]==i]), bins=[np.unique(c[:,0])]) or 0) return np.column_stack([result[0] for result in results]) if __name__ == '__main__': a = np.random.randint(1,900,200000) b = np.random.randint(1,10,200000) c = np.column_stack((a,b)) print '---- Simple Loop ----' start = time.time() results = simpleLoop(c) print results[0] print time.time() - start print '---- Histogram dd no looping ----' start = time.time() results = simpleHistogram(c) print results[0][0] print time.time() - start print '---- Histogram run 1 time for each item in column 1 (10 times) ----' start = time.time() results = simpleHistogram1(c) print results[0] print time.time() - start

1 0

ANN: NumPy 1.6.0
by Ralf Gommers May 16, 2011

May 16, 2011

Hi, I am pleased to announce the release of NumPy 1.6.0. This release is the result of 9 months of work, and includes many new features, performance improvements and bug fixes. Some highlights are: - Re-introduction of datetime dtype support to deal with dates in arrays. - A new 16-bit floating point type. - A new iterator, which improves performance of many functions. Sources and binaries can be found at http://sourceforge.net/projects/numpy/files/NumPy/1.6.0/<http://sourceforge.net/projects/numpy/files/NumPy/1.6.0rc2/> For release notes see below. Thank you to everyone who contributed to this release. Enjoy, The NumPy developers ========================= NumPy 1.6.0 Release Notes ========================= This release includes several new features as well as numerous bug fixes and improved documentation. It is backward compatible with the 1.5.0 release, and supports Python 2.4 - 2.7 and 3.1 - 3.2. Highlights ========== * Re-introduction of datetime dtype support to deal with dates in arrays. * A new 16-bit floating point type. * A new iterator, which improves performance of many functions. New features ============ New 16-bit floating point type ------------------------------ This release adds support for the IEEE 754-2008 binary16 format, available as the data type ``numpy.half``. Within Python, the type behaves similarly to `float` or `double`, and C extensions can add support for it with the exposed half-float API. New iterator ------------ A new iterator has been added, replacing the functionality of the existing iterator and multi-iterator with a single object and API. This iterator works well with general memory layouts different from C or Fortran contiguous, and handles both standard NumPy and customized broadcasting. The buffering, automatic data type conversion, and optional output parameters, offered by ufuncs but difficult to replicate elsewhere, are now exposed by this iterator. Legendre, Laguerre, Hermite, HermiteE polynomials in ``numpy.polynomial`` ------------------------------------------------------------------------- Extend the number of polynomials available in the polynomial package. In addition, a new ``window`` attribute has been added to the classes in order to specify the range the ``domain`` maps to. This is mostly useful for the Laguerre, Hermite, and HermiteE polynomials whose natural domains are infinite and provides a more intuitive way to get the correct mapping of values without playing unnatural tricks with the domain. Fortran assumed shape array and size function support in ``numpy.f2py`` ----------------------------------------------------------------------- F2py now supports wrapping Fortran 90 routines that use assumed shape arrays. Before such routines could be called from Python but the corresponding Fortran routines received assumed shape arrays as zero length arrays which caused unpredicted results. Thanks to Lorenz Hüdepohl for pointing out the correct way to interface routines with assumed shape arrays. In addition, f2py supports now automatic wrapping of Fortran routines that use two argument ``size`` function in dimension specifications. Other new functions ------------------- ``numpy.ravel_multi_index`` : Converts a multi-index tuple into an array of flat indices, applying boundary modes to the indices. ``numpy.einsum`` : Evaluate the Einstein summation convention. Using the Einstein summation convention, many common multi-dimensional array operations can be represented in a simple fashion. This function provides a way compute such summations. ``numpy.count_nonzero`` : Counts the number of non-zero elements in an array. ``numpy.result_type`` and ``numpy.min_scalar_type`` : These functions expose the underlying type promotion used by the ufuncs and other operations to determine the types of outputs. These improve upon the ``numpy.common_type`` and ``numpy.mintypecode`` which provide similar functionality but do not match the ufunc implementation. Changes ======= ``default error handling`` -------------------------- The default error handling has been change from ``print`` to ``warn`` for all except for ``underflow``, which remains as ``ignore``. ``numpy.distutils`` ------------------- Several new compilers are supported for building Numpy: the Portland Group Fortran compiler on OS X, the PathScale compiler suite and the 64-bit Intel C compiler on Linux. ``numpy.testing`` ----------------- The testing framework gained ``numpy.testing.assert_allclose``, which provides a more convenient way to compare floating point arrays than `assert_almost_equal`, `assert_approx_equal` and `assert_array_almost_equal`. ``C API`` --------- In addition to the APIs for the new iterator and half data type, a number of other additions have been made to the C API. The type promotion mechanism used by ufuncs is exposed via ``PyArray_PromoteTypes``, ``PyArray_ResultType``, and ``PyArray_MinScalarType``. A new enumeration ``NPY_CASTING`` has been added which controls what types of casts are permitted. This is used by the new functions ``PyArray_CanCastArrayTo`` and ``PyArray_CanCastTypeTo``. A more flexible way to handle conversion of arbitrary python objects into arrays is exposed by ``PyArray_GetArrayParamsFromObject``. Deprecated features =================== The "normed" keyword in ``numpy.histogram`` is deprecated. Its functionality will be replaced by the new "density" keyword. Removed features ================ ``numpy.fft`` ------------- The functions `refft`, `refft2`, `refftn`, `irefft`, `irefft2`, `irefftn`, which were aliases for the same functions without the 'e' in the name, were removed. ``numpy.memmap`` ---------------- The `sync()` and `close()` methods of memmap were removed. Use `flush()` and "del memmap" instead. ``numpy.lib`` ------------- The deprecated functions ``numpy.unique1d``, ``numpy.setmember1d``, ``numpy.intersect1d_nu`` and ``numpy.lib.ufunclike.log2`` were removed. ``numpy.ma`` ------------ Several deprecated items were removed from the ``numpy.ma`` module:: * ``numpy.ma.MaskedArray`` "raw_data" method * ``numpy.ma.MaskedArray`` constructor "flag" keyword * ``numpy.ma.make_mask`` "flag" keyword * ``numpy.ma.allclose`` "fill_value" keyword ``numpy.distutils`` ------------------- The ``numpy.get_numpy_include`` function was removed, use ``numpy.get_include`` instead.

3 6

Creating a 1-d structured array
by Bruce Southey May 16, 2011

May 16, 2011

Hi, How do you create a 'single' structured array using np.array()? Basically I am attempting to do something like this that does not work: a=np.array([1,2, 3,4, 5,6], dtype=np.dtype([('foo', int)])) I realize that this is essentially redundant as if A is an 1-d array then a structured array with a named field 'foo' is the same thing - A would be A['foo'], just shorter. So if that is valid then a clearer error message is required to indicate this and provide the suitable error message to address ticket 1264 (http://projects.scipy.org/numpy/ticket/1264). $ python Python 2.7 (r27:82500, Sep 16 2010, 18:02:00) [GCC 4.5.1 20100907 (Red Hat 4.5.1-3)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> import numpy as np >>> a=np.array([1,2, 3,4, 5,6], dtype=np.dtype([('foo', int)])) Traceback (most recent call last): File "<stdin>", line 1, in <module> TypeError: expected a readable buffer object >>> np.array([(1,2), (3,4), (5,6)], dtype=np.dtype([('foo', int), ('foo2', int)])) array([(1, 2), (3, 4), (5, 6)], dtype=[('foo', '<i8'), ('foo2', '<i8')]) >>> a=np.array([(1), (3), (5)], dtype=np.dtype([('foo', int)])) Traceback (most recent call last): File "<stdin>", line 1, in <module> TypeError: expected a readable buffer object >>> np.array([(1), (3), (5)], dtype=np.dtype(int)) array([1, 3, 5]) >>> np.array((1,2, 3,4, 5,6)) array([1, 2, 3, 4, 5, 6]) Thanks Bruce

4 9

Documentation errors
by Charles R Harris May 15, 2011

May 15, 2011

A number of the ufuncs in the reference doc are missing arguments, for example remainder(x1) Return element-wise remainder of division. mod(x1) Return element-wise remainder of division. fmod(x1) Return the element-wise remainder of division. absolute() Calculate the absolute value element-wise. rint() Round elements of the array to the nearest integer. sign() Returns an element-wise indication of the sign of a number. conj() Return the And this is true of binary ufuncs as well greater(x1) Return the truth value of (x1 > x2) element-wise. greater_equal(x1) Return the truth value of (x1 >= x2) element-wise. less(x1) Return the truth value of (x1 < x2) element-wise. less_equal(x1) Return the truth value of (x1 =< x2) element-wise. not_equal(x1) Return (x1 != x2) element-wise. equal(x1) Return (x1 == x2) element-wise. These are generated by the .. autosummary:: command, so the error probably lies there. Also, what is the easiest way to generate the reference guide for testing purposes? Can the autosummary command deal with classes that are only generated when modules are imported? I want to document the new polynomial classes but they are generated on the fly from templates. Chuck

3 6

RHE 4 and test failure with complex numbers (tickets: 1323, 1324 and 1325)
by Bruce Southey May 15, 2011

May 15, 2011

Hi, I think tickets 1323, 1324 and 1325 have a common source of error and, thus, are duplicates. http://projects.scipy.org/numpy/ticket/1323 http://projects.scipy.org/numpy/ticket/1324 http://projects.scipy.org/numpy/ticket/1325 These tickets could be closed just because Red Hat Enterprise Linux (RHE) 4 has been replaced by RHE 5 (first released 2007-03-14), which in turn is being replaced by RHE 6 (2010-11-10). I could not replicate this with numpy version 1.5.0 and 1.5.1 with Python 2.5.5 build with GCC 4.5.1 20100924 (Red Hat 4.5.1-4). As indicated by ticket 1325, this might be the gcc version used since RHE 4 used GCC 3.4 (RHE 5 is GCC 4.1 and RHE 6 is GCC 4.4 - http://www.redhat.com/rhel/compare/). The compiler versions used to build Python and numpy are needed unless someone can replicate these failures using either RHE 4 or older gcc version with the current numpy release. Bruce

2 1

AUTO: Jon K Peck is out of the office (returning 05/18/2011)
by Jon K Peck May 12, 2011

May 12, 2011

I am out of the office until 05/18/2011. I am out of the office traveling Wed - Thursday, May 11-12 and Saturday-Tuesday, May 14-17. I will have limited access to email during this time, so I will be delayed in responding. Note: This is an automated response to your message "NumPy-Discussion Digest, Vol 56, Issue 25" sent on 5/12/11 11:00:02. This is the only notification you will receive while this person is away.

1 0

ufunc 's order of execution [relevant when output overlaps with input]
by srean May 12, 2011

May 12, 2011

Hi, is there a guarantee that ufuncs will execute left to right and in sequential order ? For instance is the following code standards compliant ? >>> import numpy as n >>> a=n.arange(0,5) array([0, 1, 2, 3, 4]) >>> n.add(a[0:-1], a[1:], a[0:-1]) array([1, 3, 5, 7]) The idea was to reuse and hence save space. The place where I write to is not accessed again. I am quite surprised that the following works correctly. >>>n.add.accumulate(a,out=a) I guess it uses a buffer and rebinds `a` to that buffer at the end. It is also faster than >>> n.add(a[0:-1], a[1:], a[1:]) --sean

2 2