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.

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.

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)

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.

Understanding Zookeeper protocols

Credit : zookeeper.apache.org

Zookeeper is the open sourced library of cluster membership. It is a centralized service for maintaining configuration information, naming, providing distributed synchronization, and providing group services used by company like twitter.

Though Zookeeper is most used library, but it has no concrete documentation of low-level tcp, wire protocol. The code is itself is the documentation, so to know insight code has to read. Result, there are only client bindings is available for a number of languages, but very few with pure implementation in theirs native code.

After reading code of zookeeper client in Python implementation by hannosch, (as I find java a bit verbose. But for java people see  the jute file used for RPC). Also with the help of wireshark, to inspect zookeeper packets, able to understand protocols quite well.

For Scala, I have written simple client over twitter's finagle codec to demonstrate the working of zookeeper connection.