Mailman 3 January 2013 - NumPy-Discussion

natural alignment
by Henry Gomersall Jan. 9, 2013

Jan. 9, 2013

Further to my previous emails about getting SIMD aligned arrays, I've noticed that numpy arrays aren't always naturally aligned either. For example, numpy.float96 arrays are not always aligned on 12-byte boundaries under 32-bit linux/gcc. Indeed, .alignment on the array always seems to return 4 (with 64-bit, .alignment returns 4, 8, and 16 for float32, float64 and longdouble respectively). Can I assume _anything_ in general about the alignment of a numpy array? (I mean, based on what all implementations of the underlying malloc etc will return). Should I rely on what is returned from .alignment? cheers, Henry

1 0

Bug with ufuncs made with frompyfunc
by OKB (not okblacke) Jan. 9, 2013

Jan. 9, 2013

A bug causing errors with using methods of ufuncs created with frompyfunc was mentioned on the list over a year ago: http://mail.scipy.org/pipermail/numpy-discussion/2011- September/058501.html Is there any word on the status of this bug? I wasn't able to find a ticket in the bug tracker.

2 1

Embedded NumPy LAPACK errors
by mw＠eml.cc Jan. 8, 2013

Jan. 8, 2013

Hiall, I am trying to embed numerical code in a mexFunction, as called by MATLAB, written as a Cython function. NumPy core functions and BLAS work fine, but calls to LAPACK function such as SVD seem to be made against to MATLAB's linked MKL, and this generates MKL errors. When I try this with Octave, it works fine, presumably because it is compiled against the same LAPACK as the NumPy I am embedding. Assuming I haven't made big mistakes up to here, I have the following questions: Is there a way to request numpy.linalg to use a particular LAPACK library, e.g. /usr/lib/liblapack.so ? If not, is there a reasonable way to build numpy.linalg such that it interfaces with MKL correctly ? thanks in advance for any help, Marmaduke [1] The Cython code in question : https://gist.github.com/4433635 Please see the mexFunction at the bottom.

3 3

Re: [Numpy-discussion] Which Python to use for Mac binaries
by Robert Lupton the Good Jan. 7, 2013

Jan. 7, 2013

I am sympathetic with this attitude ("Avoid using system Python for anything"), but I don't think it's the right one. For example, the project I'm working on (HSC/LSST for astrofolk) is using python/C++ for astronomical imaging, and we expect to have the code running on a significant number of end-user laptops. If the instructions start out with: 0. Install a new version of python it's a significant barrier. What if they've already involved other packages into the system python? Python is a central part of modern operating systems, and people should not have to manage two versions of python to use numpy. It's tempting to say, "First install g++ 4.7 so we can use C++11 features" it's simply not viable, and I think that saying, "first install a new python" is comparable. Yes, I know that you can have more than one python/compiler suite installed simultaneously, but that's not something for casual users to have to get involved in. R > On Sun, Jan 6, 2013 at 2:04 AM, Ralf Gommers <ralf.gommers(a)gmail.com> wrote: >>> Which exact Python do we need to use on Mac? Do we need to use the >>> binary installer from python.org? ... > >> Avoid using system Python for anything. The first thing to do on any new OS >> X system is install Python some other way, preferably from python.org. > > +1

2 1

ANN: HDF5 for Python (h5py) 2.1.1
by Andrew Collette Jan. 7, 2013

Jan. 7, 2013

Announcing HDF5 for Python (h5py) 2.1.1 ======================================= HDF5 for Python 2.1.1 is now available! This bugfix release also marks a number of changes for the h5py project intended to make the development process more responsive, including a move to GitHub and a switch to a rapid release model. Development has moved over to GitHub at http://github.com/h5py/h5py. We welcome bug reports and pull requests from anyone interested in contributing. Releases will now be made every 4-6 weeks, in order to get bugfixes and new features out to users quickly while still leaving time for testing. * New main website: http://www.h5py.org * Mailing list: http://groups.google.com/group/h5py What is h5py? ============= The h5py package is a Pythonic interface to the HDF5 binary data format. It lets you store huge amounts of numerical data, and easily manipulate that data from NumPy. For example, you can slice into multi-terabyte datasets stored on disk, as if they were real NumPy arrays. Thousands of datasets can be stored in a single file, categorized and tagged however you want. H5py uses straightforward NumPy and Python metaphors, like dictionary and NumPy array syntax. For example, you can iterate over datasets in a file, or check out the .shape or .dtype attributes of datasets. You don't need to know anything special about HDF5 to get started. In addition to the easy-to-use high level interface, h5py rests on a object-oriented Cython wrapping of the HDF5 C API. Almost anything you can do from C in HDF5, you can do from h5py. Best of all, the files you create are in a widely-used standard binary format, which you can exchange with other people, including those who use programs like IDL and MATLAB. What's new in 2.1.1? ==================== This is a bugfix release. The most substantial changes were: * Fixed a memory leak related to variable-length strings (Thanks to Luke Campbell for extensive testing and bug reports) * Fixed a threading deadlock related to the use of H5Aiterate * Fixed a double INCREF memory leak affecting Unicode variable-length strings * Fixed an exception when taking the repr() of objects with non-ASCII names

