week-8 More ufunc overrides

Comments

This week was spent implementing the API for overriding ufuncs. After a lengthy discussion, I decided that it would be best to proceed with implementing things using a duck typing mechanism, similar to the one Python uses, instead of multi-methods. Nathaniel Smith and Pauli Virtanen had a lot of great suggestions on the interface which has been totally changed since last week. __ufunc_override__ has become __numpy_ufunc__ and is now a method which takes arguments like:

def __numpy_ufunc__(self, ufunc, method, i, inputs, kwargs):
    ...

Where:

  • ufunc is the ufunc object that was called.
  • method is a string indicating which Ufunc method was called (one of "__call__", "reduce", "reduceat", "accumulate", "outer", "inner").
  • i is the index of self in inputs.
  • inputs is a tuple of the input arguments to the ufunc
  • kwargs is a dictionary containing the optional input arguments of the ufunc. The out argument is always contained in kwargs, if given.

The implementation is best described by the NEP:

  • If one of the input arguments implements __numpy_ufunc__ it is executed instead of the Ufunc.

  • If more than one of the input arguments implements __numpy_ufunc__, they are tried in the following order: subclasses before superclasses, otherwise left to right. The first __numpy_ufunc__ method returning something else than NotImplemented determines the return value of the Ufunc. This is to accommodate Python's own MRO behavior.

  • If all __numpy_ufunc__ methods of the input arguments return NotImplemented, a TypeError is raised.

  • If a __numpy_ufunc__ method raises an error, the error is propagated immediately.

Demo

There have been several models working in a limited fashion so far, but here is an example. I added this to scipy/sparse/base.py

def __numpy_ufunc__(*args, **kwargs): self = args[0] ufunc = args[1] method = args[2] i = args[3] ufunc_args = args[4] if ufunc == np.multiply and method == "__call__": other_args = tuple(j for j in ufunc_args if j is not ufunc_args[i]) return self.multiply(*other_args, **kwds) else: return NotImplemented

So __numpy_ufunc__ is now a method on all sparse matrices.

In [21]:
import scipy.sparse as sp
import numpy as np

a = np.random.randint(5, size=(3,3))
b = np.random.randint(5, size=(3,3))

asp = sp.csr_matrix(a)
bsp = sp.csr_matrix(b)
In [22]:
np.multiply(a,b)
Out[22]:
array([[ 2,  0, 12],
       [ 0,  0, 12],
       [ 3,  0,  9]])
In [23]:
np.multiply(asp, bsp).todense()
Out[23]:
matrix([[ 2,  0, 12],
        [ 0,  0, 12],
        [ 3,  0,  9]], dtype=int64)
In [24]:
np.multiply(asp,b)
Out[24]:
matrix([[ 2,  0, 12],
        [ 0,  0, 12],
        [ 3,  0,  9]], dtype=int64)
In [25]:
np.multiply(a, bsp)
Out[25]:
matrix([[ 2,  0, 12],
        [ 0,  0, 12],
        [ 3,  0,  9]], dtype=int64)

week-7 Adding Ufunc Overrides to NumPy

Comments

This week was spent adding the Ufunc override functionality to NumPy. There is a pull request, a lot of review still needs to be done. The current implementation is still subject to change, and docs and tests need to be added. But I feel much more comfortable with ufuncs inner workings and the Python/NumPy C API now.

Implementation Changes

The code path of interest here is in numpy/core/src/umath/ufunc_object.c when ufunc_generic_call is called. Here there is a new function which checks the arguments for an overriding function. If it is found, it is called with the arguments.

A few thing which could be changed in the implementation:

  • The name of the attribute. Since numpy.dot is not a ufunc but we still want to override it. __ufunc_override__ might not be the best name.

  • The overriding function should be moved into the array API or ufunc API.

  • The keying of __ufunc_override__. Currently this dictionary is keyed with the __name__ of the ufnucs, however this probably is not unique. I'm currently changing the implementation to be keyed with the callable ufunc function, I'll push this if it is working out alright because it seems more reliable.

  • Coercing the args before passing them to the overriding function. Currently args are passed unchanged to the overriding function, however this can be a problem if the overriding function is a method that expects self to be the first argument. A possible solution is to just change the args to have the overriding arg come first. But this might be problematic for ufuncs which are not associative.

  • ??? I don't have a lot of feedback from the NumPy community yet. There could be a bigger problem I'm missing.

