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