1 0

high dimensional array -> python scalar/index
by Sebastian Berg Jan. 6, 2013

Jan. 6, 2013

Question for everyone, is this really reasonable: >>> import numpy as np >>> from operator import index >>> index(np.array([[5]])) 5 >>> int(np.array([[5]])) 5 >>> [0,1,2,3][np.array([[2]])] 2 To me, this does not make sense, why should we allow to use a high dimensional object like a normal scalar (its ok for 0-d arrays I guess)? Personally I would be for deprecating these usages, even if that (probably) means you cannot reshape your array with a matrix (as it is 2D) ;-): >>> np.arange(10).reshape(np.matrix([5,-1]).T) array([[0, 1], [2, 3], [4, 5], [6, 7], [8, 9]])

2 2

Rank-0 arrays - reprise
by Matthew Brett Jan. 6, 2013

Jan. 6, 2013

Hi, Following on from Nathaniel's explorations of the scalar - array casting rules, some resources on rank-0 arrays. The discussion that Nathaniel tracked down on "rank-0 arrays"; it also makes reference to casting. The rank-0 arrays seem to have been one way of solving the problem of maintaining array dtypes other than bool / float / int: http://mail.scipy.org/pipermail/numpy-discussion/2002-September/001612.html Quoting from an email from Travis in that thread, replying to an email from Tim Hochberg: http://mail.scipy.org/pipermail/numpy-discussion/2002-September/001647.html <quote> > Frankly, I have no idea what the implimentation details would be, but > could we get rid of rank-0 arrays altogether? I have always simply found > them strange and confusing... What are they really neccesary for > (besides holding scalar values of different precision that standard > Pyton scalars)? With new coercion rules this becomes a possibility. Arguments against it are that special rank-0 arrays behave as more consistent numbers with the rest of Numeric than Python scalars. In other words they have a length and a shape and one can right N-dimensional code that works the same even when the result is a scalar. Another advantage of having a Numeric scalar is that we can control the behavior of floating point operations better. e.g. if only Python scalars were available and sum(a) returned 0, then 1 / sum(a) would behave as Python behaves (always raises error). while with our own scalars 1 / sum(a) could potentially behave however the user wanted. </quote> There seemed then to be some impetus to remove rank-0 arrays and replace them with Python scalar types with the various numpy precisions : http://mail.scipy.org/pipermail/numpy-discussion/2002-September/013983.html Travis' recent email hints at something that seems similar, but I don't understand what he means: http://mail.scipy.org/pipermail/numpy-discussion/2012-December/064795.html <quote> Don't create array-scalars. Instead, make the data-type object a meta-type object whose instances are the items returned from NumPy arrays. There is no need for a separate array-scalar object and in fact it's confusing to the type-system. I understand that now. I did not understand that 5 years ago. </quote> Travis - can you expand? I remember rank-0 arrays being confusing in that I sometimes get a python scalar and sometimes a numpy scalar, and I may want a python scalar, and have to special-case the rank-0 array, but I don't remember precisely why I needed the python scalar. Any other comments / records of rank-0 arrays being confusing? Best, Matthew

6 13

Re: [Numpy-discussion] Invalid value encoutered : how to, prevent numpy.where to do this?
by Eric Emsellem Jan. 5, 2013

Jan. 5, 2013

Thanks! This makes sense of course. And yes the operation I am trying to do is rather complicated so I need to rely on a prior selection. Now I would need to optimise this for large arrays and the code does go through these command line many many times. When I have to operate on the two different parts of the array, I guess just using the following is the fastest way (as you indicated) : result = np.empty_like(data) mask = (data == 0) result[mask] = 0.0 result[~mask] = 1.0/data[~mask] But if I only need to do this on one side of the selection, I guess I would just do: result = np.empty_like(data) mask = (data != 0) result[mask] += 1.0 / data[mask] I have tried using three version of "mask = " with the rest of the code being the same: 1- mask = where(data != 0) 2- mask = np.where(data != 0) 3- mask = (data != 0) and it looks like #3 is the fastest, then #2 (20% slower) then #1 (50% slower than #3). I am not sure why, but Is that making sense? Or is there even a faster way (for large data arrays, and complicated operations)? thanks Eric > If your operation doesn't factor like this though then you can always > use something more cumbersome like > result = np.empty_like(data) > mask = (data == 0) > result[mask] = 0 > result[~mask] = 1.0/data[~mask] > > Or in 1.7 this could be written > result = np.zeros_like(data) > np.divide(1.0, data, where=(data != 0), out=result) > > -n >