New Ufunc Overhead

Now _find_ufunc_override is called for every ufunc call. For standard cases where the ufunc is operating on ndarrays or scalars the overhead comes from calling:

if (PyArray_CheckExact(obj) || PyArray_IsAnyScalar(obj)) {
    continue;}

On each argument. This adds about 5μs to 10μs, or 4% to 5% decrease in speed. (Tested with %timeit and multiply and add). Either way this seems to be acceptably small. I'm going to try Arink Verma's setup to get a more accurate benchmarks.

week-6 Ufunc Overrides

Comments

This week marked the beginning of the next stage of my proposal ambiguously named "NumPy Interactions". The first item scheduled is "Modify NumPy universal functions (ufuncs) to aware of sparse matrices (1 week)" The first thing I realized this wee is that my time allotment was too short. This is going to take a several weeks.

I'll explain the problem sparse matrices have with sparse matrices first, then describe my proposed solution.

Problem

Ufunc don't play well with objects that are not derived from ndarrays. There is some limited functionality such as __array_prepare__ and __array_wrap__ but these methods still prohibit you from accessing the array data, or changing it shape (not very useful for coercing sparse matrices into something ufuncable)

Sparse matrices cannot be used with ufuncs reliably, for example:

In [1]:
import numpy as np
import scipy.sparse as sp

a = np.random.randint(5, size=(3,3))
b = np.random.randint(5, size=(3,3))

asp = sp.csr_matrix(a)
bsp = sp.csr_matrix(b)
In [2]:
a, b
Out[2]:
(array([[0, 4, 4],
       [1, 3, 2],
       [1, 3, 1]]),
 array([[0, 1, 0],
       [0, 0, 1],
       [4, 0, 1]]))
In [3]:
np.multiply(a, b) # The right answer
Out[3]:
array([[0, 4, 0],
       [0, 0, 2],
       [4, 0, 1]])
In [4]:
np.multiply(asp, bsp).todense() # calls __mul__ which does matrix multiplication
Out[4]:
matrix([[16,  0,  8],
        [ 8,  1,  5],
        [ 4,  1,  4]], dtype=int64)
In [5]:
np.multiply(a, bsp) # Returns NotImplemented to user, bad!
Out[5]:
NotImplemented

Returning NotImplemented to user should not happen. I'm not sure if the blame for this lies in scipy.sparse or numpy, but it should be fixed.

In [6]:
np.multiply(asp, b)

/home/blake/.virtualenvs/scipy3/lib/python3.2/site-packages/scipy/sparse/compressed.py:350: SparseEfficiencyWarning: Comparing sparse matrices using >= and <= is inefficient, try using <, >, or !=, instead.
  "try using <, >, or !=, instead.", SparseEfficiencyWarning)
Out[6]:
array([[ <3x3 sparse matrix of type '<class 'numpy.int64'>'
	with 8 stored elements in Compressed Sparse Row format>,
        <3x3 sparse matrix of type '<class 'numpy.int64'>'
	with 8 stored elements in Compressed Sparse Row format>,
        <3x3 sparse matrix of type '<class 'numpy.int64'>'
	with 8 stored elements in Compressed Sparse Row format>],
       [ <3x3 sparse matrix of type '<class 'numpy.int64'>'
	with 8 stored elements in Compressed Sparse Row format>,
        <3x3 sparse matrix of type '<class 'numpy.int64'>'
	with 8 stored elements in Compressed Sparse Row format>,
        <3x3 sparse matrix of type '<class 'numpy.int64'>'
	with 8 stored elements in Compressed Sparse Row format>],
       [ <3x3 sparse matrix of type '<class 'numpy.int64'>'
	with 8 stored elements in Compressed Sparse Row format>,
        <3x3 sparse matrix of type '<class 'numpy.int64'>'
	with 8 stored elements in Compressed Sparse Row format>,
        <3x3 sparse matrix of type '<class 'numpy.int64'>'
	with 8 stored elements in Compressed Sparse Row format>]], dtype=object)

