BUG: where= in gufuncs is slightly broken with out= and casts

EDIT (seberg): The original behavior change described here seems OK. However, another issue was identified: #18700 (comment)

When the out parameter on ufunc is given, numpy 1.19 (and older versions) would always pass the values of the out array to the kernel. In numpy 1.20, that changed and a fresh uninitialized array is given when the type does not match exactly.

In other words... given a no-op ufunc (kernel does nothing) that takes 1 input and 1 output. And, it is called as:

out = np.arange(10)
no_op_ufunc(inp, out=out)

Is the output guaranteed to retain the original values?

Reproducing code example:

Reproducer from numba/numba#6864

import numba
import numpy as np

@numba.guvectorize(["void(float64[:], uint8[:])"], "(n)->(n)", nopython=True)
def func(x, out):
    
    for i in range(x.size):
        # set every fourth element to 1
        if i % 4 == 0:
            out[i] = 1

x = np.random.rand(150,150)
out = np.zeros_like(x, dtype=np.int8)  # dtype does not match expected

func(x, out)

with Numba 0.53 and Numpy 1.20.1, the result is:

array([[  1,  49,  29, ...,  96,   1,   2],
       [  1,   0, -80, ..., -38,   1,  97],
       [  1,   2,   0, ...,   0,   1,   0],
       ...,
       [  1,  -1,  -1, ...,   0,   1,   0],
       [  1,   0,   0, ...,   0,   1,   0],
       [  1,   0,   0, ..., -34,   1,  97]], dtype=int8)

In np1.20, the skipped slots are containing random values.

with Numba 0.53 and Numpy 1.19.5, the result is:

array([[1, 0, 0, ..., 0, 1, 0],
       [1, 0, 0, ..., 0, 1, 0],
       [1, 0, 0, ..., 0, 1, 0],
       ...,
       [1, 0, 0, ..., 0, 1, 0],
       [1, 0, 0, ..., 0, 1, 0],
       [1, 0, 0, ..., 0, 1, 0]], dtype=int8)

In np1.19, the skipped slots are retaining the original zero values.

Error message:

No error message. The problem is a change in behavior.

NumPy/Python version information:

>>> import sys, numpy; print(numpy.__version__, sys.version)
1.20.1 3.9.2 (default, Mar  3 2021, 11:58:52)
[Clang 10.0.0 ]