2 1

Re: [Numpy-discussion] Invalid value encoutered : how to, prevent numpy.where to do this?
by Eric Emsellem Jan. 5, 2013

Jan. 5, 2013

Thanks! This makes sense of course. And yes the operation I am trying to do is rather complicated so I need to rely on a prior selection. Now I would need to optimise this for large arrays and the code does go through these command line many many times. When I have to operate on the two different parts of the array, I guess just using the following is the fastest way (as you indicated) : result = np.empty_like(data) mask = (data == 0) result[mask] = 0.0 result[~mask] = 1.0/data[~mask] But if I only need to do this on one side of the selection, I guess I would just do: result = np.empty_like(data) mask = (data != 0) result[mask] += 1.0 / data[mask] I have tried using three version of "mask = " with the rest of the code being the same: 1- mask = where(data != 0) 2- mask = np.where(data != 0) 3- mask = (data != 0) and it looks like #3 is the fastest, then #2 (20% slower) then #1 (50% slower than #3). I am not sure why, but Is that making sense? Or is there even a faster way (for large data arrays, and complicated operations)? thanks Eric > If your operation doesn't factor like this though then you can always > use something more cumbersome like > result = np.empty_like(data) > mask = (data == 0) > result[mask] = 0 > result[~mask] = 1.0/data[~mask] > > Or in 1.7 this could be written > result = np.zeros_like(data) > np.divide(1.0, data, where=(data != 0), out=result) > > -n >

1 0

Re: [Numpy-discussion] Invalid value encoutered : how to, prevent numpy.where to do this?
by Eric Emsellem Jan. 5, 2013

Jan. 5, 2013

Thanks! This makes sense of course. And yes the operation I am trying to do is rather complicated so I need to rely on a prior selection. Now I would need to optimise this for large arrays and the code does go through these command line many many times. When I have to operate on the two different parts of the array, I guess just using the following is the fastest way (as you indicated) : result = np.empty_like(data) mask = (data == 0) result[mask] = 0.0 result[~mask] = 1.0/data[~mask] But if I only need to do this on one side of the selection, I guess I would just do: result = np.empty_like(data) mask = (data != 0) result[mask] += 1.0 / data[mask] I have tried using three version of "mask = " with the rest of the code being the same: 1- mask = where(data != 0) 2- mask = np.where(data != 0) 3- mask = (data != 0) and it looks like #3 is the fastest, then #2 (20% slower) then #1 (50% slower than #3). I am not sure why, but Is that making sense? Or is there even a faster way (for large data arrays, and complicated operations)? thanks Eric > If your operation doesn't factor like this though then you can always > use something more cumbersome like > result = np.empty_like(data) > mask = (data == 0) > result[mask] = 0 > result[~mask] = 1.0/data[~mask] > > Or in 1.7 this could be written > result = np.zeros_like(data) > np.divide(1.0, data, where=(data != 0), out=result) > > -n >

1 0

Howto bisect old commits correctly
by Sebastian Berg Jan. 5, 2013

Jan. 5, 2013

Hey, this is probably just because I do not have any experience with bisect and the like, but when I try running a bisect keep running into: ImportError: /home/sebastian/.../lib/python2.7/site-packages/numpy/core/multiarray.so: undefined symbol: PyDataMem_NEW or: RuntimeError: module compiled against API version 8 but this version of numpy is 7 I am sure I am missing something simple, but I have no idea where to look. Am I just forgetting to delete some things and my version is not clean!? Regards, Sebastian

2 2

ANN: Spyder v2.1.13
by Pierre Raybaut Jan. 5, 2013

Jan. 5, 2013