I'm not sure what happened here either, but I think raising TypeError would be preferable.

Proposed Solution

Basically, I want to implement a mechanism that will let arguments override ufuncs. So any object that ufuncs do not work with can define the attribute __ufunc_override__ as a dictionary keyed the name of the ufunc to be overridden, and valued with the callable function which should override it.

Each time a ufunc is run it will check each argument for the __ufunc_override__ attribute. So this will add the overhead of one attribute check for each argument. The logic of calling the correct replacement function will come after this check. So it won't impact performance for normal operations that don't want to override the ufunc.

I prototyped this behavior in python. ufunc_override.py (link) defines a decorator that preforms the logic described above. Then all the ufuncs in numpy are decorated with it. So import ufunc_override imports a wrapped version of numpy.

I'm also working with a slightly modified version of scipy.sparse. With this added to the base class.

__ufunc_override__ = {'multiply':multiply}

The relavant code on this branch is here.

A little demo.

In [7]:
import ufunc_override as uo
In [8]:
asp.__ufunc_override__
Out[8]:
{'multiply': <function scipy.sparse.base.multiply>}
In [9]:
uo.multiply(asp, bsp).todense()
Out[9]:
matrix([[0, 4, 0],
        [0, 0, 2],
        [4, 0, 1]], dtype=int64)
In [10]:
uo.multiply(a, bsp)
Out[10]:
matrix([[0, 4, 0],
        [0, 0, 2],
        [4, 0, 1]], dtype=int64)
In [11]:
uo.multiply(asp, b)
Out[11]:
matrix([[0, 4, 0],
        [0, 0, 2],
        [4, 0, 1]], dtype=int64)

It works! How about we try with something a little more basic too, a bare bones class that should override the add ufunc.

In [12]:
class DemoClass(object):
    def replacement_func(*args):
        return 42 # The answer.
    __array_priority__= 66 # Some number that's > 10.1
    __ufunc_override__={'add':replacement_func}
    
foo = DemoClass()
In [13]:
uo.add(a, foo)
Out[13]:
42
In [14]:
uo.add(asp, foo)
Out[14]:
42

backporting

Comments

I recently sent a pull request for adding sparse matrix boolean comparisons to SciPy. However when comparing a sparse matrix with dense matrix, with the dense matrix on the LHS, the behavior was unpredictable. This is because the __array_priority__ was not being respected. Which was fixed in NumPy by this pull request #3324 (Thanks nouiz!). This will be included in NumPy 1.8, but I would like to have it in the current stable version of NumPy, so that I can use it in my work with SciPy. So I am backporting it.

I thought this backporting process was a insightful demonstration of slightly more advanced git topics. So I'm writing this walkthrough.

Basically we make a branch off of numpy/maintenance/1.7.x, cherry pick a commit from numpy/master then submit a pull request back to numpy/maintenance/1.7.x.

  1. Assuming you already have a fork of Numpy on github. We need to update it from upstream.

    # Add upstream.
git remote add upstream https://github.com/numpy/numpy.git

# Get the latest updates.
git fetch upstream

# Make sure you are on master.
git checkout master

# Apply the updates locally.
git rebase upstream/master

# Push the updated code to your github repo.
git push origin

  2. Next we need to make the branch we will work on. This needs to be based on the older version of numpy (not master).

    # Make a new branch based of numpy/maintenance/1.7.x, 
# backport-3324 is our new name for the branch.
git checkout -b backport-3324 upstream/maintenance/1.7.x

  3. Now we need to apply the changes from master to this branch using cherry-pick

    # This pull request included commits aa7a047 to c098283 (inclusive)
# so use the .. syntax, the ^ makes the range inclusive.
git cherry-pick aa7a047^..c098283

...

# Fix any conflicts then if needed.
git cherry-pick --continue

  4. I ran into some conflicts cherry picking here. These are resolved the same as with merge/rebase conflicts. Except here you can use git blame to see the difference between master and the backported branch to make sure nothing gets screwed up.

  5. Push your new branch to your git hub repo.

    git push -u origin backport-3324

Finally make your pull request using github.com. My pull request is here.

I wonder how long it will take for Travis ci to update it version of NumPy for testing SciPy?