|
243 | 243 | }, |
244 | 244 | { |
245 | 245 | "cell_type": "code", |
246 | | - "execution_count": 303, |
| 246 | + "execution_count": 1, |
247 | 247 | "metadata": {}, |
248 | 248 | "outputs": [], |
249 | 249 | "source": [ |
250 | | - "import random, numpy as np\n", |
251 | | - "\n", |
| 250 | + "import random, numpy as np" |
| 251 | + ] |
| 252 | + }, |
| 253 | + { |
| 254 | + "cell_type": "code", |
| 255 | + "execution_count": 21, |
| 256 | + "metadata": {}, |
| 257 | + "outputs": [], |
| 258 | + "source": [ |
| 259 | + "# %load -s random_walk1D, src-diffu/random_walk.py\n", |
252 | 260 | "def random_walk1D(x0, N, p, random=random):\n", |
253 | | - " \"\"\"1D random walk with 1 particle.\"\"\"\n", |
| 261 | + " \"\"\"1D random walk with 1 particle and N moves.\"\"\"\n", |
| 262 | + " # random is argument so we can use np.random instead\n", |
| 263 | + " # and use it for testing\n", |
| 264 | + "\n", |
254 | 265 | " # Store position in step k in position[k]\n", |
255 | | - " position = np.zeros(N)\n", |
| 266 | + " position = np.zeros(N+1)\n", |
256 | 267 | " position[0] = x0\n", |
257 | 268 | " current_pos = x0\n", |
258 | | - " for k in range(N-1):\n", |
| 269 | + " for k in range(N):\n", |
259 | 270 | " r = random.uniform(0, 1)\n", |
260 | 271 | " if r <= p:\n", |
261 | 272 | " current_pos -= 1\n", |
262 | 273 | " else:\n", |
263 | 274 | " current_pos += 1\n", |
264 | 275 | " position[k+1] = current_pos\n", |
265 | | - " return position" |
| 276 | + " return position\n" |
266 | 277 | ] |
267 | 278 | }, |
268 | 279 | { |
|
287 | 298 | }, |
288 | 299 | { |
289 | 300 | "cell_type": "code", |
290 | | - "execution_count": 295, |
| 301 | + "execution_count": 22, |
291 | 302 | "metadata": {}, |
292 | 303 | "outputs": [], |
293 | 304 | "source": [ |
|
306 | 317 | }, |
307 | 318 | { |
308 | 319 | "cell_type": "code", |
309 | | - "execution_count": 296, |
| 320 | + "execution_count": 23, |
310 | 321 | "metadata": {}, |
311 | 322 | "outputs": [], |
312 | 323 | "source": [ |
|
315 | 326 | }, |
316 | 327 | { |
317 | 328 | "cell_type": "code", |
318 | | - "execution_count": 304, |
| 329 | + "execution_count": 24, |
319 | 330 | "metadata": {}, |
320 | 331 | "outputs": [], |
321 | 332 | "source": [ |
| 333 | + "# %load -s random_walk1D_vec, src-diffu/random_walk.py\n", |
322 | 334 | "def random_walk1D_vec(x0, N, p, random=random):\n", |
323 | | - " \"\"\"1D random walk with 1 particle.\"\"\"\n", |
| 335 | + " \"\"\"Vectorized version of random_walk1D.\"\"\"\n", |
324 | 336 | " # Create and initialise r with dimension x_d\n", |
325 | 337 | " x_d = Dimension(name='x_d', spacing=Constant('h_x'))\n", |
326 | 338 | " r = TimeFunction(name='r', dimensions=(x_d,), shape=(N+1,))\n", |
|
334 | 346 | " eq = Eq(r.forward, r + steps)\n", |
335 | 347 | " op = Operator(eq)\n", |
336 | 348 | " op.apply()\n", |
337 | | - " return r.data" |
| 349 | + " return r.data\n" |
338 | 350 | ] |
339 | 351 | }, |
340 | 352 | { |
|
671 | 683 | }, |
672 | 684 | { |
673 | 685 | "cell_type": "code", |
674 | | - "execution_count": 283, |
| 686 | + "execution_count": 25, |
675 | 687 | "metadata": {}, |
676 | 688 | "outputs": [], |
677 | 689 | "source": [ |
678 | | - "from devito import Grid\n", |
679 | | - "\n", |
680 | | - "def random_walks1D(x0, N, p, num_walks=1, num_times=1, random=random):\n", |
| 690 | + "# %load -s random_walks1D_vec, src-diffu/random_walk.py\n", |
| 691 | + "def random_walks1D_vec(x0, N, p, num_walks=1, num_times=1, random=random):\n", |
681 | 692 | " \"\"\"Vectorized version of random_walks1D.\"\"\"\n", |
682 | | - " position = np.zeros(N+1) # Accumulated positions\n", |
| 693 | + " position = np.zeros(N+1) # Accumulated positions\n", |
683 | 694 | " position2 = np.zeros(N+1) # Accumulated positions**2\n", |
684 | 695 | " # Histogram at num_times selected time points\n", |
685 | 696 | " pos_hist = np.zeros((num_walks, num_times))\n", |
|
711 | 722 | " \n", |
712 | 723 | " pos_hist[:,:] = np.transpose(r.data[pos_hist_times,:])\n", |
713 | 724 | " \n", |
714 | | - " return position, position2, pos_hist, np.array(pos_hist_times)" |
| 725 | + " return position, position2, pos_hist, np.array(pos_hist_times)\n" |
715 | 726 | ] |
716 | 727 | }, |
717 | 728 | { |
|
754 | 765 | }, |
755 | 766 | { |
756 | 767 | "cell_type": "code", |
757 | | - "execution_count": 297, |
| 768 | + "execution_count": 26, |
758 | 769 | "metadata": {}, |
759 | 770 | "outputs": [], |
760 | 771 | "source": [ |
|
779 | 790 | }, |
780 | 791 | { |
781 | 792 | "cell_type": "code", |
782 | | - "execution_count": 298, |
| 793 | + "execution_count": 27, |
783 | 794 | "metadata": {}, |
784 | 795 | "outputs": [ |
785 | 796 | { |
|
814 | 825 | "run, and it can be very difficult to spot bugs through averaged statistical\n", |
815 | 826 | "quantities.\n", |
816 | 827 | "\n", |
| 828 | + "A full test suite for this notebook, along with the rest of the functions for random walks, can be found in [`random_walk.py`](src-diffu/random_walk.py).\n", |
| 829 | + "\n", |
817 | 830 | "## Demonstration of multiple walks\n", |
818 | 831 | "\n", |
819 | 832 | "Assuming now that the code works, we can just scale up the number of\n", |
|
1084 | 1097 | }, |
1085 | 1098 | { |
1086 | 1099 | "cell_type": "code", |
1087 | | - "execution_count": 299, |
| 1100 | + "execution_count": null, |
1088 | 1101 | "metadata": {}, |
1089 | 1102 | "outputs": [], |
1090 | 1103 | "source": [ |
| 1104 | + "# %load -s random_walk2D, src-diffu/random_walk.py\n", |
1091 | 1105 | "def random_walk2D(x0, N, p, random=random):\n", |
1092 | 1106 | " \"\"\"2D random walk with 1 particle and N moves: N, E, W, S.\"\"\"\n", |
1093 | 1107 | " # Store position in step k in position[k]\n", |
1094 | 1108 | " d = len(x0)\n", |
1095 | 1109 | " position = np.zeros((N+1, d))\n", |
1096 | | - " position[0,:] = x0\n", |
| 1110 | + " position[0, :] = x0\n", |
1097 | 1111 | " current_pos = np.array(x0, dtype=float)\n", |
1098 | 1112 | " for k in range(N):\n", |
1099 | 1113 | " r = random.uniform(0, 1)\n", |
|
1105 | 1119 | " current_pos += np.array([0, -1]) # Move south\n", |
1106 | 1120 | " else:\n", |
1107 | 1121 | " current_pos += np.array([-1, 0]) # Move west\n", |
1108 | | - " position[k+1,:] = current_pos\n", |
1109 | | - " return position" |
| 1122 | + " position[k+1, :] = current_pos\n", |
| 1123 | + " return position\n" |
1110 | 1124 | ] |
1111 | 1125 | }, |
1112 | 1126 | { |
|
1144 | 1158 | }, |
1145 | 1159 | { |
1146 | 1160 | "cell_type": "code", |
1147 | | - "execution_count": 301, |
| 1161 | + "execution_count": null, |
1148 | 1162 | "metadata": {}, |
1149 | 1163 | "outputs": [], |
1150 | 1164 | "source": [ |
| 1165 | + "# %load -s random_walkdD, src-diffu/random_walk.py\n", |
1151 | 1166 | "def random_walkdD(x0, N, p, random=random):\n", |
1152 | 1167 | " \"\"\"Any-D (diagonal) random walk with 1 particle and N moves.\"\"\"\n", |
1153 | 1168 | " # Store position in step k in position[k]\n", |
1154 | 1169 | " d = len(x0)\n", |
1155 | 1170 | " position = np.zeros((N+1, d))\n", |
1156 | | - " position[0,:] = x0\n", |
| 1171 | + " position[0, :] = x0\n", |
1157 | 1172 | " current_pos = np.array(x0, dtype=float)\n", |
1158 | 1173 | " for k in range(N):\n", |
1159 | 1174 | " for i in range(d):\n", |
|
1162 | 1177 | " current_pos[i] -= 1\n", |
1163 | 1178 | " else:\n", |
1164 | 1179 | " current_pos[i] += 1\n", |
1165 | | - " position[k+1,:] = current_pos\n", |
1166 | | - " return position" |
| 1180 | + " position[k+1, :] = current_pos\n", |
| 1181 | + " return position\n" |
1167 | 1182 | ] |
1168 | 1183 | }, |
1169 | 1184 | { |
|
1175 | 1190 | }, |
1176 | 1191 | { |
1177 | 1192 | "cell_type": "code", |
1178 | | - "execution_count": 305, |
| 1193 | + "execution_count": null, |
1179 | 1194 | "metadata": {}, |
1180 | 1195 | "outputs": [], |
1181 | 1196 | "source": [ |
| 1197 | + "# %load -s random_walkdD_vec, src-diffu/random_walk.py\n", |
1182 | 1198 | "def random_walkdD_vec(x0, N, p, random=random):\n", |
| 1199 | + " \"\"\"Vectorized version of random_walkdD.\"\"\"\n", |
| 1200 | + "\n", |
1183 | 1201 | " # Here, x0 is an array of initial values in each spatial dimension\n", |
1184 | 1202 | " d = len(x0)\n", |
1185 | 1203 | " \n", |
|
1197 | 1215 | " eq = Eq(r.forward, r + steps)\n", |
1198 | 1216 | " op = Operator(eq)\n", |
1199 | 1217 | " op.apply()\n", |
1200 | | - " return r.data" |
| 1218 | + " return r.data\n" |
1201 | 1219 | ] |
1202 | 1220 | }, |
1203 | 1221 | { |
|
1221 | 1239 | }, |
1222 | 1240 | { |
1223 | 1241 | "cell_type": "code", |
1224 | | - "execution_count": 225, |
| 1242 | + "execution_count": null, |
1225 | 1243 | "metadata": {}, |
1226 | 1244 | "outputs": [], |
1227 | 1245 | "source": [ |
1228 | | - "def random_walksdD_vec(x0, N, p, num_walks=1, num_times=1, random=random):\n", |
| 1246 | + "# %load -s random_walksdD_vec, src-diffu/random_walk.py\n", |
| 1247 | + "def random_walksdD_vec(x0, N, p, num_walks=1, num_times=1, random=np.random):\n", |
1229 | 1248 | " \"\"\"Vectorized version of random_walksdD.\"\"\"\n", |
1230 | 1249 | " d = len(x0)\n", |
1231 | 1250 | " position = np.zeros((N+1, d)) # Accumulated positions\n", |
|
1261 | 1280 | " for k in range(num_walks):\n", |
1262 | 1281 | " pos_hist[k,:] = r.data[pos_hist_times,k]\n", |
1263 | 1282 | " \n", |
1264 | | - " return position, position2, pos_hist, np.array(pos_hist_times)" |
| 1283 | + " return position, position2, pos_hist, np.array(pos_hist_times)\n" |
1265 | 1284 | ] |
1266 | 1285 | }, |
1267 | 1286 | { |
|
0 commit comments