lundi 6 juin 2011

Broadcasting: adding arrays of different shapes

Something very different between python-numpy and IDL is the way they are both dealing with linear algebra in case of operations on different size or shape arrays.

In IDL, trying to operate on 2 vectors with different size reduces the oprtation to the lowest size:
IDL> a=[1,2,3]
IDL> b=[10,20,30,40]
IDL> print,a+b
       11       22       33

As we can see, the latest element is purely omitted.
IDL> c=[[1,2,3,4]]
IDL> d=[[10],[20],[30],[40]]
IDL> print,c
       1       2       3       4
IDL> print,d
      10
      20
      30
      40

The shape of the first operand is conserved:
IDL> print,c+d
      11      22      33      44
IDL> print,d+c
      11
      22
      33
      44

And we can loose part of the vector:
IDL> print,a+c
       2       4       6
IDL> print,a+d
      11      22      33



On the contrary, Python-numpy is adding some information, this is the broadcasting.

import numpy as np
a=np.array([1,2,3.])
b=np.array([10,20,30.,40])

a*b
# Traceback (most recent call last):
#  File "<stdin>", line 1, in <module>
#ValueError: shape mismatch: objects cannot be broadcast to a single shape



BUT:

a = a.reshape(3,1)
a
#array([[ 1.],

#      [ 2.],
#       [ 3.]])

b
#array([ 10.,  20.,  30.,  40.])


a*b
#array([[  10.,   20.,   30.,   40.],
#       [  20.,   40.,   60.,   80.],
#       [  30.,   60.,   90.,  120.]])


Another example where Python is guessing in which dimension the array must be extended:
c= a*b
c.shape
# (3, 4)
c*a
#array([[  10.,   20.,   30.,   40.],
#       [  40.,   80.,  120.,  160.],
#       [  90.,  180.,  270.,  360.]])
c*b
#array([[  100.,   400.,   900.,  1600.],
#       [  200.,   800.,  1800.,  3200.],
#       [  300.,  1200.,  2700.,  4800.]])
b2==np.array([10,20,30.])
#array([False, False, False], dtype=bool)
d=a*b2
d.shape
#(3, 3)
d
#array([[ 1.,  2.,  3.],
#       [ 2.,  4.,  6.],
3       [ 3.,  6.,  9.]])
d*a
#array([[  1.,   2.,   3.],

#       [  4.,   8.,  12.],
#       [  9.,  18.,  27.]])
d*b2
#array([[  1.,   4.,   9.],
#       [  2.,   8.,  18.],
#       [  3.,  12.,  27.]])


More details and nice figures in:

http://www.scipy.org/EricsBroadcastingDoc

1 commentaire: