ENH: implement np.unwrap as a generalized ufunc - #31848

FYI this is what one gets with the current unwrap testing on main. Some error messages are different here versus main and naturally on main subclasses are not returned. But I'd argue for the cases that unwrap shouldn't work anyways, the error messages are a bit better now instead of overflow errors and TypeError: ufunc 'remainder' not supported

FAILED numpy/lib/tests/test_function_base.py::TestUnwrap::test_subclass_preserved[float32] - AssertionError: assert False
FAILED numpy/lib/tests/test_function_base.py::TestUnwrap::test_subclass_preserved[float64] - AssertionError: assert False
FAILED numpy/lib/tests/test_function_base.py::TestUnwrap::test_subclass_preserved[int8] - AssertionError: assert False
FAILED numpy/lib/tests/test_function_base.py::TestUnwrap::test_subclass_preserved[int32] - AssertionError: assert False
FAILED numpy/lib/tests/test_function_base.py::TestUnwrap::test_subclass_preserved[int64] - AssertionError: assert False
FAILED numpy/lib/tests/test_function_base.py::TestUnwrap::test_subclass_preserved[bool] - AssertionError: assert False
FAILED numpy/lib/tests/test_function_base.py::TestUnwrap::test_subclass_preserved[object] - AssertionError: assert False
FAILED numpy/lib/tests/test_function_base.py::TestUnwrap::test_complex_raises - TypeError: ufunc 'remainder' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting r...
FAILED numpy/lib/tests/test_function_base.py::TestUnwrap::test_unsigned_integer_period_raises[uint8] - OverflowError: Python integer -4 out of bounds for uint8
FAILED numpy/lib/tests/test_function_base.py::TestUnwrap::test_unsigned_integer_period_raises[uint32] - OverflowError: Python integer -4 out of bounds for uint32
FAILED numpy/lib/tests/test_function_base.py::TestUnwrap::test_unsigned_integer_period_raises[uint64] - OverflowError: Python integer -4 out of bounds for uint64

IMO this should be ready for a review pass

Thanks for the review! @ngoldbaum I'll address everything in the coming week as best as I can.
In general I was mostly trying to preserve behavior with main as much as I could but couldn't fully preserve the exceptions indeed. I guess because unwrap didn't really care too much about errors from how I read the imlementation and was like "yeah whatever the other functions raise". You also pointed out some real bugs in the ufunc code which I will fix obviously.

@mhvk what are your thoughts on 9a45e86 regarding the discont parameter? I asked claude to make it NEP-50 compartible and gave it the NEP.

Or I have one more idea. We get rid of the whole python promotion and passing signature in the ufunc and do the promotion with a promoter in c++. @ngoldbaum and @mhvk what do you think?

I agree with Marten it doesn't seem to make much sense to keep e.g. discont double if the input is float32; I don't think we're gaining any real precision as this is an elementwise comparison, and there is no error accumulation / propagation. I would simply cast it to the input array dtype (you could make an exception for integers to preserve boundary precision, though there could be better ways to handle the boundary). NEP 50 applies to promoting Python scalars especially when the dtype is "real" (ends up in a result for example) but this seems to be on the other end of a specific internal mechanic. In any case, I'm pretty sure this new loop registration is complicated and not right just to keep discont's type.

Because I know very little about implementing promoters, I had claude generate me one to see the idea and pushed it in another branch. This is the commit: ikrommyd@c32bc51. Do you think we should go into that direction?

I still don't think we should care about specializing precision for discont, and that branch just moves the same logic into C. I think that should probably be resolved first. To be clear, I think a C promoter is worth it if it is more simply correct (for both legacy behavior and user dtypes), not to solve that problem.

Main's output dtype is result_type(p, period) only, never discont, so folding discont into that promotion would change the return dtype based on a threshold value, which is worse than what we're fixing. But discont's own comparison precision on main already follows NEP 50 (weak scalars downcast to p's dtype, explicit ones promote), so a gufunc that always compares discont at a fixed type (double, or always T) diverges from main in one direction or the other. For example if p=[0,3] with float16 dtype, and discont=3.0001, main gives [0,-1], while fixed-double gives [0,3]. That's why I think discont needs its own resolution, separate from T. Doing it in python or c++ does not matter much.

I agree casting to T is technically a behavior change, just not sure if we shouldn't make it. I also don't think it's going to affect weak scalars (as tested), just explicit ones.

But even right now, an explicit longdouble discont would still lose precision. Then we would need at least 3 more loops to bring this to its logical completion!

I'd be fine with casting to T for floats but we can't do it for integer Ts

Is it a thing to do multiple release notes like an improvement + another one for behavior changes or should everything go into the same ‎doc/release/upcoming_changes/9959.improvement.rst?

Two different release notes are a thing, you can add a second one for the changes.

@ngoldbaum @MaanasArora I believe I have addressed the comments and have added the converged dtype logic in 86a2566
The riscv CI failure is unrelated infrastructure problem.
Regarding Nathan's comments about UB with signed integers, I think these should probably all be fixed centrally in numpy at some point and not do anything here as no numpy additions/subtractions that I saw care about UB so far.

Thank y'all for the reviews! Cheers!