Skip to content

BUG: Make Polynomial evaluation adhere to nep 50 #26550

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 9 commits into from
Jul 6, 2024

Conversation

eendebakpt
Copy link
Contributor

  • Make sure evaluation of a Polynomial respects nep50 for scalar input:
import numpy as np

p=np.polynomial.Polynomial(np.array([1, 2], dtype=np.float32))
print(type(p(2))) # np.float64
# correct, since the domain argument is the default which maps to [-1., 1.]

w=np.array([-1,1], dtype=np.float32)
p=np.polynomial.Polynomial(np.array([1, 2], dtype=np.float32), domain=w, window=w)
print(type(p(2))) # np.float32 (was float64 on main)
import numpy as np
arr = np.polydiv(1, np.float32(1))
arr.dtype # float64

The input here are polynomial coefficients, which are really an array and not a scalar. So the output type seems correct, even though 1 looks like a scalar input.

Also see #26484

@eendebakpt eendebakpt changed the title Draft: BUG: Make Polynomial evaluation adhere to nep 50 BUG: Make Polynomial evaluation adhere to nep 50 May 27, 2024
@@ -346,7 +346,8 @@ def mapdomain(x, old, new):
array([-1.0+1.j , -0.6+0.6j, -0.2+0.2j, 0.2-0.2j, 0.6-0.6j, 1.0-1.j ]) # may vary