Hi all, On the behalf of Spyder's development team (http://code.google.com/p/spyderlib/people/list), I'm pleased to announce that Spyder v2.1.13 has been released and is available for Windows XP/Vista/7, GNU/Linux and MacOS X: http://code.google.com/p/spyderlib/ This is a pure maintenance release -- a lot of bugs were fixed since v2.1.11 (v2.1.12 was released exclusively inside WinPython distribution): http://code.google.com/p/spyderlib/wiki/ChangeLog Spyder is a free, open-source (MIT license) interactive development environment for the Python language with advanced editing, interactive testing, debugging and introspection features. Originally designed to provide MATLAB-like features (integrated help, interactive console, variable explorer with GUI-based editors for dictionaries, NumPy arrays, ...), it is strongly oriented towards scientific computing and software development. Thanks to the `spyderlib` library, Spyder also provides powerful ready-to-use widgets: embedded Python console (example: http://packages.python.org/guiqwt/_images/sift3.png), NumPy array editor (example: http://packages.python.org/guiqwt/_images/sift2.png), dictionary editor, source code editor, etc. Description of key features with tasty screenshots can be found at: http://code.google.com/p/spyderlib/wiki/Features On Windows platforms, Spyder is also available as a stand-alone executable (don't forget to disable UAC on Vista/7). This all-in-one portable version is still experimental (for example, it does not embed sphinx -- meaning no rich text mode for the object inspector) but it should provide a working version of Spyder for Windows platforms without having to install anything else (except Python 2.x itself, of course). Don't forget to follow Spyder updates/news: * on the project website: http://code.google.com/p/spyderlib/ * and on our official blog: http://spyder-ide.blogspot.com/ Last, but not least, we welcome any contribution that helps making Spyder an efficient scientific development/computing environment. Join us to help creating your favourite environment! (http://code.google.com/p/spyderlib/wiki/NoteForContributors) Enjoy! -Pierre

1 0

Numpy speed ups to simple tasks - final findings and suggestions
by Raul Cota Jan. 5, 2013

Jan. 5, 2013

