Converting np.sinc into a ufunc

classic Classic list List threaded Threaded
12 messages Options
Reply | Threaded
Open this post in threaded view
|

Converting np.sinc into a ufunc

Sebastian Berg
Hi all,

there is an open PR (https://github.com/numpy/numpy/pull/12924) to
convert `np.sinc` into a ufunc. Since it should improve general
precision in `np.sinc`, I thought we could try to move that forward a
bit. We check whether this is worth it or not in the end.

However, it would also change behaviour slightly since `np.sinc(x=arr)`
will not work, as ufuncs are positional arguments only (we could wrap
`sinc`, but that hides all the nice features). Otherwise, there should
be no change except additional features of ufuncs and the move to a C
implementation.

This is mostly to see if anyone is worried about possible slight API
change here.

All the Best,

Sebastian

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpy-discussion

signature.asc (849 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Converting np.sinc into a ufunc

Nathan Goldbaum
It might be worth using BigQuery to search the github repository public dataset for usages of np.sinc with keyword arguments.

On Wed, May 22, 2019 at 1:05 PM Sebastian Berg <[hidden email]> wrote:
Hi all,

there is an open PR (https://github.com/numpy/numpy/pull/12924) to
convert `np.sinc` into a ufunc. Since it should improve general
precision in `np.sinc`, I thought we could try to move that forward a
bit. We check whether this is worth it or not in the end.

However, it would also change behaviour slightly since `np.sinc(x=arr)`
will not work, as ufuncs are positional arguments only (we could wrap
`sinc`, but that hides all the nice features). Otherwise, there should
be no change except additional features of ufuncs and the move to a C
implementation.

This is mostly to see if anyone is worried about possible slight API
change here.

All the Best,

Sebastian
_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpy-discussion

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Converting np.sinc into a ufunc

ralfgommers


On Wed, May 22, 2019 at 7:34 PM Nathan Goldbaum <[hidden email]> wrote:
It might be worth using BigQuery to search the github repository public dataset for usages of np.sinc with keyword arguments.

We spent some effort at Quansight to try different approaches to this. BigQuery turns out to be suboptimal, parsing code with ast.parse is more robust. Chris Ostrouchov just released some code for this (blog post with details to follow) and the results of running that code: https://github.com/Quansight-Labs/python-api-inspect/blob/master/data/numpy-summary.csv

np.sinc has 35 usages. to put that in perspective, np.array has ~31,000, np.dot ~2200, np.floor ~220, trace/inner/spacing/copyto are all similar to sinc.





On Wed, May 22, 2019 at 1:05 PM Sebastian Berg <[hidden email]> wrote:
Hi all,

there is an open PR (https://github.com/numpy/numpy/pull/12924) to
convert `np.sinc` into a ufunc. Since it should improve general
precision in `np.sinc`, I thought we could try to move that forward a
bit. We check whether this is worth it or not in the end.

Can you quantify the precision improvement (approximately)?


However, it would also change behaviour slightly since `np.sinc(x=arr)`
will not work, as ufuncs are positional arguments only (we could wrap
`sinc`, but that hides all the nice features).

This would give an exception, so at least it's easy to fix. As backwards compat breaks go, this is a pretty minor one I think.

Otherwise, there should
be no change except additional features of ufuncs and the move to a C
implementation.

I see this is one of the functions that uses asanyarray, so what about impact on subclass behavior?

Cheers,
Ralf


This is mostly to see if anyone is worried about possible slight API
change here.

All the Best,

Sebastian
_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpy-discussion
_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpy-discussion

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Converting np.sinc into a ufunc

Stephan Hoyer-2

On Wed, May 22, 2019 at 2:00 PM Ralf Gommers <[hidden email]> wrote:


On Wed, May 22, 2019 at 7:34 PM Nathan Goldbaum <[hidden email]> wrote:
It might be worth using BigQuery to search the github repository public dataset for usages of np.sinc with keyword arguments.

We spent some effort at Quansight to try different approaches to this. BigQuery turns out to be suboptimal, parsing code with ast.parse is more robust. Chris Ostrouchov just released some code for this (blog post with details to follow) and the results of running that code: https://github.com/Quansight-Labs/python-api-inspect/blob/master/data/numpy-summary.csv

np.sinc has 35 usages. to put that in perspective, np.array has ~31,000, np.dot ~2200, np.floor ~220, trace/inner/spacing/copyto are all similar to sinc.

Searching Google's internal code base (including open source dependencies), I found many uses of np.sinc, but no uses of the keyword argument "x".

I think it's pretty safe to go ahead here.

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Converting np.sinc into a ufunc

Marten van Kerkwijk
In reply to this post by ralfgommers

Otherwise, there should
be no change except additional features of ufuncs and the move to a C
implementation.

I see this is one of the functions that uses asanyarray, so what about impact on subclass behavior?

So, subclasses are passed on, as they are in ufuncs. In general, that should mean it is OK for subclasses. For astropy's Quantity, it would be an improvement, as `where` was never properly supported, and with a ufunc, we can handle it easily.

-- Marten

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Converting np.sinc into a ufunc

Marten van Kerkwijk
In reply to this post by Stephan Hoyer-2
On a more general note, if we change to a ufunc, it will get us stuck with sinc being the normalized version, where the units of the input have to be in the half-cycles preferred by signal-processing people rather than the radians preferred by mathematicians.

In this respect, note that there is an outstanding issue about whether to allow one to choose between the two: https://github.com/numpy/numpy/issues/13457 (which itself was raised following an inconclusive PR that tried to add a keyword argument for it).

Adding a keyword argument is much easier for a general function than for a ufunc.

-- Marten

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Converting np.sinc into a ufunc

Charles R Harris


On Wed, May 22, 2019 at 7:14 PM Marten van Kerkwijk <[hidden email]> wrote:
On a more general note, if we change to a ufunc, it will get us stuck with sinc being the normalized version, where the units of the input have to be in the half-cycles preferred by signal-processing people rather than the radians preferred by mathematicians.

In this respect, note that there is an outstanding issue about whether to allow one to choose between the two: https://github.com/numpy/numpy/issues/13457 (which itself was raised following an inconclusive PR that tried to add a keyword argument for it).

Adding a keyword argument is much easier for a general function than for a ufunc.


I'd be tempted to have two sinc functions with the different normalizations. Of course, one could say the same about trig functions in both radians and degrees. If I had to pick one, I'd choose sinc in radians, but I think that ship has sailed.

Chuck 

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Converting np.sinc into a ufunc

Joshua Wilson
Re Ralf's question:

> Can you quantify the precision improvement (approximately)?

On one level you'll get a large decrease in relative error around the
zeros of the sinc function because argument reduction is being done by
a number which is exactly representable in double precision (i.e. the
number 2) versus an irrational number. For example, consider:

>>> import numpy as np
>>> import mpmath
>>> x = 1 + 1e-8
>>> (np.sinc(x) - mpmath.sincpi(x))/mpmath.sincpi(x)

On master that will give you

mpf('-5.1753363184721223e-9')

versus

mpf('1.6543612517040003e-16')

on the new branch.

*But* there are some caveats to that answer: since we are close to the
zeros of the sinc function the condition number is large, so real
world data that has already been rounded before even calling sinc will
incur the same mathematically unavoidable loss of precision.

On Wed, May 22, 2019 at 6:24 PM Charles R Harris
<[hidden email]> wrote:

>
>
>
> On Wed, May 22, 2019 at 7:14 PM Marten van Kerkwijk <[hidden email]> wrote:
>>
>> On a more general note, if we change to a ufunc, it will get us stuck with sinc being the normalized version, where the units of the input have to be in the half-cycles preferred by signal-processing people rather than the radians preferred by mathematicians.
>>
>> In this respect, note that there is an outstanding issue about whether to allow one to choose between the two: https://github.com/numpy/numpy/issues/13457 (which itself was raised following an inconclusive PR that tried to add a keyword argument for it).
>>
>> Adding a keyword argument is much easier for a general function than for a ufunc.
>>
>
> I'd be tempted to have two sinc functions with the different normalizations. Of course, one could say the same about trig functions in both radians and degrees. If I had to pick one, I'd choose sinc in radians, but I think that ship has sailed.
>
> Chuck
> _______________________________________________
> NumPy-Discussion mailing list
> [hidden email]
> https://mail.python.org/mailman/listinfo/numpy-discussion
_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Converting np.sinc into a ufunc

ralfgommers
In reply to this post by Charles R Harris


On Thu, May 23, 2019 at 3:24 AM Charles R Harris <[hidden email]> wrote:


On Wed, May 22, 2019 at 7:14 PM Marten van Kerkwijk <[hidden email]> wrote:
On a more general note, if we change to a ufunc, it will get us stuck with sinc being the normalized version, where the units of the input have to be in the half-cycles preferred by signal-processing people rather than the radians preferred by mathematicians.

In this respect, note that there is an outstanding issue about whether to allow one to choose between the two: https://github.com/numpy/numpy/issues/13457 (which itself was raised following an inconclusive PR that tried to add a keyword argument for it).

Adding a keyword argument is much easier for a general function than for a ufunc.


I'd be tempted to have two sinc functions with the different normalizations. Of course, one could say the same about trig functions in both radians and degrees. If I had to pick one, I'd choose sinc in radians, but I think that ship has sailed.

Please let's not add more functions. This shouldn't be in numpy in the first place if we had that choice today (but the ship has sailed). I'd refer to scipy.special.sinc rather than expand the coverage here in ways that are specific to signal processing or some other domain.

Ralf



_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Converting np.sinc into a ufunc

Marten van Kerkwijk
I agree that we should not have two functions

I also am rather unsure whether a ufunc is a good idea. Earlier, while discussing other possible additions, like `erf`, the conclusion seemed to be that in numpy we should just cover whatever is in the C standard. This suggests `sinc` should not be a ufunc.

-- Marten

p.s.`scipy.special.sinc` *is* `np.sinc`

p.s.2 The accuracy argument is not in itself an argument for a ufunc, as it could be done in python too.

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Converting np.sinc into a ufunc

Robert Kern-2
On Thu, May 23, 2019 at 7:20 AM Marten van Kerkwijk <[hidden email]> wrote:
I agree that we should not have two functions

I also am rather unsure whether a ufunc is a good idea. Earlier, while discussing other possible additions, like `erf`, the conclusion seemed to be that in numpy we should just cover whatever is in the C standard. This suggests `sinc` should not be a ufunc.

That standard is for "what special functions we include in numpy, regardless of implementation", not "which special functions we already have in numpy should be implemented as ufuncs or regular functions".
 
--
Robert Kern

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Converting np.sinc into a ufunc

Sebastian Berg
In reply to this post by Marten van Kerkwijk
On Thu, 2019-05-23 at 10:17 -0400, Marten van Kerkwijk wrote:
> I agree that we should not have two functions
>
> I also am rather unsure whether a ufunc is a good idea. Earlier,
> while discussing other possible additions, like `erf`, the conclusion
> seemed to be that in numpy we should just cover whatever is in the C
> standard. This suggests `sinc` should not be a ufunc.
>

We are not adding a function though, as Robert already noted. So this
is about how much we prefer ufuncs for functions that are a perfect
fit.

The accuracy can certainly be improved, but it involves some branching,
so I think we would see at least 50% speed penalty compared to the
ufunc version in the end (right now the speed improvement is about 20%
for larger arrays, much more for smaller of course).

I do not have a perfect feeling about what the precision improvements
mean exactly here, but I posted some relative errors below as
additional stats [0].

Overall I think I would be pretty neutral if there was no gain at all
(as there is some loss of maintainability). Here it seems that we have
some decent enhancement as well though, so I am slightly in favor.

Best,

Sebastian


[0] And here maybe an additional point for better relative precision:
```
xarr = np.linspace(0, 2.1, 200000)
res = []
for x in xarr:
    res.append((np.sinc(x) - mpmath.sincpi(x))/mpmath.sincpi(x))

res = np.asarray(res, dtype=object)
print(abs(res).mean(), abs(res).max())
```

master:
1.70112344295248e-15 6.61977521930425e-11

New branch:
5.64884361112036e-17 4.01208410359583e-16

we probably should have tests if we add this.


> -- Marten
>
> p.s.`scipy.special.sinc` *is* `np.sinc`
>
> p.s.2 The accuracy argument is not in itself an argument for a ufunc,
> as it could be done in python too.
> _______________________________________________
> NumPy-Discussion mailing list
> [hidden email]
> https://mail.python.org/mailman/listinfo/numpy-discussion

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpy-discussion

signature.asc (849 bytes) Download Attachment