ENH: Make it possible to NpyIter_RemoveAxis an empty dimension #3861
Conversation
| NIT_ADVANCE_AXISDATA(axisdata, 1); | ||
| } | ||
| NIT_ITERSIZE(iter) = size; | ||
| } |
There was a problem hiding this comment.
Can we just use the second code path unconditionally? (With a comment noting that recomputing from scratch is safest because it's possible the removed axis was zero-sized.) No point in having multiple code paths where one will do.
|
Aside from that comment LGTM. Simplifies the code and makes it more general, nice. |
|
ERROR: test_0_size (test_linalg.TestInv) Traceback (most recent call last): File "/home/travis/virtualenv/python2.6/lib/python2.6/site-packages/numpy/linalg/tests/test_linalg.py", line 293, in test_0_size res = linalg.inv(a) File "/home/travis/virtualenv/python2.6/lib/python2.6/site-packages/numpy/linalg/linalg.py", line 503, in inv ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj) ValueError: On entry to ZGESV parameter number 4 had an illegal value |
|
Ping Travis. |
|
Yeah, the test failure is real. I removed the special case, and whatever is linked locally can handle the empty operation in BLAS, but on travis it can't. So I have to plug some kind of special case(s) back in there. Not sure yet on which level. |
|
@seberg Where is this PR after all the other changes in nditer? |
|
It has some code overlap (gets simpler then it is currently), but feature wise it has nothing to do with the rest. Still would have to see which gufuncs support empty arrays and which don't though... |
c899f31 to
21af68e
Compare
|
This should be fine, but since lapack does not support empty arrays (at least mine apparently), all those functions would need to be looked at... |
|
@seberg I've forgotten, if I ever knew, what this is about. Seems like a pretty small patch and if it fixes other problems I'd like it to go in. |
|
Its about gufuncs not being able to handle empty arrays. unfortunatly, lapack doesn't handle, so that module needs to be carefully checked ebfore putting this in. |
|
OK, I have updated this, so some of the things now accept 0-d such as determinant. Cholesky I am not sure, etc. should have a bit careful look over the test case and everything that was changed has a test case etc. (i.e. I am a bit worried about the linalg changes, less so the nditer ones). |
|
☔ The latest upstream changes (presumably #6115) made this pull request unmergeable. Please resolve the merge conflicts. |
|
@seberg Could you rebase this? I'm going to branch 1.10 tomorrow, and if you think this is ready I'll put it in. |
|
Probably is pretty much, but I am on vacation and do not have my compiter with me. But maybe it is better to wait for bew master anyway, in case I missed a function. |
3e48ef9 to
8a093ab
Compare
|
Rebased this, but if someone has the patience to double check the linalg stuff, maybe would not hurt. |
|
☔ The latest upstream changes (presumably #6905) made this pull request unmergeable. Please resolve the merge conflicts. |
|
☔ The latest upstream changes (presumably #7073) made this pull request unmergeable. Please resolve the merge conflicts. |
d46c121 to
3e62d56
Compare
| @@ -1010,8 +991,6 @@ def eigvalsh(a, UPLO='L'): | |||
| _assertRankAtLeast2(a) | |||
| _assertNdSquareness(a) | |||
| t, result_t = _commonType(a) | |||
| if _isEmpty2d(a): | |||
There was a problem hiding this comment.
There was a bug here for complex value, since it missed the cast to real (I think in one other function as well).
There was a problem hiding this comment.
@eric-wieser can you have a look if you have time. Since I don't know all your newer changes, and I guess it may be that the tests are not quite how they would be nice (dunno how well the rebase worked). Also this moves some of your fixups into the C-lapack wrappers.
There was a problem hiding this comment.
Looks generally good, but a couple of comments:
- A lot of lapack functions do support
n == 0, but have an argument that we currently pass asnthat we should be passing asmax(n, 1)according to the docs - The test-suite already tests some empty arrays through the
domethod, so maybe the test_size_0 could be merged, if we really care (not sure we do). At any rate, I think that adding some dtype-changing tests to the main set of tests would be sensible. Perhaps for another PR though
| gufunc = _umath_linalg.solve1 | ||
| else: | ||
| if b.size == 0: | ||
| if (a.shape[-1] == 0 and b.shape[-2] == 0) or b.shape[-1] == 0: | ||
| a = a[:,:1].reshape(a.shape[:-1] + (1,)) |
There was a problem hiding this comment.
This looks to me like it was wrong in the first place, and should have been a[...,:1]
numpy/linalg/umath_linalg.c.src
Outdated
| LAPACK(@cblas_type@getrf)(&m, &m, (void *)src, &m, pivots, &info); | ||
| if (m != 0) { | ||
| /* Lapack (at least some versions) do not handle n == 0 */ | ||
| LAPACK(@cblas_type@getrf)(&m, &m, (void *)src, &m, pivots, &info); |
There was a problem hiding this comment.
Which versions are we supporting? Even the oldest 3.0.0 version in lapack_lite we're using supports this, but the third m should be max(1,M) according to the documentation
There was a problem hiding this comment.
Dunno, my guess was I had problems on my system, or, I just did not know if it was necessary or not. The max(1, m) thing may just be the reason too I suppose.
There was a problem hiding this comment.
@eric-wieser is this the case for all of these functions? Or which one would need this kind of max handling? I will add a commit undoing it and see what happens for now, and we can still go back and forth if problems show up with some lapack or so.
There was a problem hiding this comment.
OK, it worked here, but for none of the other cases, but I guess you expected that?
There was a problem hiding this comment.
Though, I did not try to replace parameters with 1 instead of 0 there. Is it worth checking that out?
There was a problem hiding this comment.
You'd need to look at the documentation for each lapacj function in turn, which can be found within f2c_[czds]_lapack.c
There was a problem hiding this comment.
Hmmm, ok, should do some time I guess, but don't have the nerve for it right now.
| op += os_p; | ||
| } | ||
| op -= ob_p; | ||
| op += os_m; |
| } | ||
| END_OUTER_LOOP | ||
| return; | ||
| } |
There was a problem hiding this comment.
Wouldn't a better fix be to just reassign the pointers, and run the below code?
Something like
if (dn == 0) {
@typ@ ip_val = 0;
@typ@ ip2_val = 0;
ip1 = ip2 = &ip_val;
// set strides to 0
dn = 1
}There was a problem hiding this comment.
Nice hack, didn't think about it much.
There was a problem hiding this comment.
Hmm, no, the char *ip1=args[0], *ip2=args[1], *op=args[2]; line is important where it is, since the macro will manipulate args
Yeah, I would have thought you want the loop over n innermost, and then your problems would go away
Never mind then. |
|
Yeah, with the innermost loop over n, this just falls out: BEGIN_OUTER_LOOP_3
char *ip1=args[0], *ip2=args[1], *op=args[2];
for (m = 0; m < dm; m++) {
for (p = 0; p < dp; p++) {
@typ@ op_val = 0;
for (n = 0; n < dn; n++) {
@typ@ val1 = (*(@typ@ *)ip1);
B37A
@typ@ val2 = (*(@typ@ *)ip2);
op_val += val1 * val2;
ip1 += is1_n;
ip2 += is2_n;
}
*(@typ@ *)op = op_val;
ip1 -= ib1_n;
ip2 -= ib2_n;
ip2 += is2_p;
op += os_p;
}
ip2 -= ib2_p;
op -= ob_p;
ip1 += is1_m;
op += os_m;
}
END_OUTER_LOOPPerhaps this has performance/caching implications - I guess the old code iterated over I guess an optimal implementation would look at the strides and do a different loop in each case. For instance, my implementation would be more efficient at computing |
443ba37 to
6f75eb3
Compare
|
☔ The latest upstream changes (presumably #8847) made this pull request unmergeable. Please resolve the merge conflicts. |
6f75eb3 to
fe3e16c
Compare
|
This didn't previously segfault, did it? Bad rebase? |
|
One failure is in this test |
doc/release/1.13.0-notes.rst
Outdated
| @@ -120,6 +140,8 @@ behaviour is recovered if ``axis=None`` (default). | |||
| ``np.gradient`` now supports unevenly spaced data | |||
| ------------------------------------------------ | |||
| Users can now specify a not-constant spacing for data. | |||
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | |||
| Users can now specify a not-constant spacing for data. | |||
doc/release/1.13.0-notes.rst
Outdated
| @@ -131,11 +153,13 @@ In particular ``np.gradient`` can now take: | |||
|
|
|||
| This means that, e.g., it is now possible to do the following:: | |||
|
|
|||
| This means that, e.g., it is now possible to do the following: | |||
doc/release/1.13.0-notes.rst
Outdated
| @@ -227,6 +251,8 @@ More ``linalg`` operations now accept empty vectors and matrices | |||
| All of the following functions in ``np.linalg`` now work when given input | |||
| arrays with a 0 in the last two dimensions: ``det``, ``slogdet``, ``pinv``, | |||
| ``eigvals``, ``eigvalsh``, ``eig``, ``eigh``. | |||
| arrays with a 0 in the last two dimensions: `det``, ``slogdet``, ``pinv``, | |||
| ``eigvals``, ``eigvalsh``, ``eig``, ``eigh``, ``cholesky``. | |||
There was a problem hiding this comment.
Ooops, I am sorry, I am a bit overwhelmed with other things right now, not sure I will be checking on things regularly for a while :/.
There was a problem hiding this comment.
OK, that should have fixed it up, when going back to the previous version of setting to zero, I got the pointer arithmetic wrong apparently.
|
☔ The latest upstream changes (presumably #8999) made this pull request unmergeable. Please resolve the merge conflicts. |
The caller should take of not using pointers given in such a case
99ce747 to
ca8b952
Compare
The necessary fixup on the C-side of linalg has been done already (i.e. the gufuncs correctly work for these empty arrays). This also enables cholesky decomposition and fixes a small bug in pinv handling. Co-authored-by: Eric Wieser <wieser.eric@gmail.com>
The case where the result is not empty was missing the zeroing out now that it is legal to use that input.
ca8b952 to
1c4e0d4
Compare
|
Nice job! Based on the silence on the waiting list, this looks pretty safe to merge? |
|
I don't have a problem with this getting merged. Support for empty arrays is not complete and this is another step on the way, |
|
OK; lets give it a shot. |
|
On the subject of empty array kludges - what is going on here? I can't see how the iteration can be empty, yet the results not be |
|
Apparently not a single one of our tests hits that code |
|
Sounds like it was planned to catch those empty array matrix multiplies for us, but the if is incorrect. Since this is not really a reduction.... I would tend to leave it as is after this now (i.e. let the inner loop worry about it) and remove the code block (assuming, that we are not missing something really silly). |
|
Doesn't look like it's checking for the right thing - that's looking at the |
|
Yeah. |
This is the followup to that "kludge" discussion. I did not really check much if there might be something going awefully wrong, though on the other hand it is hard to imagine something too bad.
There are certainly more kludges to remove here assuming that it works. So should check for those first.