BUG: better detect memory overlap when calling gufuncs #11381

mattip · 2018-06-19T21:30:17Z

Ufuncs need to check output args for memory overlap. In PyUFunc_GeneralizedFunction, the check assumed the gufunc was operating elementwise.

Found when I started implementing matmul as a gufunc. We seem to have avoided this issue by not having an out kwarg in our current gufuncs (i.e. linalg, see issue #11380).

This may impact downstream users of gufuncs who already try to detect this condition, but brings gufuncs into alignment with ufuncs and the comment in the code "For generalized ufuncs, we can't do buffering, so must COPY or UPDATEIFCOPY".

Test added that failed before, now passes. An additional passing test added for matmul, which currently uses einsum for implementation, but that will change in the future.

mattip · 2018-06-19T21:43:13Z

Looking through scipy, astropy I find no use of PyUFunc_FromFuncAndDataAndSignature (a signature is needed to trigger this issue). The public facing Numba repo does not use it via c, and exposes it to python via _internal.fromfunc which is wrapped by build_ufunc which does not call the variant with the gufunc signature.

mattip · 2018-06-20T01:25:14Z

Found and fixed a bug in 4c75113 where ret was not precisely out in einsum if there is overlap. ret was the temporary ndarray created by writeback semantics, not the original out ndarray. If they have the same dtype the values are guaranteed to be the same, but they are not the same entity which wastes memory.

This could be a separate PR but they both come under the heading of gufuncs and memory overlap.

eric-wieser · 2018-06-21T06:49:41Z

numpy/core/src/umath/ufunc_object.c

-                      NPY_ITER_ALIGNED |
-                      NPY_ITER_OVERLAP_ASSUME_ELEMENTWISE;
+                      NPY_ITER_ALIGNED;
+        if (op[i] && PyArray_NDIM(op[i]) < 2) {


What's special about 0d / 1d arrays?

My fix was wrong. Any element-wise gufunc (one with an empty signature '(),()->()' will not hit this code path, it will not have core_enabled. Removing this flag.

mhvk

Nice 8000 catch - a copy & paste error from quite a while ago, I suspect! Looks good modulo the single-element thing you're already removing.

pv · 2018-06-21T21:10:25Z

One can also argue this is a bug in the matmul ufunc implementation, and the NPY_ITER_OVERLAP_ASSUME_ELEMENTWISE is as intended. I think there is a design decision that ufunc inner loops are supposed to be able to deal with overlapping operands/aliasing, which is certainly assumed in the scalar case. This is of course more difficult to deal with in the gufunc case (I don't recall whether the inner loop core dimensions are contiguous in memory or not).

I suspect (without trying) that matmul.accumulate and reduce don't work correctly currently, which would demonstrate it's a bug in the inner loop implementation. EDIT: I see now reduce/accumulate raise an error these operations are not defined for gufuncs, fair enough, although the reason why is not obvious at first glance.

Anyway, I'd prefer to make an explicit design decision here, rather than considering this a bugfix.

mhvk · 2018-06-21T21:26:10Z

The inner loops just get access to the arrays; the loop dimensions do not have to be contiguous.

In principle, the inner loop indeed should check whether in and out overlap.

Indeed, I just wrote code that explicitly uses the fact that in and out can be the same to speed up the loop... I hadn't quite appreciated that this change might make that much less efficient...

mattip · 2018-06-21T22:28:09Z

Closing this (I will refactor the einsum issue to a different PR), but I think that code is copy-pasted and wrong.

The documentation seems to strengthen my claim, and a couple of other comments:

We are duplicating effort. Since the flags is unconditionally set, the gufunc call tries to prevent copying, and consumers of the code who need to detect overlap will have to run that check in reverse and invent a mechanism to allocate new memory.

The C code that check the overlap states the following, which is even clearer than the documentation

       * Operands with NPY_ITER_OVERLAP_ASSUME_ELEMENTWISE enabled are not
       * considered overlapping if the arrays are exactly the same. In this
       * case, the iterator loops through them in the same order element by
       * element.  (As usual, the user-provided inner loop is assumed to be
       * able to deal with this level of simple aliasing.)

it cannot be assumed that generalized ufuncs "loop through in the same order element by element", that is the job of simple ufuncs

The only public heavy consumer of gufuncs, linalg, makes no effort to deal with overlaps, while making a large effort to deal with non-contiguous or wrongly-ordered data. On the other hand, the python wrappers for the actual _umath_linalg do not expose the out argument so the lack does not mean much
You can always add the flag to the gufunc->op_flags, but cannot remove it since the flags are only additive. This claim too is weak since there is actually no official interface to set those flags, one would have to reach inside the C struct and set them.

mattip · 2018-06-21T22:30:38Z

@mhvk can you share your gufunc? I am curious how I could turn it into a test and rewrite it to not depend on this IMO buggy behaviour.

njsmith · 2018-06-22T00:55:14Z

it cannot be assumed that generalized ufuncs "loop through in the same order element by element", that is the job of simple ufuncs

Whenever possible, I think it'd be good to move towards unifying the gufunc/ufunc semantics and code paths.

Also, I think there may be cases where gufuncs want to allow aliasing and handle it specially? For example, consider a sort gufunc with signature (n)->(n). If the input and outputs are the same, it can perform the sort in-place, which will be more efficient than allocating a temporary buffer, copying the input into it, sorting it, and then copying back into the output (= input) buffer.

It's also entirely possible that at some point we'll have simple ufuncs where handling aliasing is difficult -- e.g. think about a ufunc where the scalars are some complex multi-field dtype.

So I lean towards saying that we should handle gufuncs and ufuncs the same way here... which might include adding some flag that a loop can use to say that it doesn't handle aliasing, and ask numpy to take care of that for it using common machinery?

mhvk · 2018-06-22T01:29:29Z

@mattip - it is exactly the case of identical in and out that may need to be preserved (though your point about being able to add but not subtract a flag is taken). The PR that I mentioned is astropy/astropy#7502, but that may not help you much since it uses a rather complicated jinja2 template to help generate code for a lot of functions - but some of those functions act in-place and this is where the ability to give in=out helps. A specific example where this is checked is:

static void ufunc_loop_aper(
    char **args, npy_intp *dimensions, npy_intp* steps, void* data)
{
    npy_intp i_o;
    npy_intp n_o = *dimensions++;
    char *theta = *args++;
    npy_intp s_theta = *steps++;
    char *astrom_in = *args++;
    npy_intp s_astrom_in = *steps++;
    char *astrom = *args++;
    npy_intp s_astrom = *steps++;
    double (*_theta);
    eraASTROM (*_astrom);
    for (i_o = 0; i_o < n_o;
         i_o++, theta += s_theta, astrom += s_astrom, astrom_in += s_astrom_in) {
        _theta = ((double (*))theta);
        _astrom = ((eraASTROM (*))astrom);
        if (astrom_in != astrom) {
            memcpy(astrom, astrom_in, 1*sizeof(eraASTROM));
        }
        eraAper(*_theta, _astrom);
    }
}

mattip · 2018-06-24T23:44:15Z

there should be a point in the ufunc call that a ufunc can say "element-wise aliasing is OK", or "leave all aliasing to the inner loop" or "always copy-on-overlap"

BUG: better detect memory overlap when calling gufuncs

f93ed4f

mattip added 00 - Bug component: numpy.ufunc labels Jun 19, 2018

BUG: enforce 'ret is ouf' for 'ret = einsum(..., out=out)'

4c75113

mattip mentioned this pull request Jun 20, 2018

nditer "out" operands are not always accessible #10956

Open

eric-wieser reviewed Jun 21, 2018

View reviewed changes

mhvk reviewed Jun 21, 2018

View reviewed changes

mattip closed this Jun 24, 2018

mattip mentioned this pull request Jun 25, 2018

ENH: Break the assumption that all ufuncs and gufuncs want is element-wise loop aliasing #11416

Closed

mattip mentioned this pull request Jul 5, 2018

ENH: allocate working buffers outside ufunc's inner loop #11510

Open

This was referenced Jul 17, 2018

MAINT: refactor ufunc iter operand flags handling #11580

Merged

NEP: backwards compatibility and deprecation policy #11596

Merged

mattip deleted the gufunc-out-overlap branch October 9, 2018 17:22

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Uh oh!

BUG: better detect memory overlap when calling gufuncs #11381

BUG: better detect memory overlap when calling gufuncs #11381

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

BUG: better detect memory overlap when calling gufuncs #11381

BUG: better detect memory overlap when calling gufuncs #11381

Uh oh!

Conversation

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!