-
Notifications
You must be signed in to change notification settings - Fork 4
/
03-broadcasting.html
280 lines (275 loc) · 13.8 KB
/
03-broadcasting.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8">
<meta name="generator" content="pandoc">
<title>Software Carpentry: Advanced NumPy</title>
<link rel="shortcut icon" type="image/x-icon" href="/favicon.ico" />
<meta name="viewport" content="width=device-width, initial-scale=1.0" />
<link rel="stylesheet" type="text/css" href="css/bootstrap/bootstrap.css" />
<link rel="stylesheet" type="text/css" href="css/bootstrap/bootstrap-theme.css" />
<link rel="stylesheet" type="text/css" href="css/swc.css" />
<link rel="alternate" type="application/rss+xml" title="Software Carpentry Blog" href="http://software-carpentry.org/feed.xml"/>
<meta charset="UTF-8" />
<!-- HTML5 shim, for IE6-8 support of HTML5 elements -->
<!--[if lt IE 9]>
<script src="http://html5shim.googlecode.com/svn/trunk/html5.js"></script>
<![endif]-->
</head>
<body class="lesson">
<div class="container card">
<div class="banner">
<a href="http://software-carpentry.org" title="Software Carpentry">
<img alt="Software Carpentry banner" src="img/software-carpentry-banner.png" />
</a>
</div>
<article>
<div class="row">
<div class="col-md-10 col-md-offset-1">
<a href="index.html"><h1 class="title">Advanced NumPy</h1></a>
<h2 class="subtitle">Broadcasting</h2>
<section class="objectives panel panel-warning">
<div class="panel-heading">
<h2 id="learning-objectives"><span class="glyphicon glyphicon-certificate"></span>Learning Objectives</h2>
</div>
<div class="panel-body">
<p>After the lesson learner:</p>
<ul>
<li>Knows how to add a scalar to all elements of an array.</li>
<li>Can predict the result of additions of matrices and row or column vectors.</li>
<li>Can explain why broadcasting is better than using for loops.</li>
<li>Understands the rules of broadcasting and can predict the shape of broadcasted arrays.</li>
<li>Knows how to control broadcasting using <code>np.newaxis</code> object.</li>
</ul>
</div>
</section>
<p>It’s possible to do operations on arrays of different sizes. In some case NumPy can transform these arrays automatically so that they all have the same size: this conversion is called <strong>broadcasting</strong>.</p>
<div class="figure">
<img src="fig/numpy_broadcasting.png" title="numpy broadcasting in 2D" alt="numpy broadcasting in 2D" />
<p class="caption">numpy broadcasting in 2D</p>
</div>
<p>Let’s try to reproduce the above diagram. First, we create two one-dimensional arrays:</p>
<pre><code>>>> a = np.arange(4) * 10
>>> b = np.arange(3)
>>> a
array([ 0, 10, 20, 30])
>>> b
array([0, 1, 2])</code></pre>
<p>We can tile them in 2D using <code>np.tile</code> function:</p>
<pre><code>>>> b2 = np.tile(b, (4, 1))
>>> b2
array([[0, 1, 2],
[0, 1, 2],
[0, 1, 2],
[0, 1, 2]])</code></pre>
<p>We do the same with the second array, but we need also to transpose (exchange columns with rows) the resulting array:</p>
<pre><code>>>> a2 = np.tile(a, (3, 1))
>>> a2 = a2.T
>>> a2
array([[ 0, 0, 0],
[10, 10, 10],
[20, 20, 20],
[30, 30, 30]])</code></pre>
<p>Note that the <code>np.tile</code> function creates new arrays and copies the data. Then you can add the arrays element-wise:</p>
<pre><code>>>> a2 + b2
array([[ 0, 1, 2],
[10, 11, 12],
[20, 21, 22],
[30, 31, 32]])</code></pre>
<p>In the second example we add a one-dimensional array to a two-dimensional array. NumPy will automatically “tile” the 1D array along the missing direction:</p>
<pre><code>>>> a2 + b
array([[ 0, 1, 2],
[10, 11, 12],
[20, 21, 22],
[30, 31, 32]])</code></pre>
<p>However, in this case no copy of <code>b</code> array is involved. NumPy will instead use the same data in <code>b</code> for each row of <code>a</code> – we will cover the mechanism behind it at the end of the lesson.</p>
<p>In the third example we add a single column with a single vector. To obtain a column array from a 1D array we need to convert it to 2D array of four rows and one column. In NumPy we can add singular dimensions (dimensions of size 1) by a special object <code>np.newaxis</code>:</p>
<pre><code>>>> a.shape
(4,)
>>> a_column = a[:, np.newaxis]
>>> a_column.shape
(4, 1)
>>> a_column
array([[ 0],
[10],
[20],
[30]])</code></pre>
<p>We can add a column vector and a 1D array:</p>
<pre><code>>>> a_column + b
array([[ 0, 1, 2],
[10, 11, 12],
[20, 21, 22],
[30, 31, 32]])</code></pre>
<p>This is the same as adding a column and row vector:</p>
<pre><code>>>> b_row = b[np.newaxis, :]
>>> b_row
array([[0, 1, 2]])
>>> b_row.shape
(1, 3)
>>> a_column + b_row
array([[ 0, 1, 2],
[10, 11, 12],
[20, 21, 22],
[30, 31, 32]])</code></pre>
<section class="challenge panel panel-success">
<div class="panel-heading">
<h2 id="normalising-data"><span class="glyphicon glyphicon-pencil"></span>Normalising data</h2>
</div>
<div class="panel-body">
<p>Given the following array:</p>
<pre><code>a = np.random.rand(10, 100) </code></pre>
<p>For each column of <code>a</code> subtract its mean. Next, do the same with rows.</p>
</div>
</section>
<p>Broadcasting seems a bit magical, but it is actually quite natural to use it when we want to solve a problem whose output data is an array with more dimensions than input data. There a simple rule that allow to determine the validity of broadcasting and the shape of broadcasted arrays:</p>
<blockquote>
<p>In order to broadcast, the size of the trailing axes for both arrays in an operation must either be the same or one of them must be one.</p>
</blockquote>
<p>This does indeed work for the three addition from the figure</p>
<pre><code>a: 4 x 3
b: 4 x 3
result: 4 x 3
a: 4 x 3
b: 3
result: 4 x 3
a: 4 x 1
b: 3
result: 4 x 3</code></pre>
<p>Lets look at two more examples:</p>
<pre><code>Image (3d array): 256 x 256 x 3
Scale (1d array): 3
Result (3d array): 256 x 256 x 3
A (4d array): 8 x 1 x 6 x 1
B (3d array): 7 x 1 x 5
Result (4d array): 8 x 7 x 6 x 5</code></pre>
<section class="challenge panel panel-success">
<div class="panel-heading">
<h2 id="broadcasting-rules"><span class="glyphicon glyphicon-pencil"></span>Broadcasting rules</h2>
</div>
<div class="panel-body">
<p>Given the arrays:</p>
<pre><code>X = np.random.rand(10,3)
Y = np.random.rand(3)</code></pre>
<p>which of the following will <em>not</em> produce an error:</p>
<ol style="list-style-type: lower-alpha">
<li><p><code>X + Y</code></p></li>
<li><p><code>X[np.newaxis, :] + Y</code></p></li>
<li><p><code>X + Y[:, np.newaxis]</code></p></li>
<li><p><code>X[:, np.newaxis] + Y</code></p></li>
<li><p><code>X + Y[np.newaxis, :]</code></p></li>
<li><p><code>X[:, np.newaxis, :] + Y</code></p></li>
</ol>
<p>What will be the shapes of the final broadcasted arrays? Try to guess and then check.</p>
</div>
</section>
<section class="challenge panel panel-success">
<div class="panel-heading">
<h2 id="three-dimensional-broadcasting"><span class="glyphicon glyphicon-pencil"></span>Three-dimensional broadcasting</h2>
</div>
<div class="panel-body">
<p>Below, produce the array containing the sum of every element in <code>x</code> with every element in <code>y</code></p>
<div class="sourceCode"><pre class="sourceCode python"><code class="sourceCode python">x <span class="op">=</span> np.random.rand(<span class="dv">3</span>, <span class="dv">5</span>)
y <span class="op">=</span> np.random.randint(<span class="dv">10</span>, size<span class="op">=</span><span class="dv">8</span>)
z <span class="op">=</span> x <span class="co"># FIX THIS</span></code></pre></div>
</div>
</section>
<section class="challenge panel panel-success">
<div class="panel-heading">
<h2 id="broadcasting-indices"><span class="glyphicon glyphicon-pencil"></span>Broadcasting indices</h2>
</div>
<div class="panel-body">
<p>Predict and verify the shape of <code>y</code>:</p>
<div class="sourceCode"><pre class="sourceCode python"><code class="sourceCode python">x <span class="op">=</span> np.empty((<span class="dv">10</span>, <span class="dv">8</span>, <span class="dv">6</span>))
idx0 <span class="op">=</span> np.zeros((<span class="dv">3</span>, <span class="dv">8</span>)).astype(<span class="bu">int</span>)
idx1 <span class="op">=</span> np.zeros((<span class="dv">3</span>, <span class="dv">1</span>)).astype(<span class="bu">int</span>)
idx2 <span class="op">=</span> np.zeros((<span class="dv">1</span>, <span class="dv">1</span>)).astype(<span class="bu">int</span>)
y <span class="op">=</span> x[idx0, idx1, idx2]</code></pre></div>
</div>
</section>
<p>A lot of grid-based or network-based problems can also use broadcasting. For instance, if we want to compute the distance from the origin of points on a 10x10 grid, we can do</p>
<pre><code>>>> x = np.arange(5)
>>> y = np.arange(5)[:, np.newaxis]
>>> distance = np.sqrt(x ** 2 + y ** 2)
>>> distance
array([[ 0. , 1. , 2. , 3. , 4. ],
[ 1. , 1.41421356, 2.23606798, 3.16227766, 4.12310563],
[ 2. , 2.23606798, 2.82842712, 3.60555128, 4.47213595],
[ 3. , 3.16227766, 3.60555128, 4.24264069, 5. ],
[ 4. , 4.12310563, 4.47213595, 5. , 5.65685425]])</code></pre>
<section class="challenge panel panel-success">
<div class="panel-heading">
<h2 id="creating-a-two-dimensional-grid"><span class="glyphicon glyphicon-pencil"></span>Creating a two-dimensional grid</h2>
</div>
<div class="panel-body">
<p>What are the dimensionalities of <code>x</code>, <code>y</code> and <code>z</code> in the two cases:</p>
<pre><code>x, y = np.mgrid[:10, :5]
z = x + y</code></pre>
<p>and</p>
<pre><code>x, y = np.ogrid[:10, :5]
z = x + y</code></pre>
<p>What might be the advantage of using <code>np.ogrid</code> over <code>np.mgrid</code>?</p>
</div>
</section>
<aside class="callout panel panel-info">
<div class="panel-heading">
<h2 id="worked-example-route-66"><span class="glyphicon glyphicon-pushpin"></span>Worked example: Route 66</h2>
</div>
<div class="panel-body">
<p>Let’s construct an array of distances (in miles) between cities of Route 66: Chicago, Springfield, Saint-Louis, Tulsa, Oklahoma City, Amarillo, Santa Fe, Albuquerque, Flagstaff and Los Angeles.</p>
<pre><code>>>> mileposts = np.array([0, 198, 303, 736, 871, 1175, 1475, 1544,
... 1913, 2448])
>>> distance_array = np.abs(mileposts - mileposts[:, np.newaxis])
>>> distance_array
array([[ 0, 198, 303, 736, 871, 1175, 1475, 1544, 1913, 2448],
[ 198, 0, 105, 538, 673, 977, 1277, 1346, 1715, 2250],
[ 303, 105, 0, 433, 568, 872, 1172, 1241, 1610, 2145],
[ 736, 538, 433, 0, 135, 439, 739, 808, 1177, 1712],
[ 871, 673, 568, 135, 0, 304, 604, 673, 1042, 1577],
[1175, 977, 872, 439, 304, 0, 300, 369, 738, 1273],
[1475, 1277, 1172, 739, 604, 300, 0, 69, 438, 973],
[1544, 1346, 1241, 808, 673, 369, 69, 0, 369, 904],
[1913, 1715, 1610, 1177, 1042, 738, 438, 369, 0, 535],
[2448, 2250, 2145, 1712, 1577, 1273, 973, 904, 535, 0]])</code></pre>
<div class="figure">
<img src="fig/route66.png" alt="Distances on Route 66" />
<p class="caption">Distances on Route 66</p>
</div>
</div>
</aside>
<section class="challenge panel panel-success">
<div class="panel-heading">
<h2 id="distances"><span class="glyphicon glyphicon-pencil"></span>Distances</h2>
</div>
<div class="panel-body">
<p>Given an array of latitudes and longitudes of major European capitals calculate pairwise distances between them. Use the approximate formula:</p>
<p><span class="math display">\[D=6371.009\sqrt{(\Delta\phi)^2 + (\Delta\lambda)^2}\qquad \text{(in kilometers)},\]</span></p>
<p>where <span class="math inline">\(\Delta\phi=\phi_1-\phi_2\)</span> and <span class="math inline">\(\Delta\lambda=\lambda_1-\lambda_2\)</span> are the differences between the latitudes and longitude of two cities in radians. (<em>Hint</em>: To convert degrees to radians multiply them by <span class="math inline">\(\pi/180\)</span>).</p>
<pre><code>coords = np.array([
[ 23.71666667, 37.96666667], # Athens
[ 13.38333333, 52.51666667], # Berlin
[ -0.1275 , 51.50722222], # London
[ -3.71666667, 40.38333333], # Madrid
[ 2.3508 , 48.8567 ], # Paris
[ 12.5 , 41.9 ] # Rome
]) </code></pre>
<p>When you are done you can compare the results with a more <a href="https://en.wikipedia.org/wiki/Geographical_distance#Spherical_Earth_projected_to_a_plane">precise formula</a>:</p>
<p><span class="math display">\[D=6371.009\sqrt{(\Delta\phi)^2 + (\cos(\phi_m)\Delta\lambda)^2}\]</span></p>
<p>where <span class="math inline">\(\phi_m = (\phi_1+\phi_2) / 2\)</span> is the mean latitude.</p>
</div>
</section>
</div>
</div>
</article>
<div class="footer">
<a class="label swc-blue-bg" href="http://software-carpentry.org">Software Carpentry</a>
<a class="label swc-blue-bg" href="https://github.com/paris-swc/advanced-numpy-lesson">Source</a>
<a class="label swc-blue-bg" href="mailto:admin@software-carpentry.org">Contact</a>
<a class="label swc-blue-bg" href="LICENSE.html">License</a>
</div>
</div>
<!-- Javascript placed at the end of the document so the pages load faster -->
<script src="http://software-carpentry.org/v5/js/jquery-1.9.1.min.js"></script>
<script src="css/bootstrap/bootstrap-js/bootstrap.js"></script>
<script src='https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'></script>
</body>
</html>