Skip to content

BUG,MAINT: Fix issues with non-reduce broadcasting axes #15162

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mattip merged 11 commits into from
Jun 4, 2020
Merged

BUG,MAINT: Fix issues with non-reduce broadcasting axes #15162

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
d3bfaf8
MAINT: Simplify undoing of iter-axis-perm in nditer
seberg May 11, 2020
407f626
MAINT: (nditer) move duplicated reduce axis check into helper func
seberg May 11, 2020
f31aeef
MAINT: Simplify signature of npyiter_replace_axisdata
seberg May 11, 2020
07ea525
MAINT: Simplify strides and axis check in npyiter_new_temp_array
seberg May 11, 2020
7f8dc86
ENH: Allow specifying reduction axis in nditer/NpyIter
seberg May 11, 2020
87cb35f
MAINT: Use reduction axis marker for einsum
seberg May 11, 2020
b5c2f94
BUG: Avoid incorrect broadcasts on non-core outputs in gufuncs
seberg May 11, 2020
729b9ca
BUG,MAINT: Simplify reduction result creation using reduce axis
seberg May 11, 2020
c65acda
DOC: Slightly improve error message on incorrectly shaped reductions
seberg May 11, 2020
153376f
MAINT: Address review comments related to NPY_ITER_REDUCTION_AXIS
seberg May 30, 2020
8da4fba
BUG: The reduction output must not cause the input to be broadcast
seberg May 30, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions numpy/core/src/multiarray/common.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -342,5 +342,14 @@ NPY_NO_EXPORT PyArrayObject *
new_array_for_sum(PyArrayObject *ap1, PyArrayObject *ap2, PyArrayObject* out,
int nd, npy_intp dimensions[], int typenum, PyArrayObject **result);


/*
* Used to indicate a broadcast axis, see also `npyiter_get_op_axis` in
* `nditer_constr.c`. This may be the preferred API for reduction axes
* probably. So we should consider making this public either as a macro or
* function (so that the way we flag the axis can be changed).
*/
#define NPY_ITER_REDUCTION_AXIS(axis) (axis + (1 << (NPY_BITSOF_INT - 2)))

#endif

6 changes: 2 additions & 4 deletions numpy/core/src/multiarray/einsum.c.src
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2740,7 +2740,7 @@ PyArray_EinsteinSum(char *subscripts, npy_intp nop,
op_axes[nop][idim] = idim;
}
for (idim = ndim_output; idim < ndim_iter; ++idim) {
op_axes[nop][idim] = -1;
op_axes[nop][idim] = NPY_ITER_REDUCTION_AXIS(-1);

@eric-wieser eric-wieser May 30, 2020

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will this change make np.einsum("...i,...i->...i", np.arange(3), np.array([1])) an error? That is, should named dimensions in einsum be viewed as non-broadcastable?

@seberg seberg May 30, 2020

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, it changes that. I can undo. But einsum does not really have keepdims semantics nor is that notation einsum notation AFAIK.
Completely changing the resul behaviour by providing (a potentially wrong) output array seems undesirable to me?

@eric-wieser eric-wieser May 30, 2020 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What's the resulting error look like with this branch?

Not sure what you mean by that not being einsum notation - that notation works for me

@seberg seberg May 30, 2020

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, this is unchanged, I thought you had np.einsum("...i->...i", np.arange(3), out=np.array([1]))

If you do not provide an out array, no semantics are ever changed here.

@seberg seberg May 30, 2020

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, sorry... I do not think anything changed at all. Since opposite broadcasting should be the only thing that changed (since only np.newaxis is marked here). But einsum specifically checks for out having too many dimensions.

@mhvk mhvk May 30, 2020

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That sounds good. But as a check, would one or the other solve the no-input-but-have-output case discussed in #15140?

@mhvk mhvk May 30, 2020

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, I had a quick look and this PR does not address #15140, so no worries about it!

@seberg seberg May 30, 2020

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmmmm. I am wondering if we could mark an operand as output to indicate that it should not affect shape (or maybe using the WRITABLE flag is OK, but I fear that is a bit of a potential regression).

I suppose when there is no shape, the shape is just should be () (and if there are reduce, but keepdims, axis, they are always inserted with size 1 right now).

With your gufunc the question is whether you actually want the output to define the size in that case. Technically, that could be nice in some cases, unless you argue that the user should somehow provide shape=... in that case.
Although, for the iterator maybe that would be the same, and the gufunc machinery would just not pass the flag if there are no inputs?

@seberg seberg May 30, 2020 •
edited
Loading

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, I am conflating two issues there. With no arguments we could already simply default to a 0-D iterator without much hassle probably...
How/if the output argument can affect the iterator shape is another issue.

(edit, but for another PR.)

@mhvk mhvk May 30, 2020

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed, for my strange case I do not think the output shape has to affect the iterator; it was about automatically allocating gufunc dimensions if needed, even if there are no inputs.

}

/* Set the iterator per-op flags */
Expand All @@ -2753,13 +2753,11 @@ PyArray_EinsteinSum(char *subscripts, npy_intp nop,
op_flags[nop] = NPY_ITER_READWRITE|
NPY_ITER_NBO|
NPY_ITER_ALIGNED|
NPY_ITER_ALLOCATE|
NPY_ITER_NO_BROADCAST;
NPY_ITER_ALLOCATE;
iter_flags = NPY_ITER_EXTERNAL_LOOP|
NPY_ITER_BUFFERED|
NPY_ITER_DELAY_BUFALLOC|
NPY_ITER_GROWINNER|
NPY_ITER_REDUCE_OK|
NPY_ITER_REFS_OK|
NPY_ITER_ZEROSIZE_OK;
if (out != NULL) {
Expand Down
16 changes: 6 additions & 10 deletions numpy/core/src/multiarray/nditer_api.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -935,13 +935,8 @@ NpyIter_GetShape(NpyIter *iter, npy_intp *outshape)
if (itflags&NPY_ITFLAG_HASMULTIINDEX) {
perm = NIT_PERM(iter);
for(idim = 0; idim < ndim; ++idim) {
npy_int8 p = perm[idim];
if (p < 0) {
outshape[ndim+p] = NAD_SHAPE(axisdata);
}
else {
outshape[ndim-p-1] = NAD_SHAPE(axisdata);
}
int axis = npyiter_undo_iter_axis_perm(idim, ndim, perm, NULL);
outshape[axis] = NAD_SHAPE(axisdata);

NIT_ADVANCE_AXISDATA(axisdata, 1);
}
Expand Down Expand Up @@ -1005,16 +1000,17 @@ NpyIter_CreateCompatibleStrides(NpyIter *iter,

perm = NIT_PERM(iter);
for(idim = 0; idim < ndim; ++idim) {
npy_int8 p = perm[idim];
if (p < 0) {
npy_bool flipped;
npy_int8 axis = npyiter_undo_iter_axis_perm(idim, ndim, perm, &flipped);
if (flipped) {
PyErr_SetString(PyExc_RuntimeError,
"Iterator CreateCompatibleStrides may only be called "
"if DONT_NEGATE_STRIDES was used to prevent reverse "
"iteration of an axis");
return NPY_FAIL;
}
else {
outstrides[ndim-p-1] = itemsize;
outstrides[axis] = itemsize;
}

itemsize *= NAD_SHAPE(axisdata);
Expand Down
Loading