I am both surprised that it changed, and that it ever "retained" the old values :(.

The only way it can work to begin with is by both filling the buffer and then writing the buffer back again, which is exactly what happens since this can't possibly retain the values (and does not):

out = np.full_like(x, 9.123, dtype=np.float16)
func(x, out)

giving:

array([[1., 9., 9., ..., 9., 1., 9.],
       [1., 9., 9., ..., 9., 1., 9.],
       [1., 9., 9., ..., 9., 1., 9.],
       ...,
       [1., 9., 9., ..., 9., 1., 9.],
       [1., 9., 9., ..., 9., 1., 9.],
       [1., 9., 9., ..., 9., 1., 9.]], dtype=float16)

casting anything. So in your int example, things only works out because casting uint to int to uint happens to retain the values.

I think everything should pan out fine if the second array was flagged for read-write (instead of readonly), although I am not sure if you can do that right now. If it was flagged readwrite, the question is whether the ufunc machinery does/should already refuse to do the operation though, since the cast is not safe!

(It does point me to the fact that maybe ArrayMethod should use the flags instead of nin and nout like ufuncs do usually/currently.)

@mattip

OK, should have guessed that this was the change, I missed the point that this was limited to gufuncs a bit, I guess:

https://github.com/numpy/numpy/pull/15162/files#diff-7b4b6a4b3c102c8fb50985a8e6de5f6eec9c53007b57088c75216ec7c01a91a1R2778

Seems like I probably had the wrong impression that this READWRITE served the only purpose of fixing broadcasting. There is a comment alluring that READWRITE was set for more reasons, which I must have missed. But I also don't understand why that comment exists. I would think that the ufunc can indicate READWRITE if actually needed, rather than assuming that it is and doing completely unnecessary casts?

@mattip do you remember the reason for this comment (that mismatches currently with the code):

numpy/numpy/core/src/umath/ufunc_object.c

Lines 2394 to 2405 in a115369

    
               /* 
        
                * We don't write to all elements, and the iterator may make 
        
                * UPDATEIFCOPY temporary copies. The output arrays (unless they are 
        
                * allocated by the iterator itself) must be considered READWRITE by the 
        
                * iterator, so that the elements we don't write to are copied to the 
        
                * possible temporary array. 
        
                */ 
        
               _ufunc_setup_flags(ufunc, NPY_ITER_COPY | NPY_UFUNC_DEFAULT_INPUT_FLAGS, 
        
                                  NPY_ITER_UPDATEIFCOPY | 
        
                                  NPY_ITER_WRITEONLY | 
        
                                  NPY_UFUNC_DEFAULT_OUTPUT_FLAGS, 
        
                                  op_flags);

Maybe there are cases where the array iterator makes temporary copies for intermediate values when the output and input overlap. I think the comment is hinting that in order for that to work properly the iter must be marked READWRITE, but it was changed here in #15162. If that conjecture is true, I wonder why there are no overlapped memory test failures after #15162.

Hmmmgrrrr. I think the comment probably originates from where handling (although that bleeds a lot of nditer implmenetation details IMO). Now gufuncs don't actually support where handling, so that would explain why no tests can fail.

That doesn't really help me figure out if the change is right here or not (it seems right to me), or whether we want it even if it is correct. Also, it begs to question whether the where implementation that probably relies on this for some code paths isn't just broken to begin with.

PR #11580 that added the comment does mention where, but I think the point was issue #11416 about overlapping memory (we need to read old values and write new values in matmul).

OK, I am pretty sure you are right your note on overlapping input and output. When that happens and where is used, we will definitely run into this. That type of behaviour seems to be the origin of the comment (although it is older than overlap detection, there are definitely other ways to trigger the (updateif)copy. I am just not sure if the ufunc machinery will ever do it in its current state).

This is ridiculous, but shows the loss of precision in the masked values:

In [21]: out = np.full(8, 1.25, dtype=np.float64)

In [22]: out
Out[22]: array([1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25])

In [23]: arr = out.view(np.int64)  # overlaps

In [24]: np.add(arr, arr)
Out[24]: 
array([9216616637413720064, 9216616637413720064, 9216616637413720064,
       9216616637413720064, 9216616637413720064, 9216616637413720064,
       9216616637413720064, 9216616637413720064])

In [25]: np.add(arr, arr, out=out, where=[True, False]*4)
Out[25]: 
array([9.21661664e+18, 1.00000000e+00, 9.21661664e+18, 1.00000000e+00,
       9.21661664e+18, 1.00000000e+00, 9.21661664e+18, 1.00000000e+00])

You should be able to create slightly less ridiculous things if you add a float64 array in-place, but force dtype=float32.

It still seems to me that this change is strictly speaking correct (happy to be convinced to undo it though) and additionally that checking the WRITEMASKED flag inside nditer would be a better solution. But, that example with the loss of precision is still confusing me.
Of course it is almost impossible to trigger it...

PR #11580 that added the comment does mention where

Yeah, in particular in the original position of that comment, NPY_ITER_READWRITE was only set if where was used. But going before that diff, the generalized ufuncs still always use READWRITE on outputs, even though the comment (i.e. where) doesn't apply. I still think that was probably to make broadcasting work.

In theory, I think the correct thing would be:

Here, everything is fine. The gufunc should be declared as taking two inputs, but manually flag the second input to READWRITE. (considering that the output is garbage without out being passed in, it is probably OK to require it to be passed)
where= seems slightly broken, although that is very hard to trigger. I wonder if the correct fix would actually be to have a PyArray_ResolveWritebackIfCopyMasked, so that the iterator can copy back selectively (and doesn't have to copy the original data). (In principle, that dance is only necessary if an actual cast is necessary so loss of data can occur! In most cases no cast is necessary for the out argument – or buffering is used – so that copying everything "back and forth" is a valid option; The difference basically replacing a full copy back and forth, with a selective copy of only the result values.)

sounds like a plan

Navigation Menu

Repository navigation

BUG: `where=` in gufuncs is slightly broken with `out=` and casts #18700

Reproducing code example:

Error message:

NumPy/Python version information:

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

	/*
	* We don't write to all elements, and the iterator may make
	* UPDATEIFCOPY temporary copies. The output arrays (unless they are
	* allocated by the iterator itself) must be considered READWRITE by the
	* iterator, so that the elements we don't write to are copied to the
	* possible temporary array.
	*/
	_ufunc_setup_flags(ufunc, NPY_ITER_COPY \| NPY_UFUNC_DEFAULT_INPUT_FLAGS,
	NPY_ITER_UPDATEIFCOPY \|
	NPY_ITER_WRITEONLY \|
	NPY_UFUNC_DEFAULT_OUTPUT_FLAGS,
	op_flags);

BUG: where= in gufuncs is slightly broken with out= and casts #18700

Description

Reproducing code example:

Error message:

NumPy/Python version information:

Activity

seberg commented on Mar 30, 2021

seberg commented on Mar 30, 2021

mattip commented on Mar 30, 2021

seberg commented on Mar 30, 2021

mattip commented on Mar 30, 2021

seberg commented on Mar 30, 2021

seberg commented on Mar 30, 2021

mattip commented on Mar 31, 2021

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions

BUG: `where=` in gufuncs is slightly broken with `out=` and casts #18700