Hello, On Dec/2/2012 I sent an email about some meaningful speed problems I was facing when porting our core program from Numeric (Python 2.2) to Numpy (Python 2.6). Some of our tests went from 30 seconds to 90 seconds for example. I saw interest from some people in this list and I left the topic saying I would do a more complete profile of the program and report back anything meaningful. It took me quite a bit to get through things because I ended up having to figure out how to create a Visual Studio project that I could debug and compile from the IDE. First, the obvious, Everything that relies heavily on Numpy for speed (mid to large arrays) is pretty much the same speed when compared to Numeric. The areas that are considerably slower in Numpy Vs Numeric are the trivial tasks that we end up using either for convenience (small arrays) or because scalar types such as 'float64' propagate everywhere throughout the program and creep into several of our data structures. This email is really only relevant to people stuck with doing trivial operations with Numpy and want a meaningful speed boost. I focused on float64. * NOTE: I ended up doing everything in Numpy 1.6.2 as opposed to using the latest stuff. I am going to guess all my findings still apply but I will only be able to confirm until later. ========================================================= In this email I include, 1) Main bottlenecks I found which I list and refer to as (a), (b) and (c). 2) The benchmark tests I made and their speed ups 3) Details on the actual changes to the C code ========================================================= Summary of conclusions, - Our code is finally running as fast as it used to by doing some changes in Numpy and also some minor changes in our code. Half of our problems were caused by instantiating small arrays several times which is fairly slow in Numpy. The other half of our problems were are caused by the slow math performance of Numpy scalars. We did find a particular python function in our code that was a big candidate to be rewritten in C and just got it done. - In Numpy I did four sets of changes in the source code. I believe three of them are relevant to every one using Numpy and one of them is probably not going to be very popular. - The main speed up is in float64 scalar operations and creation of small arrays from lists or tuples. The speed up in small array operations is only marginal but I believe there is potential to get them at least twice as fast. ========================================================= 1) By profiling the program I found three generic types of bottlenecks in Numpy that were affecting our code, a) Operations that result in Python internally raising an error e.g. PyObject_GetAttrString(obj, "__array_priority__") when __array_priority__ is not an attribute of obj b) Creation / destruction of scalar array types . In some places this was happening unnecessarily . c) Ufuncs trying to figure out the proper type for an operation (e.g. if I multiply a float64 array by a float64 array, a fair amount of time is spent deciding that it should use float64) I came up with specific changes to address (a) and (b) . I gave up on (c) for now since I couldn't think of a way to speed it up without a large re-write and I really don't know the Numpy code (never saw it before this work). ========================================================= 2) The tests I did were (some are python natives for reference), 1) Array * Array 2) PyFloat * Array 3) Float64 * Array 4) PyFloat + Array 5) Float64 + Array 6) PyFloat * PyFloat 7) Float64 * Float64 8) PyFloat * Float64 9) PyFloat * vector1[1] 10) PyFloat + Float64 11) PyFloat < Float64 12) if PyFloat < Float64: 13) Create array from list 14) Assign PyFloat to all 15) Assign Float64 to all 16) Float64 * Float64 * Float64 * Float64 * Float64 17) Float64 * Float64 * Float64 * Float64 * Float64 18) Float64 ** 2 19) PyFloat ** 2 where Array -> Numpy array of float64 of two elements (vector1 = array( [2.0, 3.1] )). PyFloat -> pyFloat = 3.1 Float64 -> Numpy scalar 'float64' (scalarFloat64 = vector1[1]) Create array from list -> newVec = array([0.2, 0.3], dtype="float64") Assign PyFloat to all -> vector1[:] = pyFloat Assign Float64 to all -> vector1[:] = scalarFloat64 I ran every test 100000 and timed it in seconds. These are the base timings with the original Numpy TIME[s] TEST 1) 0.2003 Array * Array 2) 0.2502 PyFloat * Array 3) 0.2689 Float64 * Array 4) 0.2469 PyFloat + Array 5) 0.2640 Float64 + Array 6) 0.0055 PyFloat * PyFloat 7) 0.0278 Float64 * Float64 8) 0.0778 PyFloat * Float64 9) 0.0893 PyFloat * vector1[1] 10) 0.0767 PyFloat + Float64 11) 0.0532 PyFloat < Float64 12) 0.0543 if PyFloat < Float64 : 13) 0.6788 Create array from list 14) 0.0708 Assign PyFloat to all 15) 0.0775 Assign Float64 to all 16) 0.2994 Float64 * pyFloat * pyFloat * pyFloat * pyFloat 17) 0.1053 Float64 * Float64 * Float64 * Float64 * Float64 18) 0.0918 Float64 ** 2 19) 0.0156 pyFloat ** 2 - Test (13) is the operation that takes the longest overall - PyFloat * Float64 is 14 times slower than PyFloat * PyFloat By addressing bottleneck (a) I got the following ratios of time (BaseTime/NewTime) i.e. RATIO > 1 means GOOD . RATIO TEST 1) 1.1 Array * Array 2) 1.1 PyFloat * Array 3) 1.1 Float64 * Array 4) 1.1 PyFloat + Array 5) 1.2 Float64 + Array 6) 1.0 PyFloat * PyFloat 7) 1.7 Float64 * Float64 8) 2.8 PyFloat * Float64 9) 2.1 PyFloat * vector1[1] 10) 2.8 PyFloat + Float64 11) 3.3 PyFloat < Float64 12) 3.3 if PyFloat < Float64: 13) 3.2 Create array from list 14) 1.2 Assign PyFloat to all 15) 1.2 Assign Float64 to all 16) 2.9 Float64 * pyFloat * pyFloat * pyFloat * pyFloat 17) 1.7 Float64 * Float64 * Float64 * Float64 * Float64 18) 2.4 Float64 ** 2 19) 1.0 pyFloat ** 2 Speed up from Test (13) and (16) resulted in a big speed boost in our code Keeping the changes above. By addressing (b) in a way that did not change the data types of the return values I got the following ratios of time (BaseTime/NewTime) i.e. RATIO > 1 means GOOD . RATIO TEST 1) 1.1 Array * Array 2) 1.1 PyFloat * Array 3) 1.2 Float64 * Array 4) 1.1 PyFloat + Array 5) 1.2 Float64 + Array 6) 1.0 PyFloat * PyFloat 7) 1.7 Float64 * Float64 8) 4.3 PyFloat * Float64 9) 3.1 PyFloat * vector1[1] 10) 4.4 PyFloat + Float64 11) 9.3 PyFloat < Float64 12) 9.2 if PyFloat < Float64 : 13) 3.2 Create array from list 14) 1.2 Assign PyFloat to all 15) 1.2 Assign Float64 to all 16) 4.7 Float64 * pyFloat * pyFloat * pyFloat * pyFloat 17) 1.8 Float64 * Float64 * Float64 * Float64 * Float64 18) 2.4 Float64 ** 2 19) 1.0 pyFloat ** 2 - Scalar operations are quite a bit faster but PyFloat * Float64 is 2.9 times slower than PyFloat * PyFloat I decided to then tackle (b) even further by changing things like PyFloat * Float64 to return a PyFloat as opposed to a Float64. This is the change that I don't think is going to be very popular. This is what I got, 1) 1.1 Array * Array 2) 1.1 PyFloat * Array 3) 1.2 Float64 * Array 4) 1.1 PyFloat + Array 5) 1.2 Float64 + Array 6) 1.0 PyFloat * PyFloat 7) 3.2 Float64 * Float64 8) 8.1 PyFloat * Float64 9) 4.1 PyFloat * vector1[1] 10) 8.3 PyFloat + Float64 11) 9.4 PyFloat < Float64 12) 9.2 if PyFloat < Float64 : 13) 3.2 Create array from list 14) 1.2 Assign PyFloat to all 15) 1.2 Assign Float64 to all 16) 17.3 Float64 * pyFloat * pyFloat * pyFloat * pyFloat 17) 3.3 Float64 * Float64 * Float64 * Float64 * Float64 18) 2.4 Float64 ** 2 19) 1.0 pyFloat ** 2 - Test (16) shows how only one Float64 spoils the speed of trivial math. Now imagine the effect in hundreds of lines like that. - Even Test (17) got faster which uses only Float64 - Test (18) Float64 ** 2 is still returning a float64 in this run. Regarding bottleneck (c) . Deciding the type of UFunc. I hacked a version for testing purposes to check the potential speed up (some dirty changes in generate_umath.py). This version avoided the overhead of the call to the calls to find the matching ufunc. The ratio of speed up for something like Array * Array was only 1.6 . This was not too exciting so I walked away for now. ========================================================= 3) These are the actual changes to the C code, For bottleneck (a) In general, - avoid calls to PyObject_GetAttrString when I know the type is List, None, Tuple, Float, Int, String or Unicode - avoid calls to PyObject_GetBuffer when I know the type is List, None or Tuple a.1) In arrayobject.h after the line #include "npy_interrupt.h" I added a couple of #define //Check for exact native types that for sure do not //support array related methods. Useful for faster checks when //validating if an object supports these methods #define ISEXACT_NATIVE_PYTYPE(op) (PyList_CheckExact(op) || (Py_None == op) || PyTuple_CheckExact(op) || PyFloat_CheckExact(op) || PyInt_CheckExact(op) || PyString_CheckExact(op) || PyUnicode_CheckExact(op)) //Check for exact native types that for sure do not //support buffer protocol. Useful for faster checks when //validating if an object supports the buffer protocol. #define NEVERSUPPORTS_BUFFER_PROTOCOL(op) ( PyList_CheckExact(op) || (Py_None == op) || PyTuple_CheckExact(op) ) a.2) In common.c above the line if ((ip=PyObject_GetAttrString(op, "__array_interface__"))!=NULL) { I added if (ISEXACT_NATIVE_PYTYPE(op)){ ip = NULL; } else{ and close the } before the line #if !defined(NPY_PY3K) In common.c above the line if (PyObject_HasAttrString(op, "__array__")) { I added if (ISEXACT_NATIVE_PYTYPE(op)){ } else{ and close the } before the line #if defined(NPY_PY3K) In common.c above the line if (PyObject_GetBuffer(op, &buffer_view, PyBUF_FORMAT|PyBUF_STRIDES I added if ( NEVERSUPPORTS_BUFFER_PROTOCOL(op) ){ } else{ and close the } before the line #endif a.3) In ctors.c above the line if ((e = PyObject_GetAttrString(s, "__array_struct__")) != NULL) { I added if (ISEXACT_NATIVE_PYTYPE(s)){ e = NULL; } else{ and close the } before the line n = PySequence_Size(s); In ctors.c above the line attr = PyObject_GetAttrString(input, "__array_struct__"); I added if (ISEXACT_NATIVE_PYTYPE(input)){ attr = NULL; return Py_NotImplemented; } else{ and close the } before the line if (!NpyCapsule_Check(attr)) { In ctors.c above the line inter = PyObject_GetAttrString(input, "__array_interface__"); I added if (ISEXACT_NATIVE_PYTYPE(input)){ inter = NULL; return Py_NotImplemented; } else{ and close the } before the line if (!PyDict_Check(inter)) { In ctors.c above the line array_meth = PyObject_GetAttrString(op, "__array__"); I added if (ISEXACT_NATIVE_PYTYPE(op)){ array_meth = NULL; return Py_NotImplemented; } else{ and close the } before the line if (context == NULL) { In ctors.c above the line if (PyObject_GetBuffer(s, &buffer_view, PyBUF_STRIDES) == 0 || I added if ( NEVERSUPPORTS_BUFFER_PROTOCOL(s) ){ } else{ and close the } before the line #endif a.4) In multiarraymodule.c above the line ret = PyObject_GetAttrString(obj, "__array_priority__"); I added if (ISEXACT_NATIVE_PYTYPE(obj)){ ret = NULL; } else{ and close the } before the line if (PyErr_Occurred()) { For bottleneck (b) b.1) I noticed that PyFloat * Float64 resulted in an unnecessary "on the fly" conversion of the PyFloat into a Float64 to extract its underlying C double value. This happened in the function _double_convert_to_ctype which comes from the pattern, _@name@_convert_to_ctype I ended up splitting _@name@_convert_to_ctype into two sections. One for double types and one for the rest of the types where I extract the C value directly if it passes the check to PyFloat_CheckExact (It could be extended for other types). in scalarmathmodule.c.src I added, /**begin repeat * #name = double# * #Name = Double# * #NAME = DOUBLE# * #PYCHECKEXACT = PyFloat_CheckExact# * #PYEXTRACTCTYPE = PyFloat_AS_DOUBLE# */ static int _@name@_convert_to_ctype(PyObject *a, npy_@name@ *arg1) { PyObject *temp; if (@PYCHECKEXACT@(a)){ *arg1 = @PYEXTRACTCTYPE@(a); return 0; } ... The rest of this function is the implementation of the original _@name@_convert_to_ctype(PyObject *a, npy_@name@ *arg1) The original implementation of _@name@_convert_to_ctype does not include double anymore, i.e. /**begin repeat * #name = byte, ubyte, short, ushort, int, uint, long, ulong, longlong, * ulonglong, half, float, longdouble, cfloat, cdouble, clongdouble# * #Name = Byte, UByte, Short, UShort, Int, UInt, Long, ULong, LongLong, * ULongLong, Half, Float, LongDouble, CFloat, CDouble, CLongDouble# * #NAME = BYTE, UBYTE, SHORT, USHORT, INT, UINT, LONG, ULONG, LONGLONG, * ULONGLONG, HALF, FLOAT, LONGDOUBLE, CFLOAT, CDOUBLE, CLONGDOUBLE# */ static int _@name@_convert_to_ctype(PyObject *a, npy_@name@ *arg1) b.2) This is the change that may not be very popular among Numpy users. I modified Float64 operations to return a Float instead of Float64. I could not think or see any ill effects and I got a fairly decent speed boost. in scalarmathmodule.c.src I modified to this, /**begin repeat * #name=(byte,ubyte,short,ushort,int,uint,long,ulong,longlong,ulonglong)*13, * (half, float, double, longdouble, cfloat, cdouble, clongdouble)*6, * (half, float, double, longdouble)*2# * #Name=(Byte,UByte,Short,UShort,Int,UInt,Long,ULong,LongLong,ULongLong)*13, * (Half, Float, Double, LongDouble, CFloat, CDouble, CLongDouble)*6, * (Half, Float, Double, LongDouble)*2# * #oper=add*10, subtract*10, multiply*10, divide*10, remainder*10, * divmod*10, floor_divide*10, lshift*10, rshift*10, and*10, * or*10, xor*10, true_divide*10, * add*7, subtract*7, multiply*7, divide*7, floor_divide*7, true_divide*7, * divmod*4, remainder*4# * #fperr=1*70,0*50,1*10, * 1*42, * 1*8# * #twoout=0*50,1*10,0*70, * 0*42, * 1*4,0*4# * #otyp=(byte,ubyte,short,ushort,int,uint,long,ulong,longlong,ulonglong)*12, * float*4, double*6, * (half, float, double, longdouble, cfloat, cdouble, clongdouble)*6, * (half, float, double, longdouble)*2# * #OName=(Byte,UByte,Short,UShort,Int,UInt,Long,ULong,LongLong,ULongLong)*12, * Float*4, Double*6, * (Half, Float, Double, LongDouble, CFloat, CDouble, CLongDouble)*6, * (Half, Float, Double, LongDouble)*2# * #OutUseName=(Byte,UByte,Short,UShort,Int,UInt,Long,ULong,LongLong,ULongLong)*12, * Float*4, out*6, * (Half, Float, out, LongDouble, CFloat, CDouble, CLongDouble)*6, * (Half, Float, out, LongDouble)*2# * #AsScalarArr=(1,1,1,1,1,1,1,1,1,1)*12, * 1*4, 0*6, * (1, 1, 0, 1, 1, 1, 1)*6, * (1, 1, 0, 1)*2# * #RetValCreate=(PyArrayScalar_New,PyArrayScalar_New,PyArrayScalar_New,PyArrayScalar_New,PyArrayScalar_New,PyArrayScalar_New,PyArrayScalar_New,PyArrayScalar_New,PyArrayScalar_New,PyArrayScalar_New)*12, * PyArrayScalar_New*4, PyFloat_FromDouble*6, * (PyArrayScalar_New, PyArrayScalar_New, PyFloat_FromDouble, PyArrayScalar_New, PyArrayScalar_New, PyArrayScalar_New, PyArrayScalar_New)*6, * (PyArrayScalar_New, PyArrayScalar_New, PyFloat_FromDouble, PyArrayScalar_New)*2# */ #if !defined(CODEGEN_SKIP_@oper@_FLAG) static PyObject * @name@_@oper@(PyObject *a, PyObject *b) { ... Same as before and ends with... #else ret = @RetValCreate@(@OutUseName@); if (ret == NULL) { return NULL; } if (@AsScalarArr@) PyArrayScalar_ASSIGN(ret, @OName@, out); #endif return ret; } #endif /**end repeat**/ I still need to do the section for when there are two return values and the power function. I am not sure what else could be there. ========================================================= That's about it. Sorry for the long email. I tried to summarize as much as possible. Let me know if you have any questions or if you want the actual files I modified. Cheers, Raul Cota

