nanmean(), nanstd() and other "missing" functions for 1.8
Currently, I am in the process of migrating some co-workers from Matlab and IDL, and the number one complaint I get is that numpy has nansum() but no nanmean() and nanstd(). While we do have an alternative in the form of masked arrays, most of these people are busy enough trying to port their existing code over to python that this sort of stumbling block becomes difficult to explain away. Given how relatively simple these functions are, I can not think of any reason not to include these functions in v1.8. Of course, the documentation for these functions should certainly include mention of masked arrays. There is one other non-trivial function that have been discussed before: np.minmax(). My thinking is that it would return a 2xN array (where N is whatever size of the result that would be returned if just np.min() was used). This would allow one to do "min, max = np.minmax(X)". Are there any other functions that others feel are "missing" from numpy and would like to see for v1.8? Let's discuss them here. Cheers! Ben Root
On 1 May 2013 03:36, Benjamin Root <ben.root@ou.edu> wrote:
There is one other non-trivial function that have been discussed before: np.minmax(). My thinking is that it would return a 2xN array (where N is whatever size of the result that would be returned if just np.min() was used). This would allow one to do "min, max = np.minmax(X)".
I had been looking for this function in the past, I think it is a good and necessary addition. It should also come with its companion, np.argminmax. David.
could anyone on 32bit system with fresh numpy (1.7.1) test following:
wget -nc http://www.onerussian.com/tmp/data.npy ; python -c 'import numpy as np; data1 = np.load("/tmp/data.npy"); print np.sum(data1[1,:,0,1]) - np.sum(data1, axis=1)[1,0,1]'
0.0 because unfortunately it seems on fresh ubuntu raring (in 32bit build only, seems ok in 64 bit... also never ran into it on older numpy releases):
python -c 'import numpy as np; data1 = np.load("/tmp/data.npy"); print np.sum(data1[1,:,0,1]) - np.sum(data1, axis=1)[1,0,1]' -1.11022302463e-16
PS detected by failed tests of pymvpa -- Yaroslav O. Halchenko, Ph.D. http://neuro.debian.net http://www.pymvpa.org http://www.fail2ban.org Senior Research Associate, Psychological and Brain Sciences Dept. Dartmouth College, 419 Moore Hall, Hinman Box 6207, Hanover, NH 03755 Phone: +1 (603) 646-9834 Fax: +1 (603) 646-1419 WWW: http://www.linkedin.com/in/yarik
Hi, On Tue, Apr 30, 2013 at 8:08 PM, Yaroslav Halchenko <lists@onerussian.com> wrote:
could anyone on 32bit system with fresh numpy (1.7.1) test following:
wget -nc http://www.onerussian.com/tmp/data.npy ; python -c 'import numpy as np; data1 = np.load("/tmp/data.npy"); print np.sum(data1[1,:,0,1]) - np.sum(data1, axis=1)[1,0,1]'
0.0
because unfortunately it seems on fresh ubuntu raring (in 32bit build only, seems ok in 64 bit... also never ran into it on older numpy releases):
python -c 'import numpy as np; data1 = np.load("/tmp/data.npy"); print np.sum(data1[1,:,0,1]) - np.sum(data1, axis=1)[1,0,1]' -1.11022302463e-16
PS detected by failed tests of pymvpa
Reduced case on numpy 1.7.1, 32-bit Ubuntu 12.04.2 In [64]: data = np.array([[ 0.49505185, 0.47212842], [ 0.53529587, 0.04366172], [-0.13461665, -0.01664215]]) In [65]: np.sum(data[:, 0]) - np.sum(data, axis=0)[0] Out[65]: 1.1102230246251565e-16 No difference for single vector: In [4]: data1 = data[:, 0:1] In [5]: np.sum(data1[:, 0]) - np.sum(data1, axis=0)[0] Out[5]: 0.0 Puzzling to me... Cheers, Matthew
Hi, On Tue, Apr 30, 2013 at 9:16 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Tue, Apr 30, 2013 at 8:08 PM, Yaroslav Halchenko <lists@onerussian.com> wrote:
could anyone on 32bit system with fresh numpy (1.7.1) test following:
wget -nc http://www.onerussian.com/tmp/data.npy ; python -c 'import numpy as np; data1 = np.load("/tmp/data.npy"); print np.sum(data1[1,:,0,1]) - np.sum(data1, axis=1)[1,0,1]'
0.0
because unfortunately it seems on fresh ubuntu raring (in 32bit build only, seems ok in 64 bit... also never ran into it on older numpy releases):
python -c 'import numpy as np; data1 = np.load("/tmp/data.npy"); print np.sum(data1[1,:,0,1]) - np.sum(data1, axis=1)[1,0,1]' -1.11022302463e-16
PS detected by failed tests of pymvpa
Reduced case on numpy 1.7.1, 32-bit Ubuntu 12.04.2
In [64]: data = np.array([[ 0.49505185, 0.47212842], [ 0.53529587, 0.04366172], [-0.13461665, -0.01664215]])
In [65]: np.sum(data[:, 0]) - np.sum(data, axis=0)[0] Out[65]: 1.1102230246251565e-16
No difference for single vector:
In [4]: data1 = data[:, 0:1]
In [5]: np.sum(data1[:, 0]) - np.sum(data1, axis=0)[0] Out[5]: 0.0
Also true on current numpy trunk: In [2]: import numpy as np In [3]: np.__version__ Out[3]: '1.8.0.dev-a8805f6' In [4]: data = np.array([[ 0.49505185, 0.47212842], ....: [ 0.53529587, 0.04366172], ....: [-0.13461665, -0.01664215]]) In [5]: np.sum(data[:, 0]) - np.sum(data, axis=0)[0] Out[5]: 1.1102230246251565e-16 Not true on numpy 1.6.1: In [2]: np.__version__ Out[2]: '1.6.1' In [3]: data = np.array([[ 0.49505185, 0.47212842], ....: [ 0.53529587, 0.04366172], ....: [-0.13461665, -0.01664215]]) In [4]: np.sum(data[:, 0]) - np.sum(data, axis=0)[0] Out[4]: 0.0 Cheers, Matthew
On Tue, 2013-04-30 at 22:20 -0700, Matthew Brett wrote:
Hi,
On Tue, Apr 30, 2013 at 9:16 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Tue, Apr 30, 2013 at 8:08 PM, Yaroslav Halchenko <lists@onerussian.com> wrote:
could anyone on 32bit system with fresh numpy (1.7.1) test following:
wget -nc http://www.onerussian.com/tmp/data.npy ; python -c 'import numpy as np; data1 = np.load("/tmp/data.npy"); print np.sum(data1[1,:,0,1]) - np.sum(data1, axis=1)[1,0,1]'
0.0
because unfortunately it seems on fresh ubuntu raring (in 32bit build only, seems ok in 64 bit... also never ran into it on older numpy releases):
python -c 'import numpy as np; data1 = np.load("/tmp/data.npy"); print np.sum(data1[1,:,0,1]) - np.sum(data1, axis=1)[1,0,1]' -1.11022302463e-16
PS detected by failed tests of pymvpa
Reduced case on numpy 1.7.1, 32-bit Ubuntu 12.04.2
In [64]: data = np.array([[ 0.49505185, 0.47212842], [ 0.53529587, 0.04366172], [-0.13461665, -0.01664215]])
In [65]: np.sum(data[:, 0]) - np.sum(data, axis=0)[0] Out[65]: 1.1102230246251565e-16
No difference for single vector:
In [4]: data1 = data[:, 0:1]
In [5]: np.sum(data1[:, 0]) - np.sum(data1, axis=0)[0] Out[5]: 0.0
Also true on current numpy trunk:
In [2]: import numpy as np
In [3]: np.__version__ Out[3]: '1.8.0.dev-a8805f6'
In [4]: data = np.array([[ 0.49505185, 0.47212842], ....: [ 0.53529587, 0.04366172], ....: [-0.13461665, -0.01664215]])
In [5]: np.sum(data[:, 0]) - np.sum(data, axis=0)[0] Out[5]: 1.1102230246251565e-16
Not true on numpy 1.6.1:
In [2]: np.__version__ Out[2]: '1.6.1'
In [3]: data = np.array([[ 0.49505185, 0.47212842], ....: [ 0.53529587, 0.04366172], ....: [-0.13461665, -0.01664215]])
In [4]: np.sum(data[:, 0]) - np.sum(data, axis=0)[0] Out[4]: 0.0
Puzzles me, I didn't think calculation order was different in both cases and optimized for the reduction part. But maybe check the code, if it is optimized, it would calculate this more like `res += data[0]; res += data[1]; res += data[2]` (for faster memory access), which would probably kill the extended registers (I don't know this hardware stuff, so might be wrong). One simple try hinting that this may be going on would be to data fortran order. - Sebastian
Cheers,
Matthew _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On Wed, 2013-05-01 at 11:00 +0200, Sebastian Berg wrote:
On Tue, 2013-04-30 at 22:20 -0700, Matthew Brett wrote:
Hi,
On Tue, Apr 30, 2013 at 9:16 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Tue, Apr 30, 2013 at 8:08 PM, Yaroslav Halchenko <lists@onerussian.com> wrote:
could anyone on 32bit system with fresh numpy (1.7.1) test following:
wget -nc http://www.onerussian.com/tmp/data.npy ; python -c 'import numpy as np; data1 = np.load("/tmp/data.npy"); print np.sum(data1[1,:,0,1]) - np.sum(data1, axis=1)[1,0,1]'
0.0
because unfortunately it seems on fresh ubuntu raring (in 32bit build only, seems ok in 64 bit... also never ran into it on older numpy releases):
python -c 'import numpy as np; data1 = np.load("/tmp/data.npy"); print np.sum(data1[1,:,0,1]) - np.sum(data1, axis=1)[1,0,1]' -1.11022302463e-16
PS detected by failed tests of pymvpa
Reduced case on numpy 1.7.1, 32-bit Ubuntu 12.04.2
In [64]: data = np.array([[ 0.49505185, 0.47212842], [ 0.53529587, 0.04366172], [-0.13461665, -0.01664215]])
In [65]: np.sum(data[:, 0]) - np.sum(data, axis=0)[0] Out[65]: 1.1102230246251565e-16
No difference for single vector:
In [4]: data1 = data[:, 0:1]
In [5]: np.sum(data1[:, 0]) - np.sum(data1, axis=0)[0] Out[5]: 0.0
Also true on current numpy trunk:
In [2]: import numpy as np
In [3]: np.__version__ Out[3]: '1.8.0.dev-a8805f6'
In [4]: data = np.array([[ 0.49505185, 0.47212842], ....: [ 0.53529587, 0.04366172], ....: [-0.13461665, -0.01664215]])
In [5]: np.sum(data[:, 0]) - np.sum(data, axis=0)[0] Out[5]: 1.1102230246251565e-16
Not true on numpy 1.6.1:
In [2]: np.__version__ Out[2]: '1.6.1'
In [3]: data = np.array([[ 0.49505185, 0.47212842], ....: [ 0.53529587, 0.04366172], ....: [-0.13461665, -0.01664215]])
In [4]: np.sum(data[:, 0]) - np.sum(data, axis=0)[0] Out[4]: 0.0
Puzzles me, I didn't think calculation order was different in both cases and optimized for the reduction part. But maybe check the code, if it is optimized, it would calculate this more like `res += data[0]; res += data[1]; res += data[2]` (for faster memory access), which would probably kill the extended registers (I don't know this hardware stuff, so might be wrong). One simple try hinting that this may be going on would be to data fortran order.
Well I guess it is optimized and the reason I thought it was not, was because it probably was not before 1.7. so I think this is probably the reason.
- Sebastian
Cheers,
Matthew _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Thanks everyone for the feedback. Is it worth me starting a bisection to catch where it was introduced? On Wed, 01 May 2013, Sebastian Berg wrote:
so might be wrong). One simple try hinting that this may be going on would be to data fortran order.
Well I guess it is optimized and the reason I thought it was not, was because it probably was not before 1.7. so I think this is probably the reason.
-- Yaroslav O. Halchenko, Ph.D. http://neuro.debian.net http://www.pymvpa.org http://www.fail2ban.org Senior Research Associate, Psychological and Brain Sciences Dept. Dartmouth College, 419 Moore Hall, Hinman Box 6207, Hanover, NH 03755 Phone: +1 (603) 646-9834 Fax: +1 (603) 646-1419 WWW: http://www.linkedin.com/in/yarik
On 1 May 2013 08:49, "Yaroslav Halchenko" <lists@onerussian.com> wrote:
Thanks everyone for the feedback.
Is it worth me starting a bisection to catch where it was introduced?
Is it a bug, or just typical fp rounding issues? Do we know which answer is correct? -n
On Wed, 01 May 2013, Nathaniel Smith wrote:
Thanks everyone for the feedback. Is it worth me starting a bisection to catch where it was introduced?
Is it a bug, or just typical fp rounding issues? Do we know which answer is correct?
to ignorant me, even without considering 'correctness', it is just a typical regression -- results changed from one release to another (and not to the better side). -- Yaroslav O. Halchenko, Ph.D. http://neuro.debian.net http://www.pymvpa.org http://www.fail2ban.org Senior Research Associate, Psychological and Brain Sciences Dept. Dartmouth College, 419 Moore Hall, Hinman Box 6207, Hanover, NH 03755 Phone: +1 (603) 646-9834 Fax: +1 (603) 646-1419 WWW: http://www.linkedin.com/in/yarik
01.05.2013 16:01, Yaroslav Halchenko kirjoitti: [clip]
to ignorant me, even without considering 'correctness', it is just a typical regression -- results changed from one release to another (and not to the better side).
To me this seems to be a consequence of performing additions in a different order than previously. Both results are IMHO correct, so I'm not sure there is anything to fix here. Third-party code relying on a certain outcome of rounding error is likely incorrect anyway. -- Pauli Virtanen
On Wed, May 1, 2013 at 9:12 AM, Pauli Virtanen <pav@iki.fi> wrote:
01.05.2013 16:01, Yaroslav Halchenko kirjoitti: [clip]
to ignorant me, even without considering 'correctness', it is just a typical regression -- results changed from one release to another (and not to the better side).
To me this seems to be a consequence of performing additions in a different order than previously. Both results are IMHO correct, so I'm not sure there is anything to fix here. Third-party code relying on a certain outcome of rounding error is likely incorrect anyway.
Yeah, seems to just be the standard floating point indeterminism. Using Matthew's numbers and pure Python floats: In [9]: (0.49505185 + 0.53529587) + -0.13461665 Out[9]: 0.89573107 In [10]: 0.49505185 + (0.53529587 + -0.13461665) Out[10]: 0.8957310700000001 In [11]: _9 - _10 Out[11]: -1.1102230246251565e-16 Looks like a bug in pymvpa or its test suite to me. -n
On Wed, 01 May 2013, Nathaniel Smith wrote:
not sure there is anything to fix here. Third-party code relying on a certain outcome of rounding error is likely incorrect anyway.
Yeah, seems to just be the standard floating point indeterminism. Using Matthew's numbers and pure Python floats:
In [9]: (0.49505185 + 0.53529587) + -0.13461665 Out[9]: 0.89573107
In [10]: 0.49505185 + (0.53529587 + -0.13461665) Out[10]: 0.8957310700000001
In [11]: _9 - _10 Out[11]: -1.1102230246251565e-16
Looks like a bug in pymvpa or its test suite to me.
well -- sure thing we will "fix" the unittest to not rely on precise correspondence any longer since released 1.7.1 is effected. So it is not a matter of me avoiding "fixing" pymvpa's "bug". I brought it to your attention because 1. from e.g. np.sum(data[:, 0]) - np.sum(data, axis=0)[0] which presumably should be the same order of additions for 0-th column it is not clear that they do not have to be identical. 2. so far they were identical across many numpy releases 3. they are identical on other architectures (e.g. amd64) -- Yaroslav O. Halchenko, Ph.D. http://neuro.debian.net http://www.pymvpa.org http://www.fail2ban.org Senior Research Associate, Psychological and Brain Sciences Dept. Dartmouth College, 419 Moore Hall, Hinman Box 6207, Hanover, NH 03755 Phone: +1 (603) 646-9834 Fax: +1 (603) 646-1419 WWW: http://www.linkedin.com/in/yarik
HI, On Wed, May 1, 2013 at 9:09 AM, Yaroslav Halchenko <lists@onerussian.com> wrote:
On Wed, 01 May 2013, Nathaniel Smith wrote:
not sure there is anything to fix here. Third-party code relying on a certain outcome of rounding error is likely incorrect anyway.
Yeah, seems to just be the standard floating point indeterminism. Using Matthew's numbers and pure Python floats:
In [9]: (0.49505185 + 0.53529587) + -0.13461665 Out[9]: 0.89573107
In [10]: 0.49505185 + (0.53529587 + -0.13461665) Out[10]: 0.8957310700000001
In [11]: _9 - _10 Out[11]: -1.1102230246251565e-16
Looks like a bug in pymvpa or its test suite to me.
well -- sure thing we will "fix" the unittest to not rely on precise correspondence any longer since released 1.7.1 is effected. So it is not a matter of me avoiding "fixing" pymvpa's "bug".
I brought it to your attention because
1. from e.g.
np.sum(data[:, 0]) - np.sum(data, axis=0)[0]
which presumably should be the same order of additions for 0-th column it is not clear that they do not have to be identical.
I agree it's surprising, but I guess it's reasonable for numpy to reserve the right to add these guys up in whatever order it chooses, and (in this case) maybe a different order for the axis=None, axis=X cases. Also, y'all may have noticed that it is the presence of the second vector in the array which causes the difference in the sums of the first (see my first email in this thread). If this is an order effect I guess this means that the order of operations in an sum(a, axis=X) operation depends on the shape of the array. And it looks like it depends on memory layout: In [24]: data = np.array([[ 0.49505185, 0], ....: [ 0.53529587, 0], ....: [-0.13461665, 0]]) In [25]: np.sum(data[:, 0]) - np.sum(data, axis=0)[0] Out[25]: 1.1102230246251565e-16 In [26]: data_F = np.array(data, order='F') In [27]: np.sum(data_F[:, 0]) - np.sum(data_F, axis=0)[0] Out[27]: 0.0 Do we allow the results to be different for different memory layout?
2. so far they were identical across many numpy releases
3. they are identical on other architectures (e.g. amd64)
To me that is surprising. I would have guessed that the order is the same on 32 and 64 bit, but something about the precision of intermediate operations is different. I don't know enough about amd64 to guess what that could be. Bradley's suggestion seems kind of reasonable but it's strange that numpy should use intel-80 bit intermediate values differently for 32 and 64 bit. Cheers, Matthew
On Wed, May 1, 2013 at 6:24 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
HI,
On Wed, May 1, 2013 at 9:09 AM, Yaroslav Halchenko <lists@onerussian.com> wrote:
3. they are identical on other architectures (e.g. amd64)
To me that is surprising. I would have guessed that the order is the same on 32 and 64 bit, but something about the precision of intermediate operations is different. I don't know enough about amd64 to guess what that could be. Bradley's suggestion seems kind of reasonable but it's strange that numpy should use intel-80 bit intermediate values differently for 32 and 64 bit.
"numpy" isn't doing anything different between the two. numpy generates the same C code. The C compiler may be generating different machine instructions for that code on different architectures, even closely related ones like i386 and amd64. Different optimization flags and compiler versions will probably also affect this, not just the target architecture. It's possible that those are actually the source of this observation. -- Robert Kern
just for completeness... I haven't yet double checked if I have done it correctly but here is the bisected commit: aed9925a9d5fe9a407d0ca2c65cb577116c4d0f1 is the first bad commit commit aed9925a9d5fe9a407d0ca2c65cb577116c4d0f1 Author: Mark Wiebe <mwiebe@enthought.com> Date: Tue Aug 2 13:34:13 2011 -0500 ENH: ufunc: Rewrite PyUFunc_Reduce to be more general and easier to adapt to NA masks This generalizes the 'axis' parameter to accept None or a list of axes on which to do the reduction. :040000 040000 2bdd71a1ea60c0dbfe370c77f69724fab28038e1 44f54a15f480ccaf519d10e9c42032de86bd0dca M numpy bisect run success FWIW ( ;-) ): # git describe --tags aed9925a9d5fe9a407d0ca2c65cb577116c4d0f1 v0.3.0-7757-gaed9925 # git describe --tags --contains aed9925a9d5fe9a407d0ca2c65cb577116c4d0f1 v1.7.0b1~377^2~126 On Wed, 01 May 2013, Robert Kern wrote:
On Wed, May 1, 2013 at 6:24 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
HI,
On Wed, May 1, 2013 at 9:09 AM, Yaroslav Halchenko <lists@onerussian.com> wrote:
3. they are identical on other architectures (e.g. amd64)
To me that is surprising. I would have guessed that the order is the same on 32 and 64 bit, but something about the precision of intermediate operations is different. I don't know enough about amd64 to guess what that could be. Bradley's suggestion seems kind of reasonable but it's strange that numpy should use intel-80 bit intermediate values differently for 32 and 64 bit.
"numpy" isn't doing anything different between the two. numpy generates the same C code. The C compiler may be generating different machine instructions for that code on different architectures, even closely related ones like i386 and amd64. Different optimization flags and compiler versions will probably also affect this, not just the target architecture. It's possible that those are actually the source of this observation. -- Yaroslav O. Halchenko, Ph.D. http://neuro.debian.net http://www.pymvpa.org http://www.fail2ban.org Senior Research Associate, Psychological and Brain Sciences Dept. Dartmouth College, 419 Moore Hall, Hinman Box 6207, Hanover, NH 03755 Phone: +1 (603) 646-9834 Fax: +1 (603) 646-1419 WWW: http://www.linkedin.com/in/yarik
On Wed, 2013-05-01 at 15:29 -0400, Yaroslav Halchenko wrote:
just for completeness... I haven't yet double checked if I have done it correctly but here is the bisected commit:
aed9925a9d5fe9a407d0ca2c65cb577116c4d0f1 is the first bad commit commit aed9925a9d5fe9a407d0ca2c65cb577116c4d0f1 Author: Mark Wiebe <mwiebe@enthought.com> Date: Tue Aug 2 13:34:13 2011 -0500
ENH: ufunc: Rewrite PyUFunc_Reduce to be more general and easier to adapt to NA masks
This generalizes the 'axis' parameter to accept None or a list of axes on which to do the reduction.
:040000 040000 2bdd71a1ea60c0dbfe370c77f69724fab28038e1 44f54a15f480ccaf519d10e9c42032de86bd0dca M numpy bisect run success
FWIW ( ;-) ):
There really is no point discussing here, this has to do with numpy doing iteration order optimization, and you actually *want* this. Lets for a second assume that the old behavior was better, then the next guy is going to ask: "Why is np.add.reduce(array, axis=0) so much slower then reduce(array, np.add)?". This is huge speed improvement by Marks new iterator for reductions over the slow axes, so instead of trying to track "regressions" down, I think the right thing is to say kudos for doing this improvement :). Just my opinion, Sebastian
# git describe --tags aed9925a9d5fe9a407d0ca2c65cb577116c4d0f1 v0.3.0-7757-gaed9925
# git describe --tags --contains aed9925a9d5fe9a407d0ca2c65cb577116c4d0f1 v1.7.0b1~377^2~126
On Wed, 01 May 2013, Robert Kern wrote:
On Wed, May 1, 2013 at 6:24 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
HI,
On Wed, May 1, 2013 at 9:09 AM, Yaroslav Halchenko <lists@onerussian.com> wrote:
3. they are identical on other architectures (e.g. amd64)
To me that is surprising. I would have guessed that the order is the same on 32 and 64 bit, but something about the precision of intermediate operations is different. I don't know enough about amd64 to guess what that could be. Bradley's suggestion seems kind of reasonable but it's strange that numpy should use intel-80 bit intermediate values differently for 32 and 64 bit.
"numpy" isn't doing anything different between the two. numpy generates the same C code. The C compiler may be generating different machine instructions for that code on different architectures, even closely related ones like i386 and amd64. Different optimization flags and compiler versions will probably also affect this, not just the target architecture. It's possible that those are actually the source of this observation.
Hi, On Wed, May 1, 2013 at 1:01 PM, Sebastian Berg <sebastian@sipsolutions.net> wrote:
On Wed, 2013-05-01 at 15:29 -0400, Yaroslav Halchenko wrote:
just for completeness... I haven't yet double checked if I have done it correctly but here is the bisected commit:
aed9925a9d5fe9a407d0ca2c65cb577116c4d0f1 is the first bad commit commit aed9925a9d5fe9a407d0ca2c65cb577116c4d0f1 Author: Mark Wiebe <mwiebe@enthought.com> Date: Tue Aug 2 13:34:13 2011 -0500
ENH: ufunc: Rewrite PyUFunc_Reduce to be more general and easier to adapt to NA masks
This generalizes the 'axis' parameter to accept None or a list of axes on which to do the reduction.
:040000 040000 2bdd71a1ea60c0dbfe370c77f69724fab28038e1 44f54a15f480ccaf519d10e9c42032de86bd0dca M numpy bisect run success
FWIW ( ;-) ):
There really is no point discussing here, this has to do with numpy doing iteration order optimization, and you actually *want* this. Lets for a second assume that the old behavior was better, then the next guy is going to ask: "Why is np.add.reduce(array, axis=0) so much slower then reduce(array, np.add)?". This is huge speed improvement by Marks new iterator for reductions over the slow axes, so instead of trying to track "regressions" down, I think the right thing is to say kudos for doing this improvement :).
I don't believe Yarick meant his bisection to be a criticism, but as an aid to full understanding. Is it an issue that Fortran and C contiguous arrays give different rounding error for the sums? Cheers, Matthew
On Wed, 01 May 2013, Matthew Brett wrote:
There really is no point discussing here, this has to do with numpy doing iteration order optimization, and you actually *want* this. Lets for a second assume that the old behavior was better, then the next guy is going to ask: "Why is np.add.reduce(array, axis=0) so much slower then reduce(array, np.add)?". This is huge speed improvement by Marks new iterator for reductions over the slow axes, so instead of trying to track "regressions" down, I think the right thing is to say kudos for doing this improvement :).
I don't believe Yarick meant his bisection to be a criticism, but as an aid to full understanding.
Exactly right, Matthew -- thank you! And kudos to Mark! N.B. I am generally furry and kind, not fuzzy and evil -- Yaroslav O. Halchenko, Ph.D. http://neuro.debian.net http://www.pymvpa.org http://www.fail2ban.org Senior Research Associate, Psychological and Brain Sciences Dept. Dartmouth College, 419 Moore Hall, Hinman Box 6207, Hanover, NH 03755 Phone: +1 (603) 646-9834 Fax: +1 (603) 646-1419 WWW: http://www.linkedin.com/in/yarik
On Wed, 01 May 2013, Sebastian Berg wrote:
There really is no point discussing here, this has to do with numpy doing iteration order optimization, and you actually *want* this. Lets for a second assume that the old behavior was better, then the next guy is going to ask: "Why is np.add.reduce(array, axis=0) so much slower then reduce(array, np.add)?". This is huge speed improvement by Marks new iterator for reductions over the slow axes...
btw -- is there something like panda's vbench for numpy? i.e. where it would be possible to track/visualize such performance improvements/hits? -- Yaroslav O. Halchenko, Ph.D. http://neuro.debian.net http://www.pymvpa.org http://www.fail2ban.org Senior Research Associate, Psychological and Brain Sciences Dept. Dartmouth College, 419 Moore Hall, Hinman Box 6207, Hanover, NH 03755 Phone: +1 (603) 646-9834 Fax: +1 (603) 646-1419 WWW: http://www.linkedin.com/in/yarik
On Wed, 2013-05-01 at 16:37 -0400, Yaroslav Halchenko wrote:
On Wed, 01 May 2013, Sebastian Berg wrote:
There really is no point discussing here, this has to do with numpy doing iteration order optimization, and you actually *want* this. Lets for a second assume that the old behavior was better, then the next guy is going to ask: "Why is np.add.reduce(array, axis=0) so much slower then reduce(array, np.add)?". This is huge speed improvement by Marks new iterator for reductions over the slow axes...
btw -- is there something like panda's vbench for numpy? i.e. where it would be possible to track/visualize such performance improvements/hits?
Sorry if it seemed harsh, but only skimmed mails and it seemed a bit like the an obvious piece was missing... There are no benchmark tests I am aware of. You can try: a = np.random.random((1000, 1000)) and then time a.sum(1) and a.sum(0), on 1.7. the fast axis (1), is only slightly faster then the sum over the slow axis. On earlier numpy versions you will probably see something like half the speed for the slow axis (only got ancient or 1.7 numpy right now, so reluctant to give exact timings). - Sebastian
On Wed, May 1, 2013 at 1:24 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
HI,
On Wed, May 1, 2013 at 9:09 AM, Yaroslav Halchenko <lists@onerussian.com> wrote:
On Wed, 01 May 2013, Nathaniel Smith wrote:
not sure there is anything to fix here. Third-party code relying on a certain outcome of rounding error is likely incorrect anyway.
Yeah, seems to just be the standard floating point indeterminism. Using Matthew's numbers and pure Python floats:
In [9]: (0.49505185 + 0.53529587) + -0.13461665 Out[9]: 0.89573107
In [10]: 0.49505185 + (0.53529587 + -0.13461665) Out[10]: 0.8957310700000001
In [11]: _9 - _10 Out[11]: -1.1102230246251565e-16
Looks like a bug in pymvpa or its test suite to me.
well -- sure thing we will "fix" the unittest to not rely on precise correspondence any longer since released 1.7.1 is effected. So it is not a matter of me avoiding "fixing" pymvpa's "bug".
I brought it to your attention because
1. from e.g.
np.sum(data[:, 0]) - np.sum(data, axis=0)[0]
which presumably should be the same order of additions for 0-th column it is not clear that they do not have to be identical.
I agree it's surprising, but I guess it's reasonable for numpy to reserve the right to add these guys up in whatever order it chooses, and (in this case) maybe a different order for the axis=None, axis=X cases.
Also, y'all may have noticed that it is the presence of the second vector in the array which causes the difference in the sums of the first (see my first email in this thread). If this is an order effect I guess this means that the order of operations in an sum(a, axis=X) operation depends on the shape of the array. And it looks like it depends on memory layout:
In [24]: data = np.array([[ 0.49505185, 0], ....: [ 0.53529587, 0], ....: [-0.13461665, 0]])
In [25]: np.sum(data[:, 0]) - np.sum(data, axis=0)[0] Out[25]: 1.1102230246251565e-16
In [26]: data_F = np.array(data, order='F')
In [27]: np.sum(data_F[:, 0]) - np.sum(data_F, axis=0)[0] Out[27]: 0.0
Do we allow the results to be different for different memory layout?
Wasn't this the point of some of Mark Wiebe's optimization? As far as I understand he got rid of some of the C bias, and made calculations faster for Fortran contiguous arrays. I rather have speed for my Fortran arrays, then relying on float precision issues, that bite anyway when running on many different kinds of machines. (I usually have to lower some unit test precision during Debian testing of statsmodels.) Josef
2. so far they were identical across many numpy releases
3. they are identical on other architectures (e.g. amd64)
To me that is surprising. I would have guessed that the order is the same on 32 and 64 bit, but something about the precision of intermediate operations is different. I don't know enough about amd64 to guess what that could be. Bradley's suggestion seems kind of reasonable but it's strange that numpy should use intel-80 bit intermediate values differently for 32 and 64 bit.
Cheers,
Matthew _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On Tue, Apr 30, 2013 at 8:08 PM, Yaroslav Halchenko <lists@onerussian.com>wrote:
could anyone on 32bit system with fresh numpy (1.7.1) test following:
wget -nc http://www.onerussian.com/tmp/data.npy ; python -c 'import numpy as np; data1 = np.load("/tmp/data.npy"); print np.sum(data1[1,:,0,1]) - np.sum(data1, axis=1)[1,0,1]'
0.0
because unfortunately it seems on fresh ubuntu raring (in 32bit build only, seems ok in 64 bit... also never ran into it on older numpy releases):
python -c 'import numpy as np; data1 = np.load("/tmp/data.npy"); print np.sum(data1[1,:,0,1]) - np.sum(data1, axis=1)[1,0,1]' -1.11022302463e-16
Perhaps on the 32-bit system one call is using the 80-bit extended precision register for the summation and the other one is not? -Brad
On Apr 30, 2013, at 6:37 PM, Benjamin Root <ben.root@ou.edu> wrote:
I can not think of any reason not to include these functions in v1.8.
+1
Of course, the documentation for discussed before: np.minmax(). My thinking is that it would return a 2xN array
How about a tuple: (min, max)? -Chris
On Wed, May 1, 2013 at 1:13 AM, Chris Barker - NOAA Federal < chris.barker@noaa.gov> wrote:
Of course, the documentation for discussed before: np.minmax(). My thinking is that it would return a 2xN array
How about a tuple: (min, max)?
I am not familiar enough with numpy internals to know which is the better approach to implement. I kind of feel that the 2xN array approach would be more flexible in case a user wants all of this information in a single array, while still allowing for unpacking as if it was a tuple. I would rather enable unforeseen use-cases rather than needlessly restricting them. Ben Root
On Wed, May 1, 2013 at 6:52 AM, Benjamin Root <ben.root@ou.edu> wrote:
How about a tuple: (min, max)?
I am not familiar enough with numpy internals to know which is the better approach to implement. I kind of feel that the 2xN array approach would be more flexible in case a user wants all of this information in a single array, while still allowing for unpacking as if it was a tuple.
hmm, my thinking is that the min and max values really are two different results, so getting two separate arrays makes sense to me. however, you are right, Python's nifty generic sequence unpacking lets you use a (2X...) array similarly to a tuple of two arrays, so why not? Food for thought on one reason: min, max = np.minmax(arr) would result in two arrays, but they would be views on the same array, Not sure if that matters, though. -Chris -Chris -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/OR&R (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker@noaa.gov
So, to summarize the thread so far: Consensus: np.nanmean() np.nanstd() np.minmax() np.argminmax() Vague Consensus: np.sincos() No Consensus (possibly out of scope for this topic): Better constructors for complex types I can probably whip up the PR for the nanmean() and nanstd(), and can certainly help out with the minmax funcs. Cheers! Ben Root
On Wed, May 1, 2013 at 7:10 PM, Benjamin Root <ben.root@ou.edu> wrote:
So, to summarize the thread so far:
Consensus: np.nanmean() np.nanstd() np.minmax() np.argminmax()
Vague Consensus: np.sincos()
If the return of sincos (cossin?) is an array, then it could be reshaped to be exp(1j*x), which together with exp(2*pi*1j*x) would cover some pretty common cases. No Consensus (possibly out of scope for this topic):
Better constructors for complex types
I can probably whip up the PR for the nanmean() and nanstd(), and can certainly help out with the minmax funcs.
Cheers! Ben Root
Chuck
On 1 May 2013 23:12, "Charles R Harris" <charlesr.harris@gmail.com> wrote:
On Wed, May 1, 2013 at 7:10 PM, Benjamin Root <ben.root@ou.edu> wrote:
So, to summarize the thread so far:
Consensus: np.nanmean() np.nanstd() np.minmax() np.argminmax()
Vague Consensus: np.sincos()
If the return of sincos (cossin?) is an array, then it could be reshaped
to be exp(1j*x), which together with exp(2*pi*1j*x) would cover some pretty common cases. Ufuncs already have some convention for what to do with multiple output arguments, right? Presumably whatever they do is what sincos should do. (And minmax/argminmax likewise, for consistency, even if they aren't ufuncs. Though they could be generalized ufuncs, or minmax could be minimummaximum.reduce.) I haven't checked, but I assume that what multiple output argument ufuncs do is to return a tuple. You can't use a single array in the general case, because the multiple output types might not be homogenous. -n
On Thu, 2013-05-02 at 07:03 -0400, Nathaniel Smith wrote:
On 1 May 2013 23:12, "Charles R Harris" <charlesr.harris@gmail.com> wrote:
On Wed, May 1, 2013 at 7:10 PM, Benjamin Root <ben.root@ou.edu>
wrote:
So, to summarize the thread so far:
Consensus: np.nanmean() np.nanstd() np.minmax() np.argminmax()
Vague Consensus: np.sincos()
If the return of sincos (cossin?) is an array, then it could be reshaped to be exp(1j*x), which together with exp(2*pi*1j*x) would cover some pretty common cases.
Ufuncs already have some convention for what to do with multiple output arguments, right? Presumably whatever they do is what sincos should do. (And minmax/argminmax likewise, for consistency, even if they aren't ufuncs. Though they could be generalized ufuncs, or minmax could be minimummaximum.reduce.)
I think for sincos it makes sense, for an example a ufunc with multiple arguments is `np.modf`. But I doubt reductions are automatically defined for these, so minmax probably needs to be a generalized ufunc (can you have an axis argument with those?). - Sebastian
I haven't checked, but I assume that what multiple output argument ufuncs do is to return a tuple. You can't use a single array in the general case, because the multiple output types might not be homogenous.
-n
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On Thu, May 2, 2013 at 12:03 PM, Nathaniel Smith <njs@pobox.com> wrote:
On 1 May 2013 23:12, "Charles R Harris" <charlesr.harris@gmail.com> wrote:
On Wed, May 1, 2013 at 7:10 PM, Benjamin Root <ben.root@ou.edu> wrote:
So, to summarize the thread so far:
Consensus: np.nanmean() np.nanstd() np.minmax() np.argminmax()
Vague Consensus: np.sincos()
If the return of sincos (cossin?) is an array, then it could be reshaped to be exp(1j*x), which together with exp(2*pi*1j*x) would cover some pretty common cases.
It couldn't be a mere reshape, since the complex dtype requires the real and imag components to be adjacent to each other. They wouldn't be so if sincos's return type is an array (nor even the cossin alternative). It always requires a memory copy (except in the "who cares?" case of a scalar). Composition with an efficient np.tocomplex(real, imag) implementation would cover those use cases whether sincos returns tuples or arrays.
Ufuncs already have some convention for what to do with multiple output arguments, right? Presumably whatever they do is what sincos should do. (And minmax/argminmax likewise, for consistency, even if they aren't ufuncs. Though they could be generalized ufuncs, or minmax could be minimummaximum.reduce.)
I haven't checked, but I assume that what multiple output argument ufuncs do is to return a tuple. You can't use a single array in the general case, because the multiple output types might not be homogenous.
Correct. [~] |19> np.modf.nout 2 [~] |20> np.modf(np.linspace(0, 1, 5)) (array([ 0. , 0.25, 0.5 , 0.75, 0. ]), array([ 0., 0., 0., 0., 1.])) -- Robert Kern
On Thu, May 2, 2013 at 7:28 AM, Robert Kern <robert.kern@gmail.com> wrote:
On 1 May 2013 23:12, "Charles R Harris" <charlesr.harris@gmail.com> wrote:
On Wed, May 1, 2013 at 7:10 PM, Benjamin Root <ben.root@ou.edu> wrote:
So, to summarize the thread so far:
Consensus: np.nanmean() np.nanstd() np.minmax() np.argminmax()
Vague Consensus: np.sincos()
If the return of sincos (cossin?) is an array, then it could be reshaped to be exp(1j*x), which together with exp(2*pi*1j*x) would cover some
On Thu, May 2, 2013 at 12:03 PM, Nathaniel Smith <njs@pobox.com> wrote: pretty
common cases.
It couldn't be a mere reshape, since the complex dtype requires the real and imag components to be adjacent to each other. They wouldn't be so if sincos's return type is an array (nor even the cossin alternative). It always requires a memory copy (except in the "who cares?" case of a scalar). Composition with an efficient np.tocomplex(real, imag) implementation would cover those use cases whether sincos returns tuples or arrays.
I would assume the basic return type would be complex, i.e., the cos/sin adjacent. The cos/sin parts would then be real/imag views into the array.
Ufuncs already have some convention for what to do with multiple output arguments, right? Presumably whatever they do is what sincos should do. (And minmax/argminmax likewise, for consistency, even if they aren't ufuncs. Though they could be generalized ufuncs, or minmax could be minimummaximum.reduce.)
I haven't checked, but I assume that what multiple output argument ufuncs do is to return a tuple. You can't use a single array in the general case, because the multiple output types might not be homogenous.
Correct.
[~] |19> np.modf.nout 2
[~] |20> np.modf(np.linspace(0, 1, 5)) (array([ 0. , 0.25, 0.5 , 0.75, 0. ]), array([ 0., 0., 0., 0., 1.]))
-- Robert Kern _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On Thu, May 2, 2013 at 2:38 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Thu, May 2, 2013 at 7:28 AM, Robert Kern <robert.kern@gmail.com> wrote:
On Thu, May 2, 2013 at 12:03 PM, Nathaniel Smith <njs@pobox.com> wrote:
On 1 May 2013 23:12, "Charles R Harris" <charlesr.harris@gmail.com> wrote:
On Wed, May 1, 2013 at 7:10 PM, Benjamin Root <ben.root@ou.edu> wrote:
So, to summarize the thread so far:
Consensus: np.nanmean() np.nanstd() np.minmax() np.argminmax()
Vague Consensus: np.sincos()
If the return of sincos (cossin?) is an array, then it could be reshaped to be exp(1j*x), which together with exp(2*pi*1j*x) would cover some pretty common cases.
It couldn't be a mere reshape, since the complex dtype requires the real and imag components to be adjacent to each other. They wouldn't be so if sincos's return type is an array (nor even the cossin alternative). It always requires a memory copy (except in the "who cares?" case of a scalar). Composition with an efficient np.tocomplex(real, imag) implementation would cover those use cases whether sincos returns tuples or arrays.
I would assume the basic return type would be complex, i.e., the cos/sin adjacent. The cos/sin parts would then be real/imag views into the array.
You mean that the implementation of cossin (to make things easier on ourselves) would create an (N,2) contiguous array, fill it with the cos and sin results, then reshape it to return the expected (2,N) array (or 2-tuple)? How would the user then reconstitute the exp(1j*x) result efficiently? If the use case is that important, I would just make exp(1j*x) into its own ufunc and have it use the C sincos() function internally. -- Robert Kern
On Thu, May 2, 2013 at 7:47 AM, Robert Kern <robert.kern@gmail.com> wrote:
On Thu, May 2, 2013 at 2:38 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Thu, May 2, 2013 at 7:28 AM, Robert Kern <robert.kern@gmail.com>
wrote:
On Thu, May 2, 2013 at 12:03 PM, Nathaniel Smith <njs@pobox.com> wrote:
On 1 May 2013 23:12, "Charles R Harris" <charlesr.harris@gmail.com> wrote:
On Wed, May 1, 2013 at 7:10 PM, Benjamin Root <ben.root@ou.edu>
wrote:
So, to summarize the thread so far:
Consensus: np.nanmean() np.nanstd() np.minmax() np.argminmax()
Vague Consensus: np.sincos()
If the return of sincos (cossin?) is an array, then it could be reshaped to be exp(1j*x), which together with exp(2*pi*1j*x) would cover some pretty common cases.
It couldn't be a mere reshape, since the complex dtype requires the real and imag components to be adjacent to each other. They wouldn't be so if sincos's return type is an array (nor even the cossin alternative). It always requires a memory copy (except in the "who cares?" case of a scalar). Composition with an efficient np.tocomplex(real, imag) implementation would cover those use cases whether sincos returns tuples or arrays.
I would assume the basic return type would be complex, i.e., the cos/sin adjacent. The cos/sin parts would then be real/imag views into the array.
You mean that the implementation of cossin (to make things easier on ourselves) would create an (N,2) contiguous array, fill it with the cos and sin results, then reshape it to return the expected (2,N)
Just return the transpose. Chuck
On Thu, May 2, 2013 at 3:28 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Thu, May 2, 2013 at 7:47 AM, Robert Kern <robert.kern@gmail.com> wrote:
On Thu, May 2, 2013 at 2:38 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Thu, May 2, 2013 at 7:28 AM, Robert Kern <robert.kern@gmail.com> wrote:
On Thu, May 2, 2013 at 12:03 PM, Nathaniel Smith <njs@pobox.com> wrote:
On 1 May 2013 23:12, "Charles R Harris" <charlesr.harris@gmail.com> wrote:
On Wed, May 1, 2013 at 7:10 PM, Benjamin Root <ben.root@ou.edu> wrote: > > So, to summarize the thread so far: > > Consensus: > np.nanmean() > np.nanstd() > np.minmax() > np.argminmax() > > Vague Consensus: > np.sincos() >
If the return of sincos (cossin?) is an array, then it could be reshaped to be exp(1j*x), which together with exp(2*pi*1j*x) would cover some pretty common cases.
It couldn't be a mere reshape, since the complex dtype requires the real and imag components to be adjacent to each other. They wouldn't be so if sincos's return type is an array (nor even the cossin alternative). It always requires a memory copy (except in the "who cares?" case of a scalar). Composition with an efficient np.tocomplex(real, imag) implementation would cover those use cases whether sincos returns tuples or arrays.
I would assume the basic return type would be complex, i.e., the cos/sin adjacent. The cos/sin parts would then be real/imag views into the array.
You mean that the implementation of cossin (to make things easier on ourselves) would create an (N,2) contiguous array, fill it with the cos and sin results, then reshape it to return the expected (2,N)
Just return the transpose.
Yes, that's what I was getting at with that sentence. I don't doubt that that is possible. The problem I was pointing out was in the following sentence, which you snipped: "How would the user then reconstitute the exp(1j*x) result efficiently?" Please show me the code that the user would write to compute exp(1j*x) using np.cossin() without memory copies. My suspicion is that it will be non-intuitive enough that it should always be hidden away into a well-commented utility function. In that case, I think we should just provide an np.exp1j() ufunc that just uses the C sincos() function internally and let np.sincos()/np.cossin() do whatever is most natural and consistent with other nout>1 ufuncs freed from the constraints of this use case. -- Robert Kern
On Thu, May 2, 2013 at 8:40 AM, Robert Kern <robert.kern@gmail.com> wrote:
On Thu, May 2, 2013 at 3:28 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Thu, May 2, 2013 at 7:47 AM, Robert Kern <robert.kern@gmail.com>
wrote:
On Thu, May 2, 2013 at 2:38 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Thu, May 2, 2013 at 7:28 AM, Robert Kern <robert.kern@gmail.com> wrote:
On Thu, May 2, 2013 at 12:03 PM, Nathaniel Smith <njs@pobox.com>
wrote:
On 1 May 2013 23:12, "Charles R Harris" <charlesr.harris@gmail.com
wrote: > > On Wed, May 1, 2013 at 7:10 PM, Benjamin Root <ben.root@ou.edu> > wrote: >> >> So, to summarize the thread so far: >> >> Consensus: >> np.nanmean() >> np.nanstd() >> np.minmax() >> np.argminmax() >> >> Vague Consensus: >> np.sincos() >> > > If the return of sincos (cossin?) is an array, then it could be > reshaped > to be exp(1j*x), which together with exp(2*pi*1j*x) would cover some > pretty > common cases.
It couldn't be a mere reshape, since the complex dtype requires the real and imag components to be adjacent to each other. They wouldn't be so if sincos's return type is an array (nor even the cossin alternative). It always requires a memory copy (except in the "who cares?" case of a scalar). Composition with an efficient np.tocomplex(real, imag) implementation would cover those use cases whether sincos returns tuples or arrays.
I would assume the basic return type would be complex, i.e., the cos/sin adjacent. The cos/sin parts would then be real/imag views into the array.
You mean that the implementation of cossin (to make things easier on ourselves) would create an (N,2) contiguous array, fill it with the cos and sin results, then reshape it to return the expected (2,N)
Just return the transpose.
Yes, that's what I was getting at with that sentence. I don't doubt that that is possible. The problem I was pointing out was in the following sentence, which you snipped:
"How would the user then reconstitute the exp(1j*x) result efficiently?"
Please show me the code that the user would write to compute exp(1j*x) using np.cossin() without memory copies. My suspicion is that it will be non-intuitive enough that it should always be hidden away into a well-commented utility function. In that case, I think we should just provide an np.exp1j() ufunc that just uses the C sincos() function internally and let np.sincos()/np.cossin() do whatever is most natural and consistent with other nout>1 ufuncs freed from the constraints of this use case.
As you say, have two functions, one of which would use the other with a view/transpose, whatever. For instance, given exp1j(), have another function that returns the real/imag parts. The question is what the underlying function should be and for that I think exp1j() would be a good choice. Chuck
On Thu, May 2, 2013 at 3:57 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Thu, May 2, 2013 at 8:40 AM, Robert Kern <robert.kern@gmail.com> wrote:
On Thu, May 2, 2013 at 3:28 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Thu, May 2, 2013 at 7:47 AM, Robert Kern <robert.kern@gmail.com> wrote:
On Thu, May 2, 2013 at 2:38 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Thu, May 2, 2013 at 7:28 AM, Robert Kern <robert.kern@gmail.com> wrote:
On Thu, May 2, 2013 at 12:03 PM, Nathaniel Smith <njs@pobox.com> wrote: > On 1 May 2013 23:12, "Charles R Harris" > <charlesr.harris@gmail.com> > wrote: >> >> On Wed, May 1, 2013 at 7:10 PM, Benjamin Root <ben.root@ou.edu> >> wrote: >>> >>> So, to summarize the thread so far: >>> >>> Consensus: >>> np.nanmean() >>> np.nanstd() >>> np.minmax() >>> np.argminmax() >>> >>> Vague Consensus: >>> np.sincos() >>> >> >> If the return of sincos (cossin?) is an array, then it could be >> reshaped >> to be exp(1j*x), which together with exp(2*pi*1j*x) would cover >> some >> pretty >> common cases.
It couldn't be a mere reshape, since the complex dtype requires the real and imag components to be adjacent to each other. They wouldn't be so if sincos's return type is an array (nor even the cossin alternative). It always requires a memory copy (except in the "who cares?" case of a scalar). Composition with an efficient np.tocomplex(real, imag) implementation would cover those use cases whether sincos returns tuples or arrays.
I would assume the basic return type would be complex, i.e., the cos/sin adjacent. The cos/sin parts would then be real/imag views into the array.
You mean that the implementation of cossin (to make things easier on ourselves) would create an (N,2) contiguous array, fill it with the cos and sin results, then reshape it to return the expected (2,N)
Just return the transpose.
Yes, that's what I was getting at with that sentence. I don't doubt that that is possible. The problem I was pointing out was in the following sentence, which you snipped:
"How would the user then reconstitute the exp(1j*x) result efficiently?"
Please show me the code that the user would write to compute exp(1j*x) using np.cossin() without memory copies. My suspicion is that it will be non-intuitive enough that it should always be hidden away into a well-commented utility function. In that case, I think we should just provide an np.exp1j() ufunc that just uses the C sincos() function internally and let np.sincos()/np.cossin() do whatever is most natural and consistent with other nout>1 ufuncs freed from the constraints of this use case.
As you say, have two functions, one of which would use the other with a view/transpose, whatever. For instance, given exp1j(), have another function that returns the real/imag parts. The question is what the underlying function should be and for that I think exp1j() would be a good choice.
I don't see why we would bother. Just implement them both as ufuncs that use a C sincos() function internally and be done with it. Implementing one in terms of the other requires that the other is not a ufunc (a minor irritation) and always returns non-contiguous arrays (a more substantial irritation). There isn't anything to be gained by implementing one in terms of the other. -- Robert Kern
I have created a PR for the first two (and got np.nanvar() for free). https://github.com/numpy/numpy/pull/3297 Cheers! Ben Root
On 1 May 2013 03:36, Benjamin Root <ben.root@ou.edu> wrote:
Are there any other functions that others feel are "missing" from numpy and would like to see for v1.8? Let's discuss them here.
I would like to have sincos, to compute sin and cos of the same number faster. According to some benchmarks, it is barely slower than just computing one of them. On 1 May 2013 07:13, Chris Barker - NOAA Federal <chris.barker@noaa.gov> wrote:
Of course, the documentation for discussed before: np.minmax(). My thinking is that it would return a 2xN array
How about a tuple: (min, max)?
Consider the case of np.minmax(matrix, axis=1), you will end up with a tuple of two arrays. In that scenario, you probably want to do computations with both numbers, so having them in an array seems more convenient. If there is enough reason, we could always add a "unpack=True" flag and then return a tuple.
On 05/01/2013 04:14 PM, Daπid wrote:
On 1 May 2013 03:36, Benjamin Root <ben.root@ou.edu> wrote:
Are there any other functions that others feel are "missing" from numpy and would like to see for v1.8? Let's discuss them here. I would like to have sincos, to compute sin and cos of the same number faster. According to some benchmarks, it is barely slower than just computing one of them.
Where does this `sincos` function come from?
The sincos function is in the c standard library in math.h. On May 1, 2013 7:56 AM, "Juan Luis Cano" <juanlu001@gmail.com> wrote:
On 05/01/2013 04:14 PM, Daπid wrote:
On 1 May 2013 03:36, Benjamin Root <ben.root@ou.edu> wrote:
Are there any other functions that others feel are "missing" from numpy and would like to see for v1.8? Let's discuss them here. I would like to have sincos, to compute sin and cos of the same number faster. According to some benchmarks, it is barely slower than just computing one of them.
Where does this `sincos` function come from? _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On Wed, May 1, 2013 at 4:06 PM, Zachary Ploskey <zploskey@gmail.com> wrote:
The sincos function is in the c standard library in math.h.
I don't think it's part of the C99 standard. It appears to be provided in glibc as a non-standard extension. We would have to provide our own copy, but one is available in the Cephes library. -- Robert Kern
On Wed, May 1, 2013 at 10:14 AM, Daπid <davidmenhur@gmail.com> wrote:
On 1 May 2013 03:36, Benjamin Root <ben.root@ou.edu> wrote:
Are there any other functions that others feel are "missing" from numpy and would like to see for v1.8? Let's discuss them here.
I would like to have sincos, to compute sin and cos of the same number faster. According to some benchmarks, it is barely slower than just computing one of them.
+1 Warren
Of course, the documentation for discussed before: np.minmax(). My
On 1 May 2013 07:13, Chris Barker - NOAA Federal <chris.barker@noaa.gov> wrote: thinking is that it would return a 2xN array
How about a tuple: (min, max)?
Consider the case of np.minmax(matrix, axis=1), you will end up with a tuple of two arrays. In that scenario, you probably want to do computations with both numbers, so having them in an array seems more convenient.
If there is enough reason, we could always add a "unpack=True" flag and then return a tuple. _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On Wed, May 1, 2013 at 3:36 AM, Benjamin Root <ben.root@ou.edu> wrote:
Are there any other functions that others feel are "missing" from numpy and would like to see for v1.8? Let's discuss them here.
As I mentioned before, I think numpy should have some equations for dealing with n-dimensional vectors (but would also work with complex dtypes). This would include n-dimensional equivalents of np.abs and np.angle, as well as a function to go back to an n-dimensional vector from the length and angle. Considering how critical vector analysis is to signal processing, I am surprised these don't already exist. There aren't even functions that work with 2-dimensional vectors, you have to construct a complex array first (which isn't that easy to do). Speaking of which, I think there should be a function to construct a complex array out of two identically-shaped floating-point arrays, as well as perhaps an np.i class that converts a real array to an imaginary one (using __mul__ and such).
On 1 May 2013 17:13, Todd <toddrjen@gmail.com> wrote:
Speaking of which, I think there should be a function to construct a complex array out of two identically-shaped floating-point arrays, as well as perhaps an np.i class that converts a real array to an imaginary one (using __mul__ and such).
np.i would be exactly the same as array * 1j, or am I missing anything? The same goes for constructing a complex, real + imag * 1j
On Wed, May 1, 2013 at 5:22 PM, Daπid <davidmenhur@gmail.com> wrote:
On 1 May 2013 17:13, Todd <toddrjen@gmail.com> wrote:
Speaking of which, I think there should be a function to construct a complex array out of two identically-shaped floating-point arrays, as well as perhaps an np.i class that converts a real array to an imaginary one (using __mul__ and such).
np.i would be exactly the same as array * 1j, or am I missing anything?
The same goes for constructing a complex, real + imag * 1j
it would always produce a numpy array. So array*1j and array*np.i (or np.j if you prefer, perhaps both) would be the same, but list*1j and list*np.i would not. The function version would also probably allow you to specify the dtype, which 1j does not.
On Wed, May 1, 2013 at 4:22 PM, Daπid <davidmenhur@gmail.com> wrote:
On 1 May 2013 17:13, Todd <toddrjen@gmail.com> wrote:
Speaking of which, I think there should be a function to construct a complex array out of two identically-shaped floating-point arrays, as well as perhaps an np.i class that converts a real array to an imaginary one (using __mul__ and such).
np.i would be exactly the same as array * 1j, or am I missing anything?
I don't think we have a ufunc loop for multiply() that takes a float64 and a complex128 and returns a complex128. We just have a (complex128,complex128)->complex128. `x * 1j` first converts `x` to a complex128 array with the value in the real component, then multiplies that with 1j to move that value over to the imag component. A single operation that takes a float64 array and just makes a complex128 array with the values put in the imag component will reduce a temporary.
The same goes for constructing a complex, real + imag * 1j
Similarly, we can eliminate two temporaries here. Both of the cases are probably best addressed by a single function. The syntactic sugar of an np.i object is unnecessary, IMO. imag = np.tocomplex(0.0, x) z = np.tocomplex(x, y) -- Robert Kern
participants (16)
-
Benjamin Root -
Bradley M. Froehle -
Charles R Harris -
Chris Barker - NOAA Federal -
Daπid -
josef.pktd＠gmail.com -
Juan Luis Cano -
Matthew Brett -
Nathaniel Smith -
Pauli Virtanen -
Robert Kern -
Sebastian Berg -
Todd -
Warren Weckesser -
Yaroslav Halchenko -
Zachary Ploskey