Skip to content

NEP: Add NEP 50 draft about fixing scalar promotion rules #21103

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.

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
merged 20 commits into from
May 25, 2022
Merged

Conversation

seberg
Copy link
Member

@seberg seberg commented Feb 21, 2022

Note, this really is a very early draft. For myself, I have not
even 100% made up my mind that this is actually the best approach.

However, to solicitate feedback and input, I would like to put it
out there (and have a nice rendered version).

For those mainly interested in the possible different approaches
(i.e. making decisions), the "Alternatives" section is probably most
interesting.

Please make a note of blatend errors or unnecessary repitition, but
don't be too surprised if you see them, I will not be surprised if
this needs to be redone from scratch eventually.


Artifact link: https://25381-908607-gh.circle-artifacts.com/0/doc/neps/_build/html/nep-0050-scalar-promotion.html (Updated 2022-03-16)

(early) branch for testing: https://github.com/seberg/numpy/tree/weak-scalars (Note that this branch currently does not attempt to give warnings, or implement the uint8_arr + 1000 -> error behavior, which is more complicated. It should implement the general "weak" logic though.

@jni
Copy link
Contributor

jni commented Feb 22, 2022

For what it's worth, your summary from the mailing list post:

Roughly speaking, I have made a proposal under the 3 points:

  • NumPy scalars and NumPy arrays always behave the same.
  • A NumPy array always respects the dtype
  • A Python scalar is "weak" so that uint8_arr + 3 returns a uint8_arr

reflects exactly how I voted, and I think is the best approach, specifically because (as far as I can tell) Python scalars are indistinguishable from literals and I would expect uint8_array + 3 to remain uint8. So a big 👍 from me! Thanks @seberg for tackling the Hard Problems! 😊

@seberg
Copy link
Member Author

seberg commented Feb 22, 2022

Thanks, I do encourage you to have a brief look at the alternatives (if you have not done before), because some stuff may also just be a way to try for fewer backcompat breaks, even if we don't think it is the purest approach...

Note that the examples apply just as well for operations like multiplication,
addition, or comparisons and the corresponding functions like `np.multiply`.

This NEPs proposes to refactor the behaviour around two guiding principles:
Copy link
Contributor

Choose a reason for hiding this comment

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

NEPs -> NEP I think

2. NumPy scalars or 0-D arrays must always lead to the same behaviour as
their N-D counterparts.

We propose to removes all value-based logic and add special handling for
Copy link
Contributor

Choose a reason for hiding this comment

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

removes -> remove

We propose to removes all value-based logic and add special handling for
Python scalars to preserve some of the convenience that it provided.

This changes also apply to ``np.can_cast(100, np.int8)``, however, we expect
Copy link
Contributor

Choose a reason for hiding this comment

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

This -> These maybe

two examples. For the second one, the results will be the following::

np.arange(10, dtype=np.int8) + 256 # raises a TypeError
nparray([1., 2.], dtype=np.float32) * 1e200 # warning and returns infinity
Copy link
Contributor

Choose a reason for hiding this comment

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

np.array

Copy link
Contributor

Choose a reason for hiding this comment

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

Not sure how much value there is in it, but for the second changed case that returns infinity, are there examples from other frameworks/programming languages that do something similar in a matching scenario?

I don't know that we'd necessarily want to mimic what other major languages/projects do, though citing them as examples may help make a stronger case if needed. Perhaps in this case one could just quickly verify in C or C++ or whatever.

Copy link
Member Author

Choose a reason for hiding this comment

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