2 9

Scalar casting rules use-case reprise
by Matthew Brett Jan. 4, 2013

Jan. 4, 2013

Hi, Reading the discussion on the scalar casting rule change I realized I was hazy on the use-cases that led to the rule that scalars cast differently from arrays. My impression was that the primary use-case was for lower-precision floats. That is, when you have a large float32 arr, you do not want to double your memory use with: >>> large_float32 + 1.0 # please no float64 here Probably also: >>> large_int8 + 1 # please no int32 / int64 here. That makes sense. On the other hand these are more ambiguous: >>> large_float32 + np.float64(1) # really - you don't want float64? >>> large_int8 + np.int32(1) # ditto I wonder whether the main use-case was to deal with the automatic types of Python floats and scalars? That is, I wonder whether it would be worth considering (in the distant long term), doing fancy guess-what-you-mean stuff with Python scalars, on the basis that they are of unspecified dtype, and make 0 dimensional scalars follow the array casting rules. As in: >>> large_float32 + 1.0 # no upcast - we don't know what float type you meant for the scalar >>> large_float32 + np.float64(1) # upcast - you clearly meant the scalar to be float64 In any case, can anyone remember the original use-cases well enough to record them for future decision making? Best, Matthew

4 4

still need DLR support
by William Johnston Jan. 4, 2013

