Summer of 2013 with Numpy and Google

GSoC is a learning process - we learn to work together to solve a big problem by breaking it down into smaller chunks.

This year, 2013, I got chance to work as student contributor at Numpy. And the experience was worth paying the amount of effort. I have joined Numpy, the community that's quite friendly and helpful.

I passed final evaluation but I value experience more than consequence (off-course passed makes me happy). Experience that I had for being a part of this program. Nothing could be more inspiring and motivational than working with people from almost every part of globe connected together just because of there love and passion for code and to develop amazing things.  Also got chance to learn from such experienced and enthusiastic developers.

I have written the summary of project into two post : Profiling and Bottle-necks. All together with my mentors Charles R Harris and whole community successfully achieved what we planned. Project information and final code can be seen here [https://google-melange.appspot.com/gsoc/project/google/gsoc2013/arinkverma/79001]

What I gained this year?

The most challenging aspect is to get work with people across different geographical reasons. People from different region thinks and act differently as they might have different concern. Coming out with product which is accepted by everyone is itself a strenuous task. To work in a team with huge diversity in timezone and thinking was colorful experience.

GSoC'13 Project Summary-2 : Numpy's Bottlenecks and its Removal

In last post, I have mentioned the tools which I used for profiling numpy. Call-graph, made my life not easier but a bit simpler to detect time consuming methods. But time taken just indicate the possibility of bottlenecks as many of time consuming methods are highly optimized.  Hence, critically reading code along with call-graph is the trick

Identifying bottlenecks

Following types bottlenecks were observed

Overhead for smaller arrays

• Numpy used to release Global Interpreter Lock or GIL for all two or one operand loops. But for short length array, it used to produce relative overhead instead. 
• So, releasing GIL for smaller operations was not benefiting at all.

Redundant code doing extra work

• PyArrayCanCastArrayTo used to check evenif newtype is Null. In PyArrayFromArray, It was found that if argument newtype is Null, it get value of oldtype. So, technically it has to check if casting is possible for same type, which is useless.
• Numpy used to converts the Python scalar into its matching scalar (e.g. PyLong -> int32) and then extract the C value from the NumPy scalar.
• Clearing the error flags at every function call, then checking it. This happens unconditional even-if there is no need to do.

Improper structure

• For every single operation calls, numpy has to extract value of buffersize, errormask and name to pack and build error object. These two functions, _extract_pyvals and PyUFunc_GetPyValues together use significant time. It take useless time because all this time is spent on to look up entries in a python dict, extract them, and convert them into C level data. Not once but doing that again and again on every operation. Also which remain unused if no error occurs.
• loop selection method for scalar operation is inefficient and consume time. It check through all associated dtypes of function one by one from types array. 

Not so effective memory allocation

• When allocating the return value, numpy allocate memory twice. One for the array object itself, and a second time for the array data.
• Also, shapes + strides together have 2*ndim elements, but to hold them numpy allocate a memory region sized to hold 3*ndim elements.

GSoC'13 Project Summary-1 : Numpy's profiling

Small numpy arrays are very similar to Python scalars but numpy incurs a fair amount of extra overhead for simple operations. For large arrays this doesn't matter, but for code that manipulates a lot of small pieces of data, it can be a serious bottleneck.
For example

In [1]: x = 1.0

In [2]: numpy_x = np.asarray(x)

In [3]: timeit x + x
10000000 loops, best of 3: 61 ns per loop

In [4]: timeit numpy_x + numpy_x
1000000 loops, best of 3: 1.66 us per loop

This project involved

• profiling simple operations like the above
• determining possible bottlenecks
• devising improved algorithms to solve them, with the goal of getting the numpy time as close as possible to the Python time.

Profiling tools

The very first objective to find bottleneck is profiling for time or space. During project I have used few tools for profiling and visualizing data of numpy execution flow.

Google profiling tool

This is the suit of different tools provided by Google. It Includes TCMalloc, heap-checker, heap-profiler and cpu-profiler. As need of project was to reduce time, so CPU-Profiler was used.

Setting up Gperftools

Following are the steps used to setup python C level profiler on Ubuntu 13.04. (For any other system, options see [1])

1. Make sure to build it from source. Clone svn repository from http://gperftools.googlecode.com/svn/trunk/
2. In order to build gperftools checked out from subversion repository you need to have autoconf, automake and libtool installed.
3. First, run ./autogen.sh script which generate ./configure and other files. Then run ./configure
4. 'make check', to run any self-tests that come with the package. Check is optional but recommended to use
5. After all test gets passed, type 'sudo make install' to install the programs and any data files and documentation.

Speedup UFUNC_CHECK_STATUS by avoiding heavy clearing flag

UFUNC_CHECK_STATUS is just single macro which do both checking clearing the error flags. It clear error flags every time after checking. We should avoid clear operation if not needed, as it is a bit expensive and take significant amount of time.

The way numpy detect divide-by-zero, overflow, underflow, etc., is that before each ufunc loop it clear the FP error flags, and then after the ufunc loop we see if any have become set. And clear again. I have avoided clear if not needed to save time.

Improvement

Before each ufunc loop when PyUFunc_clearfperr() flag error is checked, then clearing them if necessary. Now, checking results in macro doesn't get ignored unlike before. Earlier time taken by PyUFunc_clearfperr() and PyUFunc_getfperr() combined was around 10%, which is now dropped to 1%, for operation which don't raise any error.

callgraph comparing performance
x = np.asarray([1]); x+x;

More

PR for this change is #3739

Improvement in _extract_pyvals

Instead of *errobj = Py_BuildValue("NO", PyBytes_FromString(name), retval); in _extract_pyvals. It alone takes 10% of time, with every operations.  Its better to avoid packing of string instead assign pointer to struct and stash them to thread local storage. Hence no allocations would be need to for errorobj every-time, which almost goes unused.

Implementation

Error object include name of ufunc and callback pointer.

typedef struct {
    PyObject_HEAD
    char *name;
    PyObject* retval;
} PyErrObject;

Fast path for consecutive same operations, it better to use caching

errvalues = PyDict_GetItem(thedict, PyUFunc_PYVALS_ERROR);    
    if (errvalues == NULL) {
        errvalues = PyObject_New(PyErrObject, &PyErrObject_Type);
        errvalues->name = name;
        errvalues->retval = retval;
        PyDict_SetItem(thedict, PyUFunc_PYVALS_ERROR, (PyObject*)errvalues);
    }
    
    errvalues->name = name;
    errvalues->retval = retval;
    *errobj = errvalues;
    Py_INCREF(errvalues);

Improvement

Callgraph PyErrObject
x = np.asarray([5.0,1.0]); x+x

Time consumption of _extract_pyvals drop to 2% from 9.3%.

More

PR for this enhancement is #3686

Scope for improvement in _extract_pyvals

I my last post, I mentioned how two functions, _extract_pyvals and PyUFunc_GetPyValues together use >12% of time. But major culprit is *errobj = Py_BuildValue("NO", PyBytes_FromString(name), retval); in  _extract_pyvals. It alone takes 10% of time, with every operations.

Improvement

Caching Py_BuildValue

First approach I take, was to cached Py_BuildValue with Thread storage or dict. With this time consumption of _extract_pyvals dropped to 4% from 12%. 

errorobj caching with PyThreadState_GetDict

But since, TLS is a bit unreliable and risky. So @juliantaylor advised that it should be avoided. Even after many fixes, commit for this didnt managed to pass all test cases.

Almost 12% of useless time consumption by PyUFunc_GetPyValues

For every single operation calls, numpy has to extract value of buffersize, errormask and name to pack and build error object. These two functions, _extract_pyvals and PyUFunc_GetPyValues together use >12% of time.

Callgraph showing contribution of PyUFunc_GetPyValues

I think it take useless time because all this time is spent on to look up entries in a python dict, extract them, and convert them into C level data. Not once but doing that again and again on every operation. Instead these values should be converted once, at time of loading or when they set in first place.

Test cases to check Integer's behvaiour

As I wrote in my last post that's, its quicker to just extract the value directly from the Python scalar. But for numpy has different scenario to handle integers based on OS 32/64 bit. For example, There are basically two standards for long on 64 bit os, Microsoft uses long = int (32 bits), Linux uses long = long long (64 bits).

Hence, before getting into speedup modifications mentioned here, there is need to have much test case to ensure behavior of integer remains same.

Test to ensure errors are raised as expected

    def test_int_raise_behaviour(self):

        def Overflow_error_func(dtype): 
            res = np.typeDict[dtype](np.iinfo(dtype).max + 1)

        for code in 'lLqQ':
            assert_raises(OverflowError, Overflow_error_func, code)

Test to check size of long as per different OS

    def test_long_os_behaviour(self):
       long_iinfo = np.iinfo('l')
        ulong_iinfo = np.iinfo('L')
        if (sys.platform == "win32" or sys.platform == "win64" or
                platform.architecture()[0] == "32bit"):
            assert_equal(long_iinfo.max, 2**31-1)
            assert_equal(long_iinfo.min, -2**31)            
            assert_equal(ulong_iinfo.max, 2**32-1)
            assert_equal(ulong_iinfo.min, 0)
        elif platform.architecture()[0] == "64bit":
            assert_equal(long_iinfo.max, 2**63-1)
            assert_equal(long_iinfo.min, -2**63)
            assert_equal(ulong_iinfo.max, 2**64-1)
            assert_equal(ulong_iinfo.min, 0)

Bottleneck in conversion of integer to NumPy Scalar

What's wrong 

For scalar operations Numpy first try to extract the underlying C value from a Python Integers. It causes bottleneck because it first converts the Python scalar into its matching NumPy scalar (e.g. PyLong -> int32) and then it extracts the C value from the NumPy scalar.

Avoiding conversion

Hence avoiding this conversation improve speed significantly.  I have avoided conversion for known integer type but extracting its value directly.

For byte, short, int, long

#if PY_VERSION_HEX >= 0x03000000
    if(PyLong_CheckExact(a)){
        *arg1 = PyLong_AsLong(a);
        return 0;        
    }
#else
    if (PyInt_CheckExact(a)){
        *arg1 = PyInt_AS_LONG(a);
        return 0;
    }
#endif

Replacement for inefficient loop selection

What's wrong

It is evident that, loop selection method for scalar operation is inefficient and consume almost 4.2% of time. It check through all associated dtypes of function one by one from types array. There is scope to make this much faster and better.

Replacing loop by specialized conditions

Most of the function share identical signature. E.g These sets (add, subtracts) , (arccos, arcsin, arctan, arcsinh, arccosh) share same signature array. As if know, there are only 32 distinct signature arrays. I make code generator to identity and make specialized condition for each distinct signature arrays. Hence, improvement of 4%.

Implementation

1. Most of functions have uniform arguments, so it will better to look them first.
2. For each distinct signature array, auto-gen lookup function having if-else condition which check and return index. E.g following code is auto-generated to quickly return innerloop index of add function

   /** Warning this file is autogenerated!!!
   
       Please make changes to the code generator program (numpy/core/code_generators/generate_umath.py)
   **/ 
   static  int type21_id3_index(int x, int y){ 
     if(x==y){ 
       if(x==NPY_HALF){ return 11;}
       if(x==NPY_TIMEDELTA){ return 19;}
       if(x>=NPY_BOOL && x<=NPY_ULONGLONG){ return x+(0);}
       if(x>=NPY_FLOAT && x<=NPY_CLONGDOUBLE){ return x+(1);}
       if(x==NPY_OBJECT){ return 21;} 
     }
     if(x==NPY_DATETIME && y==NPY_TIMEDELTA){ return 18;}
     if(x==NPY_TIMEDELTA && y==NPY_DATETIME){ return 20;} 
     return -1;
   }
   

3. Encapsulates logic, with ((PyUFuncObject *)f)->sig_index(arg1, args2)

Effect of GIL release on Numpy array operation

Release of GIL

At present numpy release Global Interpreter Lock or GIL for all two or one operand loops. Macro NPY_BEGIN_THREADS is used to save the Python state and releases the GIL. Hence it can be placed right before code that does not need the Python interpreter. Like in ufunc_object.c trivial_three_operand_loop and trivial_two_operand_loop use it for innerloop.

Not so good for small ones

But for short length array, it produces relative overhead instead. Releasing GIL for smaller operations doesn't benefit at all.
Here, Nathaniel has mentioned few things as

• Vast majority of numpy code is single-threaded, so dropping the GIL is pure overhead.
• Dropping the GIL for microseconds at a time probably produces no benefit even for multi-threaded code, since by the time the other thread gets started and starts producing useful work, the numpy loop is done.
• Most numpy code calls + a lot more than it calls sin or even ** or /.

Numpy: Improvement in PyArray_FromArray

Call-graph of Numpy scalar array addition, shows that get_ufunc_arguments contribute to almost 18% in cumulative time. 

Callgraph for x = np.asarray(1.0). 

Problem

Tracing the execution path get_ufunc_arguments under for x = numpy.asarray(1.0); x + x , which flow as PyArrayFromAny, PyArrayFromArray, PyArrayCanCastArrayTo, can_cast_scalar_to.

Finding bottleneck in Python/Numpy

Callgraph of numpy array addtion See full image

For my GSoC project, I need to profile and find bottleneck in numpy code which is mostly written in c. The quick way is to see the callgraph and identify the path which consume most of the time. I am using Google profiling tool for profiling numpy's operations.

Setting up Gperftools

Following are the steps used to setup python c level profiler on Ubuntu 13.04. For any other system, options see [1]

1. Make sure to build it from source. Clone svn repository from http://gperftools.googlecode.com/svn/trunk/
2. In order to build gperftools checked out from subversion repository you need to have autoconf, automake and libtool installed.
3. First, run ./autogen.sh script which generate ./configure and other files. Then run ./configure
4. 'make check', to run any self-tests that come with the package. Check is optional but recommended to use
5. After all test gets passed, type 'sudo make install' to install the programs and any data files and documentation.

Performance parity between numpy arrays and Python scalars

Small numpy arrays are very similar to Python scalars but numpy incurs a fair amount of extra overhead for simple operations. For large arrays this doesn't matter, but for code that manipulates a lot of small pieces of data, it can be a serious bottleneck.

For example:

 
  In [1]: x = 1.0

  In [2]: numpy_x = np.asarray(x)

  In [3]: timeit x + x
  10000000 loops, best of 3: 61 ns per loop

  In [4]: timeit numpy_x + numpy_x
  1000000 loops, best of 3: 1.66 us per loop

I tried to introduced, a short path (at present) for integer and float addition of numpy array. In umath/ufunc_type_resolution.c , ufunc lookup loop find best data types based on input operands types. In short path, rather than going to loop again for addition operation, it return the best known data type.

 
/* Short path for addition of int + int */
    int key = 0;
    for (j = 0; j < nargs; ++j) {
        key = (key<<5) + dtypes[j]->type_num;       
    }
    NPY_UF_DBG_PRINT1("key is %d\n",key);
    
    if(strcmp(ufunc_name,"add")==0){
        int rent = -1;
        if(key == 7399)
            rent = 7;
        else if(key == 12684)
            rent = 13;

        NPY_UF_DBG_PRINT1("rent is %d\n",rent);
        if(rent > 0){
            *out_innerloop = ufunc->functions[rent];
            *out_innerloopdata = ufunc->data[rent];
            NPY_UF_DBG_PRINT1("type @ hashposition %d\n",rent);
            return 0;
        }
    }   

Following are the benchmark result based on vbench for numpy