"""
x = np.asanyarray(x)
if not np.isscalar(x):
Copy link
Member

@seberg seberg May 27, 2024

Choose a reason for hiding this comment

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

Lets make this type(x) not in (int, float, complex), may be odder, but something we have in a few places, and np.isscalar is definitely terrible.

EDIT: Hmmm, isscalar isn't as bad as I thought, but it's still not really clearly defined what it does.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

That would exclude subclasses of int and float. We could maybe use instanceof(x, (int, float, complex))?

I will check later what the behaviour of other methods is. In NEP50 I could not find anything about subclasses.

Copy link
Member

Choose a reason for hiding this comment

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

TBH, it might be a mixed bag, with ufuncs letting them pass, we should check it. But subclasses are tedious to deal with and I am not sure we should, since the Python literals are really what is special.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Some facts:

  • ufuncs (tested np.atan2 and np.add) work with subclasses of a Python int according to nep50:
In [17]: import numpy as np
    ...: 
    ...: class MyInt(int):
    ...:     pass
    ...: 
    ...: x=np.float16(2.)
    ...: y=MyInt(1)
    ...: z=np.add(x, y)
    ...: print(z, type(z))
3.0 <class 'numpy.float16'>
  • For performance reasons isinstance is not that bad, np.isscalar makes a difference though:
In [13]: z = np.array([1., 2.])
    ...: %timeit np.isscalar(z)
    ...: %timeit type(z) in (int, float, complex)
    ...: %timeit isinstance(z, (int, float, complex))
417 ns ± 9.47 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)
52.8 ns ± 0.846 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)
65.1 ns ± 1.27 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)

In [14]: z = 4.
    ...: %timeit np.isscalar(z)
    ...: %timeit type(z) in (int, float, complex)
    ...: %timeit isinstance(z, (int, float, complex))
153 ns ± 52.3 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)
46.2 ns ± 1.57 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)
49.7 ns ± 1.01 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)

In [15]: z = MyInt(4.)
    ...: %timeit np.isscalar(z)
    ...: %timeit type(z) in (int, float, complex)
    ...: %timeit isinstance(z, (int, float, complex))
395 ns ± 5.07 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)
52.6 ns ± 0.689 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)
36.7 ns ± 1.19 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)

Three options to continue:

  • Leave as it. That would not be terrible, but there is an issue reported so apparently some people notice and care about it
  • Keep the current PR. This only handles pure Python int and float instances.
  • Use isinstance. Subclasses are handled, which seems more consistent with the ufuncs.

I have a small preference for the last option, but also think it is a corner case that is not very important in this case.

The choice deserves some thought as there might be more cases like this that will require more work to handle. If we adopt this PR (or reject) users might expect other cases to be handled similarly. (for example arr = np.polydiv(1, np.float32(1)))

Copy link
Member

Choose a reason for hiding this comment

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

Yeah, you are right and I am regretting a little I didn't change it again in ufuncs.

Over time, I am leaning more towards NEP 50 being really all about Python literals, and there are few reasons to replace those with subclasses.
In fact, most subclasses might not be "proper" floats or integers at all (proper in the sense of the substitution principle that you are really sure it is always OK to demote them to normal Pytho integers/floats).

So yes, ufuncs use that and I guess I should consider it a bug unless others lean strongly towards supporting subclasses :( (I appreciate you said you do).
But I suspect most "python side" places use the type() check right now (e.g. _array_converter, select, and quantile/percentile`).

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I put the topic on the agenda for the triage meeting.

Copy link
Member

Choose a reason for hiding this comment

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

We discussed this and prefer moving towards the strict type check unless anyone has a use-case for the opposite.
One thing to remember is that putting things into a NumPy array is invariant by nature and NEP 50 is really about making Python literals convenient and not so much about odd subclasses also.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Thanks for the feedback. I updated the PR. I left out bool as one of the native Python types to allow for (like in your suggestion). In this case including it would not change results as bool is converted to an array with dtype bool.

window : (2,) array_like, optional
Window, see `domain` for its use. The default value is [-1, 1].
Window, see `domain` for its use. The default value is [-1., 1.].
Copy link
Member

Choose a reason for hiding this comment

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

These are unrelated docstring fixes, right? (I don't mind, just checking)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

A tiny bit related: the window and domain type impacts the output type of the Polynomial evaluation.

@eendebakpt
Copy link
Contributor Author

@seberg I included a second fast path to address the performance regression from #26843. The change is on the same line as is changed in this PR, so it was logical to combine them. If needed I can split it in another PR.

@charris
Copy link
Member

charris commented Jul 5, 2024

close/reopen

@charris charris closed this Jul 5, 2024
@charris charris reopened this Jul 5, 2024
@seberg seberg merged commit a153fb2 into numpy:main Jul 6, 2024
119 of 124 checks passed
@seberg
Copy link
Member

seberg commented Jul 6, 2024

Thanks @eendebakpt that generic check looks fine. Wouldn't backport this, but it probably isn't enough of an issue to backport anyway.

eagunn pushed a commit to eagunn/numpy that referenced this pull request Jul 9, 2024
* Make sure evaluation of a Polynomial respects nep50 for scalar input:
```
import numpy as np

p=np.polynomial.Polynomial(np.array([1, 2], dtype=np.float32))
print(type(p(2))) # np.float64
# correct, since the domain argument is the default which maps to [-1., 1.]

w=np.array([-1,1], dtype=np.float32)
p=np.polynomial.Polynomial(np.array([1, 2], dtype=np.float32), domain=w, window=w)
print(type(p(2))) # np.float32 (was float64 on main)
```

* Update documentation of the various polynomial classes for the updated domain and window (was changed in numpy#24568)
* Not addressed:
```
import numpy as np
arr = np.polydiv(1, np.float32(1))
arr.dtype # float64
```
The input here are polynomial coefficients, which are really an array and not a scalar. So the output type seems correct, even though `1` looks like a scalar input.
ArvidJB pushed a commit to ArvidJB/numpy that referenced this pull request Nov 1, 2024
* Make sure evaluation of a Polynomial respects nep50 for scalar input:
```
import numpy as np

p=np.polynomial.Polynomial(np.array([1, 2], dtype=np.float32))
print(type(p(2))) # np.float64

w=np.array([-1,1], dtype=np.float32)
p=np.polynomial.Polynomial(np.array([1, 2], dtype=np.float32), domain=w, window=w)
print(type(p(2))) # np.float32 (was float64 on main)
```

* Update documentation of the various polynomial classes for the updated domain and window (was changed in numpy#24568)
* Not addressed:
```
import numpy as np
arr = np.polydiv(1, np.float32(1))
arr.dtype # float64
```
The input here are polynomial coefficients, which are really an array and not a scalar. So the output type seems correct, even though `1` looks like a scalar input.
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.

3 participants