broacasting question
I'm trying to convert some IDL code to python/numpy and i'm having some trouble understanding the rules for boradcasting during some operations. example: given the following arrays: a = array((2048,3577), dtype=float) b = array((256,25088), dtype=float) c = array((2048,3136), dtype=float) d = array((2048,3136), dtype=float) do: a = b * c + d In IDL, the computation is done without complaint and all array sizes are preserved. In ptyhon I get a value error concerning broadcasting. I can force it to work by taking slices, but the resulting size would be a = (256x3136) rather than (2048x3577). I admit that I don't understand IDL (or python to be honest) well enough to know how it handles this to be able to replicate the result properly. Does it only operate on the smallest dimensions ignoring the larger indices leaving their values unchanged? Can someone explain this to me? -- Thomas K. Gamble tkgamble@windstream.net
2011/6/30 Thomas K Gamble <tkgamble@windstream.net>
I'm trying to convert some IDL code to python/numpy and i'm having some trouble understanding the rules for boradcasting during some operations. example:
given the following arrays: a = array((2048,3577), dtype=float) b = array((256,25088), dtype=float) c = array((2048,3136), dtype=float) d = array((2048,3136), dtype=float)
do: a = b * c + d
In IDL, the computation is done without complaint and all array sizes are preserved. In ptyhon I get a value error concerning broadcasting. I can force it to work by taking slices, but the resulting size would be a = (256x3136) rather than (2048x3577). I admit that I don't understand IDL (or python to be honest) well enough to know how it handles this to be able to replicate the result properly. Does it only operate on the smallest dimensions ignoring the larger indices leaving their values unchanged? Can someone explain this to me?
-- Thomas K. Gamble tkgamble@windstream.net _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
It should be all explained here: http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html -=- Olivier
On 30.06.2011, at 7:32PM, Thomas K Gamble wrote:
I'm trying to convert some IDL code to python/numpy and i'm having some trouble understanding the rules for boradcasting during some operations. example:
given the following arrays: a = array((2048,3577), dtype=float) b = array((256,25088), dtype=float) c = array((2048,3136), dtype=float) d = array((2048,3136), dtype=float)
do: a = b * c + d
In IDL, the computation is done without complaint and all array sizes are preserved. In ptyhon I get a value error concerning broadcasting. I can force it to work by taking slices, but the resulting size would be a = (256x3136) rather than (2048x3577). I admit that I don't understand IDL (or python to be honest) well enough to know how it handles this to be able to replicate the result properly. Does it only operate on the smallest dimensions ignoring the larger indices leaving their values unchanged? Can someone explain this to me?
If IDL does such operations silently I'd probably rather be troubled about it... Assuming you actually meant something like "a = np.ndarray((2048,3577))" (because np.array((2048,3577), dtype=float) would simply be the 2-vector [ 2048. 3577.]), the shape of a indeed matches in no way the other ones. While b,c,d do have the same total size, thus something like b.reshape((2048,3136) * c + d will work, meaning the first 8 rows of b b[:8] would be concatenated to the first row of the output, and so on... Since the total size is still smaller than a, I could only venture something is done like np.add(b.reshape(2048,3136) * c, d, out=a[:,:3136]) But to say whether this is really the equivalent result to what IDL does, one would have to study the IDL manual in detail or directly compare the output (e.g. check what happens to the values in a[:,3136:]...) Cheers, Derek
On 30.06.2011, at 7:32PM, Thomas K Gamble wrote:
I'm trying to convert some IDL code to python/numpy and i'm having some trouble understanding the rules for boradcasting during some operations. example:
given the following arrays: a = array((2048,3577), dtype=float) b = array((256,25088), dtype=float) c = array((2048,3136), dtype=float) d = array((2048,3136), dtype=float)
do: a = b * c + d
In IDL, the computation is done without complaint and all array sizes are preserved. In ptyhon I get a value error concerning broadcasting. I can force it to work by taking slices, but the resulting size would be a = (256x3136) rather than (2048x3577). I admit that I don't understand IDL (or python to be honest) well enough to know how it handles this to be able to replicate the result properly. Does it only operate on the smallest dimensions ignoring the larger indices leaving their values unchanged? Can someone explain this to me?
If IDL does such operations silently I'd probably rather be troubled about it... Assuming you actually meant something like "a = np.ndarray((2048,3577))" (because np.array((2048,3577), dtype=float) would simply be the 2-vector [ 2048. 3577.]), the shape of a indeed matches in no way the other ones. While b,c,d do have the same total size, thus something like
b.reshape((2048,3136) * c + d
will work, meaning the first 8 rows of b b[:8] would be concatenated to the first row of the output, and so on... Since the total size is still smaller than a, I could only venture something is done like
np.add(b.reshape(2048,3136) * c, d, out=a[:,:3136])
But to say whether this is really the equivalent result to what IDL does, one would have to study the IDL manual in detail or directly compare the output (e.g. check what happens to the values in a[:,3136:]...)
Cheers, Derek
Your post gave me the cluse I needed. I had my shapes slightly off in the example I gave, but if I try: a = reshape(b.flatten('F') * c.flatten('F') + d.flatten('F'), b.shape, order='F') I get a result in line with the IDL result. Another example with different total size arrays: b = np.ndarray((2048,3577), dtype=float) c = np.ndarray((256,25088), dtype=float) a= reshape(b.flatten('F')[:size(c)]/c.flatten('F'), c.shape, order='F') This also gives a result like that of IDL. Thanks for the help. -- Thomas K. Gamble tkgamble@windstream.net Receive my instruction, and not silver; and knowledge rather than choice gold. For wisdom is better than rubies; and all the things that may be desired are not to be compared to it. (Proverbs 8:10,11)
On 30.06.2011, at 11:57PM, Thomas K Gamble wrote:
np.add(b.reshape(2048,3136) * c, d, out=a[:,:3136])
But to say whether this is really the equivalent result to what IDL does, one would have to study the IDL manual in detail or directly compare the output (e.g. check what happens to the values in a[:,3136:]...)
Cheers, Derek
Your post gave me the cluse I needed.
I had my shapes slightly off in the example I gave, but if I try:
a = reshape(b.flatten('F') * c.flatten('F') + d.flatten('F'), b.shape, order='F')
I get a result in line with the IDL result.
Another example with different total size arrays:
b = np.ndarray((2048,3577), dtype=float) c = np.ndarray((256,25088), dtype=float)
a= reshape(b.flatten('F')[:size(c)]/c.flatten('F'), c.shape, order='F')
This also gives a result like that of IDL.
Right, I forgot to point out that there are at least 2 ways to bring the arrays into compatible shapes (that's the reason broadcasting does not work here, because numpy only does automatic broadcasting if there is an unambiguous way to do so). So the IDL arrays being Fortran-ordered is the essential bit of information here. Just two remarks: I. Assigning a = reshape(b.flatten('F')[:size(c)]/c.flatten('F'), c.shape, order='F') as above will create a new array of shape c.shape - if you wanted to put your results into an existing array of shape(2048,3577), you'd still have to explicitly say a[:,:3136] = ... II. The flatten() operations and the assignment above all create full copies of the arrays, thus the np.add ufunc above together with simple reshape operations might improve performance somewhat - however keeping the Fortran order also requires some costly transpositions, as for your last example a = np.divide(b.T[:3136].reshape(c.T.shape).T, c, out=a) so YMMV... Cheers, Derek
Right, I forgot to point out that there are at least 2 ways to bring the arrays into compatible shapes (that's the reason broadcasting does not work here, because numpy only does automatic broadcasting if there is an unambiguous way to do so). So the IDL arrays being Fortran-ordered is the essential bit of information here. Just two remarks: I. Assigning a = reshape(b.flatten('F')[:size(c)]/c.flatten('F'), c.shape, order='F') as above will create a new array of shape c.shape - if you wanted to put your results into an existing array of shape(2048,3577), you'd still have to explicitly say a[:,:3136] = ... II. The flatten()
That was the error in my example I refered to.
operations and the assignment above all create full copies of the arrays, thus the np.add ufunc above together with simple reshape operations might improve performance somewhat - however keeping the Fortran order also requires some costly transpositions, as for your last example
a = np.divide(b.T[:3136].reshape(c.T.shape).T, c, out=a)
so YMMV...
Cheers, Derek
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion -- Thomas K. Gamble tkgamble@windstream.net The fruit of the righteous is a tree of life; and he who wins souls is wise. (Proverbs 11:30)
On 30.06.2011, at 11:57PM, Thomas K Gamble wrote:
np.add(b.reshape(2048,3136) * c, d, out=a[:,:3136])
But to say whether this is really the equivalent result to what IDL does, one would have to study the IDL manual in detail or directly compare the output (e.g. check what happens to the values in a[:,3136:]...)
Cheers,
Derek
Your post gave me the cluse I needed.
I had my shapes slightly off in the example I gave, but if I try:
a = reshape(b.flatten('F') * c.flatten('F') + d.flatten('F'), b.shape, order='F')
I get a result in line with the IDL result.
Another example with different total size arrays:
b = np.ndarray((2048,3577), dtype=float) c = np.ndarray((256,25088), dtype=float)
a= reshape(b.flatten('F')[:size(c)]/c.flatten('F'), c.shape, order='F')
This also gives a result like that of IDL.
Right, I forgot to point out that there are at least 2 ways to bring the arrays into compatible shapes (that's the reason broadcasting does not work here, because numpy only does automatic broadcasting if there is an unambiguous way to do so). So the IDL arrays being Fortran-ordered is the essential bit of information here. Just two remarks: I. Assigning a = reshape(b.flatten('F')[:size(c)]/c.flatten('F'), c.shape, order='F') as above will create a new array of shape c.shape - if you wanted to put your results into an existing array of shape(2048,3577), you'd still have to explicitly say a[:,:3136] = ... II. The flatten()
That was the error in my first example I refered to. I confused it with the second 'divide' example and probably should have used different variable names to avoid confusing things further. Sorry.
operations and the assignment above all create full copies of the arrays, thus the np.add ufunc above together with simple reshape operations might improve performance somewhat - however keeping the Fortran order also requires some costly transpositions, as for your last example
Right now, I'm just interested in getting the right answer. Once I have that, I'll work on performance. Unfortunately the order does seem to make a difference.
a = np.divide(b.T[:3136].reshape(c.T.shape).T, c, out=a)
Interesting.
so YMMV...
Cheers, Derek
Thanks for your help. -- Thomas K. Gamble tkgamble@windstream.net LANL employee waiting out the Las Conchas fire.
On Thu, Jun 30, 2011 at 11:32 AM, Thomas K Gamble <tkgamble@windstream.net>wrote:
I'm trying to convert some IDL code to python/numpy and i'm having some trouble understanding the rules for boradcasting during some operations. example:
given the following arrays: a = array((2048,3577), dtype=float) b = array((256,25088), dtype=float) c = array((2048,3136), dtype=float) d = array((2048,3136), dtype=float)
do: a = b * c + d
In IDL, the computation is done without complaint and all array sizes are preserved. In ptyhon I get a value error concerning broadcasting. I can force it to work by taking slices, but the resulting size would be a = (256x3136) rather than (2048x3577). I admit that I don't understand IDL (or python to be honest) well enough to know how it handles this to be able to replicate the result properly. Does it only operate on the smallest dimensions ignoring the larger indices leaving their values unchanged? Can someone explain this to me?
I don't see a problem In [1]: datetime64('now') Out[1]: numpy.datetime64('2011-07-01T07:18:35-0600') In [2]: a = array((2048, 3577), float) In [3]: b = array((256, 25088), float) In [4]: c = array((2048, 3136), float) In [5]: d = array((2048, 3136), float) In [6]: a = b*c + d In [7]: a Out[7]: array([ 526336., 78679104.]) What is the '*' in your expression supposed to mean? Chuck
On Thu, Jun 30, 2011 at 11:32 AM, Thomas K Gamble
<tkgamble@windstream.net>wrote:
I'm trying to convert some IDL code to python/numpy and i'm having some trouble understanding the rules for boradcasting during some operations. example:
given the following arrays: a = array((2048,3577), dtype=float) b = array((256,25088), dtype=float) c = array((2048,3136), dtype=float) d = array((2048,3136), dtype=float)
do: a = b * c + d
In IDL, the computation is done without complaint and all array sizes are preserved. In ptyhon I get a value error concerning broadcasting. I can force it to work by taking slices, but the resulting size would be a = (256x3136) rather than (2048x3577). I admit that I don't understand IDL (or python to be honest) well enough to know how it handles this to be able to replicate the result properly. Does it only operate on the smallest dimensions ignoring the larger indices leaving their values unchanged? Can someone explain this to me?
I don't see a problem
In [1]: datetime64('now') Out[1]: numpy.datetime64('2011-07-01T07:18:35-0600')
In [2]: a = array((2048, 3577), float)
In [3]: b = array((256, 25088), float)
In [4]: c = array((2048, 3136), float)
In [5]: d = array((2048, 3136), float)
In [6]: a = b*c + d
In [7]: a Out[7]: array([ 526336., 78679104.])
What is the '*' in your expression supposed to mean?
My apologies for the errors in my example. It should have been: a = numpy.ndarray((2048,3577), dtype=float) b = numpy.ndarray((256,25088), dtype=float) c = numpy.ndarray((2048,3136), dtype=float) d = numpy.ndarray((2048,3136), dtype=float) The numbers are the array dimensions. Data values are not provided in the example. e = b * c + d f = a / b Both of these expressions result in value errors in python but IDL handles them without complaint. The * is a multiplication operator. IDL also stores its data in Fortran/column-major order, which causes some other issues.
Chuck
-- Thomas K. Gamble tkgamble@windstream.net
On Fri, Jul 1, 2011 at 6:38 PM, Thomas K Gamble <tkgamble@windstream.net>wrote:
On Thu, Jun 30, 2011 at 11:32 AM, Thomas K Gamble
<tkgamble@windstream.net>wrote:
I'm trying to convert some IDL code to python/numpy and i'm having some trouble understanding the rules for boradcasting during some operations. example:
given the following arrays: a = array((2048,3577), dtype=float) b = array((256,25088), dtype=float) c = array((2048,3136), dtype=float) d = array((2048,3136), dtype=float)
do: a = b * c + d
In IDL, the computation is done without complaint and all array sizes are preserved. In ptyhon I get a value error concerning broadcasting. I can force it to work by taking slices, but the resulting size would be a = (256x3136) rather than (2048x3577). I admit that I don't understand IDL (or python to be honest) well enough to know how it handles this to be able to replicate the result properly. Does it only operate on the smallest dimensions ignoring the larger indices leaving their values unchanged? Can someone explain this to me?
I don't see a problem
In [1]: datetime64('now') Out[1]: numpy.datetime64('2011-07-01T07:18:35-0600')
In [2]: a = array((2048, 3577), float)
In [3]: b = array((256, 25088), float)
In [4]: c = array((2048, 3136), float)
In [5]: d = array((2048, 3136), float)
In [6]: a = b*c + d
In [7]: a Out[7]: array([ 526336., 78679104.])
What is the '*' in your expression supposed to mean?
My apologies for the errors in my example. It should have been:
a = numpy.ndarray((2048,3577), dtype=float) b = numpy.ndarray((256,25088), dtype=float) c = numpy.ndarray((2048,3136), dtype=float) d = numpy.ndarray((2048,3136), dtype=float)
The numbers are the array dimensions. Data values are not provided in the example.
Ah, what that does in make (2,) arrays with the given elements and '*' and '+' are element-wise multiplication and addition. To get arrays with the dimensions you need something like In [19]: a = numpy.zeros((2048,3577), dtype=float) In [20]: b = numpy.zeros((256,25088), dtype=float) In [21]: c = numpy.zeros((2048,3136), dtype=float) In [22]: d = numpy.zeros((2048,3136), dtype=float) However, broadcasting b and c won't work, it isn't enough that 256 divides 2048, 256 must actually equal 1, which you can omit since it is a leading index. Same with 3136 and 25088 except the one won't automagically be added. So you can do things like In [24]: numpy.zeros((5,7))*numpy.zeros((7,)) Out[24]: array([[ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.]]) second array will be broadcast over the first. Chuck
participants (4)
-
Charles R Harris -
Derek Homeier -
Olivier Delalleau -
Thomas K Gamble