Jan. 4, 2013

Hello, I posted some time ago that I need Numpy for .NET for a C# DLR app. Has anyone made any progress on this? May I suggest this as a project? Thank you. Sincerely, William Johnston

1 0

test failures when embedded (in matlab)
by Robin Jan. 4, 2013

Jan. 4, 2013

Hi All, When using Numpy from an embedded Python (Python embedded in a Matlab mex function) I get a lot of test failures (see attached log). I am using CentOS 6.3, distribution packaged Python (2.6) and Numpy (1.4.1). Running numpy tests from a normal Python interpreter results in no errors or failures. Most of the failures look to be to do with errors calling Fortran functions - is it possible there is some linking / ABI problem? Would there be any way to overcome it? I get similar errors using EPD 7.3. Any advice appreciated, Cheers Robin

2 1

3D array problem in Python
by Happyman Jan. 2, 2013

Jan. 2, 2013

Hello I have 3 dimensional array which I want to calculate in a huge process. Everything is working well if I use ordinary way which is unsuitable in Python like the following: nums=32 rows=120 cols=150 for k in range(0,nums): for i in range(0,rows): for j in range(0,cols): if float ( R[ k ] [ i ] [ j ] ) == 0.0: val11 [ i ] =0.0 else: val11[ i ] [ j ], val22[ i ][ j ] = integrate.quad( lambda x : F1(x)*F2(x) , 0 , pi) But, this calculation takes so long time, let's say about 1 hour (theoretically)... Is there any better way to easily and fast calculate the process such as [ F( i ) for i in xlist ] or something like that rather than using for loop?

5 8

Conversion functions
by Schönberger Johannes Jan. 1, 2013

Jan. 1, 2013

Hello everyone, I recently opened a new pull request which adds the functionality to convert between degrees and degrees, minutes and seconds (https://github.com/numpy/numpy/pull/2869). The discussion is about whether such conversion functionality should be integrated into numpy at all or whether this belongs to scipy. I suggest to move the most common conversion functions (deg2rad, rad2deg, deg2dms, dms2deg and some more could be added) to a separate file `conversion.py` file in `numpy/lib`. I could implement this in a new pull request, if the general consensus is in favor of it? What are your thoughts? Regards, Johannes

3 2