I added an admonition, C normally does this, but of course you can set the evironment to raise an error (just like you can in NumPy, although casts currently don't check that normally!).

I decided to promise to fix that it gives a warning right now. pytorch seems to just overflow silently, I guess they probably don't give warnings reliably, because some hardware they care about doesn't really support it well.


* `JAX promotion`_ also uses the weak-scalar concept. However, it makes use
of it also for most functions. JAX further stores the "weak-type" information
on the array: ``jnp.array(1)`` remains weakly typed.
Copy link
Contributor

Choose a reason for hiding this comment

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

Heh, the JAX promotion docs are pretty detailed. I was thinking we might want to briefly mention if our new/proposed direction was consistent with some other major languages/frameworks, though not sure if that changes much.

Copy link
Member Author

@seberg seberg Feb 22, 2022

Choose a reason for hiding this comment

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

Yeah, JAX uses the weak logic, I think I should be reporting some bugs to them, but I don't remember (I believe if you swapped the inputs on some of their example, things got wrong, and JAX tags the "was scalar" to the 0-D array, so it can be passed around, which I am not sure always turns out right, but is easier for them because their arrays are immutable.

I believe most use "weak" logic (torch definitely does), although all of these skip scalars entirely probably. EDIT: Sorry, to be clear, they just don't have "scalars" and at least torch doesn't care to warn users about bad conversions.

@jakevdp
Copy link
Contributor

jakevdp commented Feb 24, 2022

Hi all - @shoyer pointed me to this PR. I wanted to mention that I recently published a detailed design doc outlining the decisions behind JAX's type promotion system, which you can read here: https://jax.readthedocs.io/en/latest/design_notes/type_promotion.html

```
* Certain operations are expected to start failing:
```python3
np.array([1], np.uint8) * 1000
Copy link
Contributor

@jakevdp jakevdp Feb 24, 2022

Choose a reason for hiding this comment

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

Can you clarify, this fails because 1000 is not representable in uint8, correct? But am I understanding correctly that np.array([1], np.uint8) * 100 would silently overflow, because 100 is representable in uint8?

My main area of confusion from the set of examples here is how to predict which operations will silently overflow and which will error.

Copy link
Member Author

Choose a reason for hiding this comment

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

Yes, that is what I have in mind currently. I am not at all set on it though – it would actually be much simpler to not do that!
Basically, my rough plan would be that the above would fail, but np.uint8(1000) would not for now, so that you could still write np.uint8(-1) if you must for some reason – because... backward compatibility mostly.

That is painful, so a simpler middle ground could also be that both np.uint8(1000) and np.uint8(1) + 1000 would give the same warning.

In general, NumPy does not give reliable warnings for integer overflows. We only give warnings for scalar python operations (and even then randomly, because scalar operator calls are poorly organized). So np.uint8(10) + 250 will give a warning, but that example would not (it works with an array).

That is one argument for why things could be a smoother transition if we limit to Python operators: For those (and scalars which are currently the place that upcasts a lot), we could pull off to ensure warnings are given. (Super painful because it requires cleaning up the scalar mess, but a clear path.)

One thing that I am unsure about what JAX does, is that I don't really think we should tag on the "weak" property to the array as JAX does. One, because I am worried it may be brittle/surprising to carry it around for long and second because unlike JAX our arrays are mutable, so there is an additional complexity. (We can still tag on something as an implementation detail.)

Copy link
Contributor

@jakevdp jakevdp Feb 24, 2022

Choose a reason for hiding this comment

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

Thanks - regarding JAX & the weak type attribute: the reason we used that approach is that in JAX's JIT-compiled code, everything is represented by a JAX array (scalar inputs as well). So to maintain the same behavior between JIT and non-JIT code paths, we need a way to mark a JAX array as behaving similarly to a Python scalar: the mechanism we chose is a weak_type attribute on the array's abstract value. Since NumPy doesn't have the same requirement, it probably makes sense not to allow for weakly-typed numpy arrays.

That said, one benefit of JAX's approach is that calling jnp.asarray on a function's inputs does not change the results even if the input is a scalar, because asarray maintains the weakness of the input scalar. In numpy's case np.asarray(1.0) would no longer have Python scalar semantics.

On the other hand, allowing zero-dimensional arrays to be weak readily leads to cases where larger arrays can be weak as well (e.g. jnp.broadcast_to(jnp.array(1.0), (100, 100)) returns a weakly-typed matrix as of the current JAX release) so weak types can no longer be thought of as strictly synonymous with scalars. That's a bit surprising, and is one reason we're currently actively discussing design considerations around the semantics of weakly-typed values and how the weak_type attribute propagates.

Copy link
Member Author

Choose a reason for hiding this comment

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

Ah, that is interesting (and I should mention the weak attribute in alternatives, even if I currently would need some convincing to consider it more seriously – and add 1-2 more examples here for clarity of the initial point, of course).

The broadcast example is nicely clear about how things may get confusing when dragging on the flag. Long ago, I had early plans to support, in theory, even [1, 2, 3] to be considered "weak" – e.g. for ufuncs – but I had abandoned them: it seemed just to much and also way too many moving parts at the time.

In the end, it would be nice to have a clear concept of how far the "weakness" reaches. I could half imagine an np.asarray(x, allow_weak=True) for library use, but I somewhat imagine libraries are better of handling it up-front with np.result_type (or by not calling asarray at all) – possibly asarray_unless_literal(x) if need be.

That's a bit surprising, and is one reason we're currently actively discussing design considerations around the semantics of weakly-typed values.

It will be interesting to hear about that! NumPy might end up being a bit more limited in options due to backcompat, but it would be nice if we can manage to align some concepts, i.e. at which step "weakness" is forgotten.

Copy link
Contributor

@mhvk mhvk left a comment

Choose a reason for hiding this comment

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

This is a great start! I very much like removing distinctions between scalars and 0 or any D arrays.

I also agree that there is least surprised by interpreting python values as literals for which, effectively, one does something like np.array(python_literal) and, if possible, does .astype(self.dtype, casting='same_kind') (i.e., allow cast from int to float and float to complex, but not the other way around).

It may be good to be explicit about what will happen for:

np.arange(10, dtype='i2') * 1.
np.arange(10, dtype='f4') * 1j

Surely, the answer should be float and complex, but should there be an attempt to stay as close as possible to the dtype? I.e., float32 and complex64? Especially the latter would seem nice (though my use case would be better served by having an expj ufunc!).


This would be the simplest solution, however, it would lead to many upcasts when
working with arrays of ``float32`` or ``int16``, etc. The solution for these cases
would be to rely on in-place operations.
Copy link
Contributor

Choose a reason for hiding this comment

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

It may be good to write it as an explicit advantage of the proposed scheme that operations with a python scalar give the same result whether or not it is done in-place.

@mhvk
Copy link
Contributor

mhvk commented Mar 1, 2022

An alternative perhaps worth discussing (in the considered but rejected section) is that for python scalars values do matter, and that they always are interpreted as the smallest dtype that can hold them. I think that can be rejected given the immediate problems with numbers like 0.1 - making that a float16 would be pretty bad as np.arange(10.) + 0.1 would become equivalent to

np.arange(10.)+np.float16(0.1)
array([0.09997559, 1.09997559, 2.09997559, 3.09997559, 4.09997559,
       5.09997559, 6.09997559, 7.09997559, 8.09997559, 9.09997559])

while on the other hand one would not want to just go for np.float64 unconditionally. So, sticking to the dtype of the array would make more sense.

@mhvk
Copy link
Contributor

mhvk commented Mar 1, 2022

p.s. I realize some of my examples are well answered in @jakevdp's write-up (very nice!). I suggest adding something like the lattice there, as it really helps to visualize the promotions.

@seberg
Copy link
Member Author

seberg commented Mar 8, 2022

Btw, I have an early testing branch here: https://github.com/seberg/numpy/tree/weak-scalars

I think I covered the central parts, but it is fairly minimal and the branch does not attempt to:

  • Give an optional warning when a different result occurs (or is very likely).
  • Raise an error for np.uint8(1) + 1000 because it will use unsafe/same-kind casting on the 1000. (The most tricky part probably!)
  • There may be some corner cases missing for scalars (the scalars are a mess that need cleanup, but I can't find a broken example right now)

As expected, the scipy test failures do not seem intimidating (a huge chunk is in sparse, I expect they are not weird). But I am not surprised on that front, since I had tested it before.
In general, the impact is far more interesting for end-users than libraries I think. Aside maybe pandas because it effectively has a lot of "parallel" infrastructure.

(The vast majority of failures are in sparse, sparse needs to "realign" with NumPy)

@rgommers
Copy link
Member

As expected, the scipy test failures do not seem intimidating (a huge chunk is in sparse, I expect they are not weird).

@seberg do you expect it'll be useful to look into those failues? If a good fraction of them point at things that could be improved in SciPy or simply rewritten to avoid the failures, do you think that helps either your effort here or SciPy itself?

@seberg
Copy link
Member Author

seberg commented Mar 10, 2022

@rgommers I did not dig deep yet, but I think the vast majority of failures fall into three categories:

  • The tests simply need adjusting: they are thorough and check for result types, but result types change. I think (I am not quite sure!) that the sparse tests fall into this category (at least mainly).
  • A few floating point equality (or approximate equality) checks fail because np.float32(...) == 123.32123 casts the Python float to float32 (it is considered "weak") rather than vice-versa, and that requires adaptation of test sensitivity.
  • E.g. integrators sometimes pass on scalars, if the input is float32 the inner function would have ended up using float64 anyway (due to the upcast of scalars). If that changes, that means some results have reduced precision and tests fail.

Right now I expect most of these are tedious, but not particularly problematic to resolve for SciPy (and many other libraries). The impact seems much more on the end-user side rather than for libraries.

However, I am not sure if SciPy or other libraries will want to support Python scalars as arguments (without calling np.asarray(argument) indiscriminately). If libraries want that, this may be a very important user-story to figure out and propose a solution for in the NEP. (I don't think it is hard to find a solution, but I am not sure the exact direction to go right now.)

@seberg seberg changed the title NEP: Add early draft for fixing scalar promotion rules NEP: Add NEP 50 draft about fixing scalar promotion rules Mar 15, 2022
@rgommers
Copy link
Member

However, I am not sure if SciPy or other libraries will want to support Python scalars as arguments (without calling np.asarray(argument) indiscriminately).

As arguments where? There are certainly APIs that take a Python scalar and use it I'd think. First example I found: stats.iqr uses a scale keyword with is a Python scalar, and it's used to scale the output array with.

@seberg
Copy link
Member Author

seberg commented Mar 15, 2022

OK, I made a fairly big update to the first part (and some additional chunks to the end/back compat that were brought up here). Including adding schema and examples (because it is probably super confusing that the "integral"/"inexact" boundary matters so much, but I think that is an important boundary, because it is what we currently have.)

As arguments where? There are certainly APIs that take a Python scalar and use it I'd think.

Well, I am thinking of "ufunc-like" functions mainly, I guess. The point is that:

np.add(np.uint8(2), 1)

would now use a "weak" 1 (EDIT: returning np.uint8(3)). But if you would wrap that into a "ufunc like" function:

def my_add(a, b):
    a = np.asarray(a)
    b = np.asarray(b)
    return a + b

my_add(np.uint8(2), 1)

Would return an int64! So if such a function should "respect" the weak Python scalar it would have to remove the asarray call, or call np.result_type up-front, or something else?

@rgommers
Copy link
Member

Well, I am thinking of "ufunc-like" functions mainly, I guess.

Ah okay, then I can't think of any right now - such functions accept scalars only as array-like and call asarray on them I expect.

Copy link
Contributor

@rossbar rossbar left a comment

Choose a reason for hiding this comment

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

Thanks for putting this together @seberg . I left some inline comments/suggestions from a first readthrough. I'm going to let it marinate then reread it again!

We propose to remove all value-based logic and add special handling for
Python scalars to preserve some of the convenience that it provided.
This will mean that Python scalars are considered "weakly" typed.
A Python integer will always be converted to the NumPy dtype retaining the
Copy link
Contributor

Choose a reason for hiding this comment

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

Should this statement be generalized to all Python scalars, i.e.

Suggested change
A Python integer will always be converted to the NumPy dtype retaining the
A Python scalar will always be converted to the NumPy dtype retaining the

:ref:`NEP 42 <NEP42>`.
Currently, this leads to a dual implementation of a new and an old (value
sensitive) system. Fixing this will greatly simplify the internal logic
and make results more consistent.
Copy link
Contributor

Choose a reason for hiding this comment

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

You mention the existence of parallel systems, but it might also be worth including language (either here or in another bullet) about improving consistency between "array logic" and "scalar logic". This is alluded to below, but I think adding a bullet to be more explicit could be an improvement (if this makes sense)

Copy link
Member Author

Choose a reason for hiding this comment

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

Well, technically the current system is consistent between arrays and scalars. It is inconsistent between 0-D and not 0-D, so not sure...


np.array([100], dtype=np.uint8) + 200 == np.array([44], dtype=np.uint8)

Considering this, we believe that the proposal follows the "principle of least
Copy link
Contributor

Choose a reason for hiding this comment

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

❤️

Copy link
Contributor

Choose a reason for hiding this comment

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

I don't understand this, as the table above shows that this behaviour is maintained?

Silent integer overflow is a viciously surreptitious source of bugs, and if I read the NEP correctly, we're increasing the surface for that by defaulting to lower integer types - e.g. uint8(1) + 2: int64(3) (old) -> uint8(3) (new)


But the following will then be a surprise::

np.array([100], dtype=np.uint8) + 200 == np.array([44], dtype=np.uint8)
Copy link
Contributor

Choose a reason for hiding this comment

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

This highlights an important point that could maybe be emphasized further (potentially higher up in the document): the value that is inspected to determine the output type is not the result of the operation but the scalar input. IMO this example really drives home one of the biggest motivations for moving away from value-based behavior!

potential changes in behavior in most relevant cases.


Impact on ``can_cast``
Copy link
Contributor

Choose a reason for hiding this comment

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

It might be worth reordering the subsections in order of "most impactful", e.g. move the can_cast subsection to the end.

OTOH maybe it's better to have the simpler case(s) with the demonstrable behavior change first. Just thinking out loud...

Comment on lines +196 to +307
This removes currently surprising cases. For example::

np.arange(10, dtype=np.uint8) + np.int64(1)
# and:
np.add(np.arange(10, dtype=np.uint8), np.int64(1))
Copy link
Contributor

Choose a reason for hiding this comment

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

This also might be improved with a tabular format, e.g. current -> np.uint8, new -> np.int64

more strict: This means that the results will be unchanged in the first
two examples. For the second one, the results will be the following::

np.arange(10, dtype=np.int8) + 256 # raises a TypeError
Copy link
Contributor

Choose a reason for hiding this comment

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

I think it's worth expanding on this example a bit too in the text below, just because this strikes me as the change which might be the most disruptive for users (in terms of breaking existing code that used to "work")

Copy link
Contributor

Choose a reason for hiding this comment

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

Agreed with @rossbar. This is likely to break a ton of essentially harmless code. I think it would be great to raise on integer overflow (consistently, not just with python scalars), but I think it would be vastly preferable to only do that when we actually fall out of (u)int64.

It's laudable to be more precise on this (i.e. if someone really wants to make sure their int8 stays an int8), but I think it's going to be a very sharp edge for the large majority of users who do not care about the underlying representation of their ints, but will be hit with errors much more often.

Comment on lines 327 to 438
As such implicit conversion to ``object`` should be rare and the work around
is clear, we expect that the backwards compatibility concerns are fairly small.
Copy link
Contributor

Choose a reason for hiding this comment

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

You could also mention that implicit conversion to object is generally something that we want to avoid; you could reference the ragged array nep.

Copy link
Contributor

@h-vetinari h-vetinari left a comment

Choose a reason for hiding this comment

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

I have a queasy feeling about this. I cannot speak about the benefits to improvements of implementation complexity of this proposal (though I don't doubt they are substantial), but from a user experience point it seems to be going into the wrong direction IMO.

If anything, I'd argue we should be fighting silent integer overflow tooth and nail, not introducing more surface where it might appear (e.g. by tending to default more to smaller integer types). I consider the breaking changes around this to be a non-negligible risk to numpy's reputation.

Moving away from the (value-based) promotion through the integer hierarchy is most likely similarly disadvantageous for the large majority of users who do not care about the machine representation of their integers.

I also tend to think that the distinction between 0-D arrays and scalars should be sharpened rather than weakened. IMO arrays should behave consistently (regardless of dimensionality), but the APIs between numpy arrays and python scalars are so different that silently ending up with the other one also poses several sharp edges.


np.array([100], dtype=np.uint8) + 200 == np.array([44], dtype=np.uint8)

Considering this, we believe that the proposal follows the "principle of least
Copy link
Contributor

Choose a reason for hiding this comment

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

I don't understand this, as the table above shows that this behaviour is maintained?

Silent integer overflow is a viciously surreptitious source of bugs, and if I read the NEP correctly, we're increasing the surface for that by defaulting to lower integer types - e.g. uint8(1) + 2: int64(3) (old) -> uint8(3) (new)

Comment on lines +91 to +106
.. figure:: _static/nep-0050-promotion-no-fonts.svg
:figclass: align-center
Copy link
Contributor

Choose a reason for hiding this comment

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

The horizontal lines in this graphic should have the same thickness and colour as the diagonal ones. At first glance it looks like uint8 can only promote to float16, not uint16.

grafik

Copy link
Member Author

Choose a reason for hiding this comment

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

Ah, got some alpha wrong or so, will fix, thanks.

Copy link
Contributor

Choose a reason for hiding this comment

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

The minimum Python scalar precision indicator is not very clear (both visually, and in concept)

Comment on lines 283 to 284
Will return an int64 array because the type of ``np.int64(1)`` is strictly
honoured.
Copy link
Contributor

Choose a reason for hiding this comment

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

It's not clear from context ("currently surprising cases") if this "will" is describing current or future behaviour.

more strict: This means that the results will be unchanged in the first
two examples. For the second one, the results will be the following::

np.arange(10, dtype=np.int8) + 256 # raises a TypeError
Copy link
Contributor

Choose a reason for hiding this comment

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

Agreed with @rossbar. This is likely to break a ton of essentially harmless code. I think it would be great to raise on integer overflow (consistently, not just with python scalars), but I think it would be vastly preferable to only do that when we actually fall out of (u)int64.

It's laudable to be more precise on this (i.e. if someone really wants to make sure their int8 stays an int8), but I think it's going to be a very sharp edge for the large majority of users who do not care about the underlying representation of their ints, but will be hit with errors much more often.

Comment on lines +319 to +348
Overflowing in the conversion rather than raising an error is a choice;
it is one that is the default in most C setups (similar to NumPy C can be
set up to raise an error due to the overflow, however).
Copy link
Contributor

Choose a reason for hiding this comment

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

I strongly believe that we should never silently overflow, but promote to the end of the integer hierarchy before doing that (so if there's ever a portable and universally available int128, we should keep going IMO).

I realise that this is incompatible with the proposed move away from value-based promotions, but to me that seems a step in the wrong direction (not talking about the internals and their complexity, but from a UX POV).

Comment on lines +348 to +378
This can lead to integer overflows and thus incorrect results beyond precision.
In many cases this may be silent, although NumPy usually gives warnings for the
scalar operators.
Copy link
Contributor

Choose a reason for hiding this comment

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

This sound like a really bad idea to me. I'm sorry to be so negative on this, but this could well lead to silent and catastrophic failures & data corruption, and IMO poses a substantial risk of blemishing numpy's reputation if not handled with extreme care.

Copy link
Member Author

Choose a reason for hiding this comment

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

Note that I am planning to add at least an optional warning tthat will notify of places where reduced precision occurs (I am not sure it would be 100%, but it should be good enough in practice).

Copy link
Contributor

Choose a reason for hiding this comment

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

Why will these warnings be optional? This is the case I am by far the most worried about.

Copy link
Member Author

Choose a reason for hiding this comment

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

Note that there are two warnings involved here:

  • Overflow warnings are normally given for NumPy scalar operations (which are the main affected ones here). (They are not given for array operations.)
  • Warnings about the actual change, which I was thinking should probably be opt-in.

For the second, I was aiming for opt-in currently, because I think they will be very noisy and cumbersome to avoid. E.g. for float32 >= 2.1 where the behaviour might change, but practically never does. Should we do the operation twice to check? Give a warning when it probably doesn't matter?

Maybe we could check whether always warnings should happen for more particularly surprising changes, right now I hoped the scalar integer overflow warnings cover the worst of it though.

* Certain operations are expected to start failing::

np.array([1], np.uint8) * 1000
np.array([1], np.uint8) == 1000 # possibly also
Copy link
Contributor

Choose a reason for hiding this comment

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

Why "possibly"? What does it depend on?

Copy link
Member Author

Choose a reason for hiding this comment

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

We could decide to return False, since we could upcast internally for comparisons. But you are right, the initial thing should fail, I think.

Comment on lines 379 to 382
np.multiple(float32_arr, 2.)
float32_arr * np.float64(2.)

Will both return a float64 rather than float32. This improves precision but
Copy link
Contributor

Choose a reason for hiding this comment

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

I thought python scalars would be treated as "weak"? Why does the first one upcast to float64?

Copy link
Member Author

Choose a reason for hiding this comment

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

I think I use "Python scalars" for int, float, complex. So the ones builtin with Python, not the NumPy scalars. But maybe I can try to avoid that term to avoid the tricky distinction.

Comment on lines +637 to +670
Such a distinction is very much possible, however, at this time NumPy will
often (and silently) convert 0-D arrays to scalars.
Copy link
Contributor

Choose a reason for hiding this comment

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

To me, it's more important that 0-D arrays behave consistently to their N-D counterparts than that they behave consistently with scalars. Some libraries such as CVXPY have tightened this distinction in recent times, and I think this is a healthy choice.

A numpy 0-D array has a completely different API than a python scalar, and liberally casting one to the other is not beneficial in my opinion. I think they should be compatible, but that types should be maintained, e.g. I expect a function operating on scalar values to maintain the input type (0-D array -> 0-D array; python scalar -> python scalar; without silently changing between them).

Copy link
Member Author

Choose a reason for hiding this comment

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

Yes, but NumPy does not currently have that choice, because we create scalars from 0-D arrays silently all the time. If scalars were distinct objects, you only got with a special operation. I might agree.

This NEP is only about the dtype of the array or scalar, not about this behaviour. I agree that it should be cleaned up, but it is another set of worms. But I also think that NumPy scalars have a dtype and because of that should behave the same as NumPy arrays for promotion.

Comment on lines +693 to +729
We reject this idea on the grounds that it will not remove the surprises
given earlier::

np.uint8(100) + 1000 == np.uint16(1100)
np.uint8(100) + 200 == np.uint8(44)
Copy link
Contributor

Choose a reason for hiding this comment

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

I'd argue that the first one is fine and the second should be np.uint16(300) (more details in the comments above).

@leofang
Copy link
Contributor

leofang commented Apr 15, 2022

cc: @kmaehashi @toslunar

@seberg
Copy link
Member Author

seberg commented Apr 20, 2022

@h-vetinari I have to catch up with this, I do not have a lot of time right now.

You seem to be mainly worried about integer overflows and the behaviour of that. The reality of it is that for array operations, we never check values. Most importantly we do not check the values of results. I do have an example somewhere that:

np.array([100], dtype=np.uint8) + 300  ->  uint16(400)
np.array([100], dtype=np.uint8) + 250  ->  uint8(94)

So we are increasing the amount of places where this happens, that is true. But the situation is much less useful currently than it may appear.
NumPy should not (and technically canont well) look at all values of two arrays. We also cannot really look at what the result would be to decide the output dtype. Our uint8 and int8 are not "integers" that are stored as int8, they are actual int8 in the C sense (adding two overflows in some way). They cannot be, since that would require to either just make the result int64 always (that does work but is not the reality [1]) or dynamically resize the dtype?

In practice, I do not think going the other way is a technically possible alternative at all. IMO NumPy has.

Should I add this as an additional sub-section to the (rejected) alternatives?


[1] In theory, we could have an int64[int8] on day. But the logic for it would be that it always returns an int64 on value-changing operations.

@h-vetinari
Copy link
Contributor

Most importantly we do not check the values of results. I do have an example somewhere that:

np.array([100], dtype=np.uint8) + 300  ->  uint16(400)
np.array([100], dtype=np.uint8) + 250  ->  uint8(94)

This to me is one of the most surprising and dangerous parts, and IMO an aspect of the design that needs to be fixed. Imagine that your code gets exercised during testing only with values falling into the first line, and once you're in production, you suddenly (and silently) hit the second one.

I think the idea to treat python scalars to only influence the type "weakly" is sound, but then it needs to take into account the value it is being added to to avoid silently causing overflow by choosing a too-small type for said python scalar in a really innocuous piece of code.

Our uint8 and int8 are not "integers" that are stored as int8, they are actual int8 in the C sense (adding two overflows in some way).

I understand that. My point is we should strive to protect users from the problems of C integers, rather than expose them even more.

We also cannot really look at what the result would be to decide the output dtype.

I think it's come up in various different shapes already, but would it make sense to expose different modes/errstates/etc. that allow users to explicitly opt into performance (no overflow check), while the default could be that we do take the performance hit to introspect enough to warn/error on overflow.

@shoyer
Copy link
Member

shoyer commented Apr 20, 2022

We also cannot really look at what the result would be to decide the output dtype.

I think it's come up in various different shapes already, but would it make sense to expose different modes/errstates/etc. that allow users to explicitly opt into performance (no overflow check), while the default could be that we do take the performance hit to introspect enough to warn/error on overflow.

I agree that support for explicitly checking for numerical overflow would be a great option to have, similar to how we handle dividing by zero (error, warn or ignore).

I'm not sure which default behavior I would prefer (maybe a warning would be a good compromise?), but any of these would be better than the current value dependent dtype promotion behavior.

@seberg seberg marked this pull request as ready for review May 24, 2022 14:23
@seberg
Copy link
Member Author

seberg commented May 24, 2022

I am sure there are still a lot of things left to discuss here and maybe change/reorganize. But should we aim to merge a draft very soon?

(I am also happy to do larger iterations here now, the question is whether the length/size of the review won't naturally stall process anyway soon.)

There are probably a few comments not addressed in depth. I will note that I did not address @h-vetinari comments about making things worse by not preventing integer overflows strictly.
I appreciate that integer overflows can be a huge user-trap, but I do not even know how the alternative would be possible to implement in a reasonable way (i.e. "improving" it rather than "making it worse").
The only realistic alternative that I see is already listed: The use of strong types for Python integers (which would prevent at least more integer overflows).

One day (maybe soon ;)), we could look into providing arbitrary Precision integers in NumPy (using tagged pointers), that "just" requires cleaning out how we do reference counting in object arrays currently...

@seberg
Copy link
Member Author

seberg commented May 24, 2022

The next larger iteration will be to create a branch (or PR) that doesn't just show the new behaviour, but also implements all the warnings I would like to have.
That may be necessary to really discuss how bad the break is for downstream.

Copy link
Contributor

@stefanv stefanv left a comment

Choose a reason for hiding this comment

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

Agreed that we should merge the draft NEP, and then start the approval process.

seberg and others added 19 commits May 25, 2022 07:24
Note, this really is a very early draft.  For myself, I have not
even 100% made up my mind that this is actually the best approach.

However, to solicitate feedback and input, I would like to put it
out there (and have a nice rendered version).

For those mainly interested in the possible different approaches
(i.e. making decisions), the "Alternatives" section is probably most
interesting.

Please make a note of blatend errors or unnecessary repitition, but
don't be too surprised if you see them, I will not be surprised if
this needs to be redone from scratch eventually.
Based on Tylers and Jake Vanderplas' comments.
…h object

This may be a bit "too much" information, but on the other hand it
is a backcompat break
…promotion

Even if I prefer not to mention the word "category" too much, integral vs. inexact
is the "magic" boundary when it comes to Python promoting an numpy
scalar to the higher category.
But I need to link to the old behaviour explenatation early on
somewhere...
And possibly it is now too much.  Plus, I am not sure the table
is actually all *that* useful
As per comment by Ross.  float64 can be put at the same precision
as "uint64".  This is incorrect, but it is the practical choice
that NumPy uses.
Co-authored-by: Stefan van der Walt <sjvdwalt@gmail.com>
@mattip mattip merged commit d985fd4 into numpy:main May 25, 2022
@mattip
Copy link
Member

mattip commented May 25, 2022

Putting this in not to end discussion but to provide a basis for further discussion

@seberg
Copy link
Member Author

seberg commented May 25, 2022

Indeed, it is just a basis for further discussion! @jakevdp, I am very curious whether there has been any further development/thoughts along this in the JAX team? Since I think you may have some important insights for NumPy as well.

@seberg seberg deleted the nep50 branch November 10, 2022 08:08
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.