1THREADING(1) User Contributed Perl Documentation THREADING(1)
2
6 PDL::Threading - Tutorial for PDL's Threading feature
7
9 One of the most powerful features of PDL is threading, which can
10 produce very compact and very fast PDL code by avoiding multiple nested
11 for loops that C and BASIC users may be familiar with. The trouble is
12 that it can take some getting used to, and new users may not appreciate
13 the benefits of threading.
14 Other vector based languages, such as MATLAB, use a subset of threading
16 techniques, but PDL shines by completely generalizing them for all
17 sorts of vector-based applications.
18
20 MATLAB typically refers to vectors, matrices, and arrays. Perl already
21 has arrays, and the terms "vector" and "matrix" typically refer to one-
22 and two-dimensional collections of data. Having no good term to
23 describe their object, PDL developers coined the term "piddle" to give
24 a name to their data type.
25 A piddle consists of a series of numbers organized as an N-dimensional
27 data set. Piddles provide efficient storage and fast computation of
28 large N-dimensional matrices. They are highly optimized for numerical
29 work.
30
32 If you have used PDL for a little while already, you may have been
33 using threading without realising it. Start the PDL shell (type
34 "perldl" or "pdl2" on a terminal). Most examples in this tutorial use
35 the PDL shell. Make sure that PDL::NiceSlice and PDL::AutoLoader are
36 enabled. For example:
37 % pdl2
39 perlDL shell v1.352
40 ...
41 ReadLines, NiceSlice, MultiLines enabled
42 ...
43 Note: AutoLoader not enabled ('use PDL::AutoLoader' recommended)
44 pdl>
46 In this example, NiceSlice was automatically enabled, but AutoLoader
48 was not. To enable it, type "use PDL::AutoLoader".
49 Let's start with a two-dimensional piddle:
51 pdl> $a = sequence(11,9)
53 pdl> p $a
54 [
55 [ 0 1 2 3 4 5 6 7 8 9 10]
56 [11 12 13 14 15 16 17 18 19 20 21]
57 [22 23 24 25 26 27 28 29 30 31 32]
58 [33 34 35 36 37 38 39 40 41 42 43]
59 [44 45 46 47 48 49 50 51 52 53 54]
60 [55 56 57 58 59 60 61 62 63 64 65]
61 [66 67 68 69 70 71 72 73 74 75 76]
62 [77 78 79 80 81 82 83 84 85 86 87]
63 [88 89 90 91 92 93 94 95 96 97 98]
64 ]
65 The "info" method gives you basic information about a piddle:
67 pdl> p $a->info
69 PDL: Double D [11,9]
70 This tells us that $a is an 11 x 9 piddle composed of double precision
72 numbers. If we wanted to add 3 to all elements in an "n x m" piddle, a
73 traditional language would use two nested for-loops:
74 # Pseudo-code. Traditional way to add 3 to an array.
76 for (x=0; x < n; x++) {
77 for (y=0; y < m; y++) {
78 a(x,y) = a(x,y) + 3
79 }
80 }
81 Note: Notice that indices start at 0, as in Perl, C and Java (and
83 unlike MATLAB and IDL).
84 But with PDL, we can just write:
86 pdl> $b = $a + 3
88 pdl> p $b
89 [
90 [ 3 4 5 6 7 8 9 10 11 12 13]
91 [ 14 15 16 17 18 19 20 21 22 23 24]
92 [ 25 26 27 28 29 30 31 32 33 34 35]
93 [ 36 37 38 39 40 41 42 43 44 45 46]
94 [ 47 48 49 50 51 52 53 54 55 56 57]
95 [ 58 59 60 61 62 63 64 65 66 67 68]
96 [ 69 70 71 72 73 74 75 76 77 78 79]
97 [ 80 81 82 83 84 85 86 87 88 89 90]
98 [ 91 92 93 94 95 96 97 98 99 100 101]
99 ]
100 This is the simplest example of threading, and it is something that all
102 numerical software tools do. The "+ 3" operation was automatically
103 applied along two dimensions. Now suppose you want to to subtract a
104 line from every row in $a:
105 pdl> $line = sequence(11)
107 pdl> p $line
108 [0 1 2 3 4 5 6 7 8 9 10]
109 pdl> $c = $a - $line
110 pdl> p $c
111 [
112 [ 0 0 0 0 0 0 0 0 0 0 0]
113 [11 11 11 11 11 11 11 11 11 11 11]
114 [22 22 22 22 22 22 22 22 22 22 22]
115 [33 33 33 33 33 33 33 33 33 33 33]
116 [44 44 44 44 44 44 44 44 44 44 44]
117 [55 55 55 55 55 55 55 55 55 55 55]
118 [66 66 66 66 66 66 66 66 66 66 66]
119 [77 77 77 77 77 77 77 77 77 77 77]
120 [88 88 88 88 88 88 88 88 88 88 88]
121 ]
122 Two things to note here: First, the value of $a is still the same. Try
124 "p $a" to check. Second, PDL automatically subtracted $line from each
125 row in $a. Why did it do that? Let's look at the dimensions of $a,
126 $line and $c:
127 pdl> p $line->info => PDL: Double D [11]
129 pdl> p $a->info => PDL: Double D [11,9]
130 pdl> p $c->info => PDL: Double D [11,9]
131 So, both $a and $line have the same number of elements in the 0th
133 dimension! What PDL then did was thread over the higher dimensions in
134 $a and repeated the same operation 9 times to all the rows on $a. This
135 is PDL threading in action.
136 What if you want to subtract $line from the first line in $a only? You
138 can do that by specifying the line explicitly:
139 pdl> $a(:,0) -= $line
141 pdl> p $a
142 [
143 [ 0 0 0 0 0 0 0 0 0 0 0]
144 [11 12 13 14 15 16 17 18 19 20 21]
145 [22 23 24 25 26 27 28 29 30 31 32]
146 [33 34 35 36 37 38 39 40 41 42 43]
147 [44 45 46 47 48 49 50 51 52 53 54]
148 [55 56 57 58 59 60 61 62 63 64 65]
149 [66 67 68 69 70 71 72 73 74 75 76]
150 [77 78 79 80 81 82 83 84 85 86 87]
151 [88 89 90 91 92 93 94 95 96 97 98]
152 ]
153 See PDL::Indexing and PDL::NiceSlice to learn more about specifying
155 subsets from piddles.
156 The true power of threading comes when you realise that the piddle can
158 have any number of dimensions! Let's make a 4 dimensional piddle:
159 pdl> $piddle_4D = sequence(11,3,7,2)
161 pdl> $c = $piddle_4D - $line
162 Now $c is a piddle of the same dimension as $piddle_4D.
164 pdl> p $piddle_4D->info => PDL: Double D [11,3,7,2]
166 pdl> p $c->info => PDL: Double D [11,3,7,2]
167 This time PDL has threaded over three higher dimensions automatically,
169 subtracting $line all the way.
170 But, maybe you don't want to subtract from the rows (dimension 0), but
172 from the columns (dimension 1). How do I subtract a column of numbers
173 from each column in $a?
174 pdl> $cols = sequence(9)
176 pdl> p $a->info => PDL: Double D [11,9]
177 pdl> p $cols->info => PDL: Double D [9]
178 Naturally, we can't just type "$a - $cols". The dimensions don't match:
180 pdl> p $a - $cols
182 PDL: PDL::Ops::minus(a,b,c): Parameter 'b'
183 PDL: Mismatched implicit thread dimension 0: should be 11, is 9
184 How do we tell PDL that we want to subtract from dimension 1 instead?
186
188 There are many PDL functions that let you rearrange the dimensions of
189 PDL arrays. They are mostly covered in PDL::Slices. The three most
190 common ones are:
191 xchg
193 mv
194 reorder
195 Method: "xchg"
197 The "xchg" method "exchanges" two dimensions in a piddle:
198 pdl> $a = sequence(6,7,8,9)
200 pdl> $a_xchg = $a->xchg(0,3)
201 pdl> p $a->info => PDL: Double D [6,7,8,9]
203 pdl> p $a_xchg->info => PDL: Double D [9,7,8,6]
204 | |
205 V V
206 (dim 0) (dim 3)
207 Notice that dimensions 0 and 3 were exchanged without affecting the
209 other dimensions. Notice also that "xchg" does not alter $a. The
210 original variable $a remains untouched.
211 Method: "mv"
213 The "mv" method "moves" one dimension, in a piddle, shifting other
214 dimensions as necessary.
215 pdl> $a = sequence(6,7,8,9) (dim 0)
217 pdl> $a_mv = $a->mv(0,3) |
218 pdl> V _____
219 pdl> p $a->info => PDL: Double D [6,7,8,9]
220 pdl> p $a_mv->info => PDL: Double D [7,8,9,6]
221 ----- |
222 V
223 (dim 3)
224 Notice that when dimension 0 was moved to position 3, all the other
226 dimensions had to be shifted as well. Notice also that "mv" does not
227 alter $a. The original variable $a remains untouched.
228 Method: "reorder"
230 The "reorder" method is a generalization of the "xchg" and "mv"
231 methods. It "reorders" the dimensions in any way you specify:
232 pdl> $a = sequence(6,7,8,9)
234 pdl> $a_reorder = $a->reorder(3,0,2,1)
235 pdl>
236 pdl> p $a->info => PDL: Double D [6,7,8,9]
237 pdl> p $a_reorder->info => PDL: Double D [9,6,8,7]
238 | | | |
239 V V v V
240 dimensions: 0 1 2 3
241 Notice what happened. When we wrote "reorder(3,0,2,1)" we instructed
243 PDL to:
244 * Put dimension 3 first.
246 * Put dimension 0 next.
247 * Put dimension 2 next.
248 * Put dimension 1 next.
249 When you use the "reorder" method, all the dimensions are shuffled.
251 Notice that "reorder" does not alter $a. The original variable $a
252 remains untouched.
253
255 Linking
256 By default, piddles are linked together so that changes on one will go
257 back and affect the original as well.
258 pdl> $a = sequence(4,5)
260 pdl> $a_xchg = $a->xchg(1,0)
261 Here, $a_xchg is not a separate object. It is merely a different way of
263 looking at $a. Any change in $a_xchg will appear in $a as well.
264 pdl> p $a
266 [
267 [ 0 1 2 3]
268 [ 4 5 6 7]
269 [ 8 9 10 11]
270 [12 13 14 15]
271 [16 17 18 19]
272 ]
273 pdl> $a_xchg += 3
274 pdl> p $a
275 [
276 [ 3 4 5 6]
277 [ 7 8 9 10]
278 [11 12 13 14]
279 [15 16 17 18]
280 [19 20 21 22]
281 ]
282 Assignment
284 Some times, linking is not the behaviour you want. If you want to make
285 the piddles independent, use the "copy" method:
286 pdl> $a = sequence(4,5)
288 pdl> $a_xchg = $a->copy->xchg(1,0)
289 Now $a and $a_xchg are completely separate objects:
291 pdl> p $a
293 [
294 [ 0 1 2 3]
295 [ 4 5 6 7]
296 [ 8 9 10 11]
297 [12 13 14 15]
298 [16 17 18 19]
299 ]
300 pdl> $a_xchg += 3
301 pdl> p $a
302 [
303 [ 0 1 2 3]
304 [ 4 5 6 7]
305 [ 8 9 10 11]
306 [12 13 14 15]
307 [16 17 18 19]
308 ]
309 pdl> $a_xchg
310 [
311 [ 3 7 11 15 19]
312 [ 4 8 12 16 20]
313 [ 5 9 13 17 21]
314 [ 6 10 14 18 22]
315 ]
316
318 Now we are ready to solve the problem that motivated this whole
319 discussion:
320 pdl> $a = sequence(11,9)
322 pdl> $cols = sequence(9)
323 pdl>
324 pdl> p $a->info => PDL: Double D [11,9]
325 pdl> p $cols->info => PDL: Double D [9]
326 How do we tell PDL to subtract $cols along dimension 1 instead of
328 dimension 0? The simplest way is to use the "xchg" method and rely on
329 PDL linking:
330 pdl> p $a
332 [
333 [ 0 1 2 3 4 5 6 7 8 9 10]
334 [11 12 13 14 15 16 17 18 19 20 21]
335 [22 23 24 25 26 27 28 29 30 31 32]
336 [33 34 35 36 37 38 39 40 41 42 43]
337 [44 45 46 47 48 49 50 51 52 53 54]
338 [55 56 57 58 59 60 61 62 63 64 65]
339 [66 67 68 69 70 71 72 73 74 75 76]
340 [77 78 79 80 81 82 83 84 85 86 87]
341 [88 89 90 91 92 93 94 95 96 97 98]
342 ]
343 pdl> $a->xchg(1,0) -= $cols
344 pdl> p $a
345 [
346 [ 0 1 2 3 4 5 6 7 8 9 10]
347 [10 11 12 13 14 15 16 17 18 19 20]
348 [20 21 22 23 24 25 26 27 28 29 30]
349 [30 31 32 33 34 35 36 37 38 39 40]
350 [40 41 42 43 44 45 46 47 48 49 50]
351 [50 51 52 53 54 55 56 57 58 59 60]
352 [60 61 62 63 64 65 66 67 68 69 70]
353 [70 71 72 73 74 75 76 77 78 79 80]
354 [80 81 82 83 84 85 86 87 88 89 90]
355 ]
356 General Strategy:
358 Move the dimensions you want to operate on to the start of your
359 piddle's dimension list. Then let PDL thread over the higher
360 dimensions.
361
363 Okay, enough theory. Let's do something a bit more interesting: We'll
364 write Conway's Game of Life in PDL and see how powerful PDL can be!
365 The Game of Life is a simulation run on a big two dimensional grid.
367 Each cell in the grid can either be alive or dead (represented by 1 or
368 0). The next generation of cells in the grid is calculated with simple
369 rules according to the number of living cells in it's immediate
370 neighbourhood:
371 1) If an empty cell has exactly three neighbours, a living cell is
373 generated.
374 2) If a living cell has less than two neighbours, it dies of
376 overfeeding.
377 3) If a living cell has 4 or more neighbours, it dies from starvation.
379 Only the first generation of cells is determined by the programmer.
381 After that, the simulation runs completely according to these rules. To
382 calculate the next generation, you need to look at each cell in the 2D
383 field (requiring two loops), calculate the number of live cells
384 adjacent to this cell (requiring another two loops) and then fill the
385 next generation grid.
386 Classical implementation
388 Here's a classic way of writing this program in Perl. We only use PDL
389 for addressing individual cells.
390 #!/usr/bin/perl -w
392 use PDL;
393 use PDL::NiceSlice;
394 # Make a board for the game of life.
396 my $nx = 20;
397 my $ny = 20;
398 # Current generation.
400 my $a = zeroes($nx, $ny);
401 # Next generation.
403 my $n = zeroes($nx, $ny);
404 # Put in a simple glider.
406 $a(1:3,1:3) .= pdl ( [1,1,1],
407 [0,0,1],
408 [0,1,0] );
409 for (my $i = 0; $i < 100; $i++) {
411 $n = zeroes($nx, $ny);
412 $new_a = $a->copy;
413 for ($x = 0; $x < $nx; $x++) {
414 for ($y = 0; $y < $ny; $y++) {
415 # For each cell, look at the surrounding neighbours.
417 for ($dx = -1; $dx <= 1; $dx++) {
418 for ($dy = -1; $dy <= 1; $dy++) {
419 $px = $x + $dx;
420 $py = $y + $dy;
421 # Wrap around at the edges.
423 if ($px < 0) {$px = $nx-1};
424 if ($py < 0) {$py = $ny-1};
425 if ($px >= $nx) {$px = 0};
426 if ($py >= $ny) {$py = 0};
427 $n($x,$y) .= $n($x,$y) + $a($px,$py);
429 }
430 }
431 # Do not count the central cell itself.
432 $n($x,$y) -= $a($x,$y);
433 # Work out if cell lives or dies:
435 # Dead cell lives if n = 3
436 # Live cell dies if n is not 2 or 3
437 if ($a($x,$y) == 1) {
438 if ($n($x,$y) < 2) {$new_a($x,$y) .= 0};
439 if ($n($x,$y) > 3) {$new_a($x,$y) .= 0};
440 } else {
441 if ($n($x,$y) == 3) {$new_a($x,$y) .= 1}
442 }
443 }
444 }
445 print $a;
447 $a = $new_a;
449 }
450 If you run this, you will see a small glider crawl diagonally across
452 the grid of zeroes. On my machine, it prints out a couple of
453 generations per second.
454 Threaded PDL implementation
456 And here's the threaded version in PDL. Just four lines of PDL code,
457 and one of those is printing out the latest generation!
458 #!/usr/bin/perl -w
460 use PDL;
461 use PDL::NiceSlice;
462 my $a = zeroes(20,20);
464 # Put in a simple glider.
466 $a(1:3,1:3) .= pdl ( [1,1,1],
467 [0,0,1],
468 [0,1,0] );
469 my $n;
471 for (my $i = 0; $i < 100; $i++) {
472 # Calculate the number of neighbours per cell.
473 $n = $a->range(ndcoords($a)-1,3,"periodic")->reorder(2,3,0,1);
474 $n = $n->sumover->sumover - $a;
475 # Calculate the next generation.
477 $a = ((($n == 2) + ($n == 3))* $a) + (($n==3) * !$a);
478 print $a;
480 }
481 The threaded PDL version is much faster:
483 Classical => 32.79 seconds.
485 Threaded => 0.41 seconds.
486 Explanation
488 How does the threaded version work?
489 There are many PDL functions designed to help you carry out PDL
491 threading. In this example, the key functions are:
492 Method: "range"
494 At the simplest level, the "range" method is a different way to select
496 a portion of a piddle. Instead of using the "$a(2,3)" notation, we use
497 another piddle.
498 pdl> $a = sequence(6,7)
500 pdl> p $a
501 [
502 [ 0 1 2 3 4 5]
503 [ 6 7 8 9 10 11]
504 [12 13 14 15 16 17]
505 [18 19 20 21 22 23]
506 [24 25 26 27 28 29]
507 [30 31 32 33 34 35]
508 [36 37 38 39 40 41]
509 ]
510 pdl> p $a->range( pdl [1,2] )
511 13
512 pdl> p $a(1,2)
513 [
514 [13]
515 ]
516 At this point, the "range" method looks very similar to a regular PDL
518 slice. But the "range" method is more general. For example, you can
519 select several components at once:
520 pdl> $index = pdl [ [1,2],[2,3],[3,4],[4,5] ]
522 pdl> p $a->range( $index )
523 [13 20 27 34]
524 Additionally, "range" takes a second parameter which determines the
526 size of the chunk to return:
527 pdl> $size = 3
529 pdl> p $a->range( pdl([1,2]) , $size )
530 [
531 [13 14 15]
532 [19 20 21]
533 [25 26 27]
534 ]
535 We can use this to select one or more 3x3 boxes.
537 Finally, "range" can take a third parameter called the "boundary"
539 condition. It tells PDL what to do if the size box you request goes
540 beyond the edge of the piddle. We won't go over all the options. We'll
541 just say that the option "periodic" means that the piddle "wraps
542 around". For example:
543 pdl> p $a
545 [
546 [ 0 1 2 3 4 5]
547 [ 6 7 8 9 10 11]
548 [12 13 14 15 16 17]
549 [18 19 20 21 22 23]
550 [24 25 26 27 28 29]
551 [30 31 32 33 34 35]
552 [36 37 38 39 40 41]
553 ]
554 pdl> $size = 3
555 pdl> p $a->range( pdl([4,2]) , $size , "periodic" )
556 [
557 [16 17 12]
558 [22 23 18]
559 [28 29 24]
560 ]
561 pdl> p $a->range( pdl([5,2]) , $size , "periodic" )
562 [
563 [17 12 13]
564 [23 18 19]
565 [29 24 25]
566 ]
567 Notice how the box wraps around the boundary of the piddle.
569 Method: "ndcoords"
571 The "ndcoords" method is a convenience method that returns an
573 enumerated list of coordinates suitable for use with the "range"
574 method.
575 pdl> p $piddle = sequence(3,3)
577 [
578 [0 1 2]
579 [3 4 5]
580 [6 7 8]
581 ]
582 pdl> p ndcoords($piddle)
583 [
584 [
585 [0 0]
586 [1 0]
587 [2 0]
588 ]
589 [
590 [0 1]
591 [1 1]
592 [2 1]
593 ]
594 [
595 [0 2]
596 [1 2]
597 [2 2]
598 ]
599 ]
600 This can be a little hard to read. Basically it's saying that the
602 coordinates for every element in $piddle is given by:
603 (0,0) (1,0) (2,0)
605 (1,0) (1,1) (2,1)
606 (2,0) (2,1) (2,2)
607 Combining "range" and "ndcoords"
609 What really matters is that "ndcoords" is designed to work together
611 with "range", with no $size parameter, you get the same piddle back.
612 pdl> p $piddle
614 [
615 [0 1 2]
616 [3 4 5]
617 [6 7 8]
618 ]
619 pdl> p $piddle->range( ndcoords($piddle) )
620 [
621 [0 1 2]
622 [3 4 5]
623 [6 7 8]
624 ]
625 Why would this be useful? Because now we can ask for a series of
627 "boxes" for the entire piddle. For example, 2x2 boxes:
628 pdl> p $piddle->range( ndcoords($piddle) , 2 , "periodic" )
630 The output of this function is difficult to read because the "boxes"
632 along the last two dimension. We can make the result more readable by
633 rearranging the dimensions:
634 pdl> p $piddle->range( ndcoords($piddle) , 2 , "periodic" )->reorder(2,3,0,1)
636 [
637 [
638 [
639 [0 1]
640 [3 4]
641 ]
642 [
643 [1 2]
644 [4 5]
645 ]
646 ...
647 ]
648 Here you can see more clearly that
650 [0 1]
652 [3 4]
653 Is the 2x2 box starting with the (0,0) element of $piddle.
655 We are not done yet. For the game of life, we want 3x3 boxes from $a:
657 pdl> p $a
659 [
660 [ 0 1 2 3 4 5]
661 [ 6 7 8 9 10 11]
662 [12 13 14 15 16 17]
663 [18 19 20 21 22 23]
664 [24 25 26 27 28 29]
665 [30 31 32 33 34 35]
666 [36 37 38 39 40 41]
667 ]
668 pdl> p $a->range( ndcoords($a) , 3 , "periodic" )->reorder(2,3,0,1)
669 [
670 [
671 [
672 [ 0 1 2]
673 [ 6 7 8]
674 [12 13 14]
675 ]
676 ...
677 ]
678 We can confirm that this is the 3x3 box starting with the (0,0) element
680 of $a. But there is one problem. We actually want the 3x3 box to be
681 centered on (0,0). That's not a problem. Just subtract 1 from all the
682 coordinates in "ndcoords($a)". Remember that the "periodic" option
683 takes care of making everything wrap around.
684 pdl> p $a->range( ndcoords($a) - 1 , 3 , "periodic" )->reorder(2,3,0,1)
686 [
687 [
688 [
689 [41 36 37]
690 [ 5 0 1]
691 [11 6 7]
692 ]
693 [
694 [36 37 38]
695 [ 0 1 2]
696 [ 6 7 8]
697 ]
698 ...
699 Now we see a 3x3 box with the (0,0) element in the centre of the box.
701 Method: "sumover"
703 The "sumover" method adds along only the first dimension. If we apply
705 it twice, we will be adding all the elements of each 3x3 box.
706 pdl> $n = $a->range(ndcoords($a)-1,3,"periodic")->reorder(2,3,0,1)
708 pdl> p $n
709 [
710 [
711 [
712 [41 36 37]
713 [ 5 0 1]
714 [11 6 7]
715 ]
716 [
717 [36 37 38]
718 [ 0 1 2]
719 [ 6 7 8]
720 ]
721 ...
722 pdl> p $n->sumover->sumover
723 [
724 [144 135 144 153 162 153]
725 [ 72 63 72 81 90 81]
726 [126 117 126 135 144 135]
727 [180 171 180 189 198 189]
728 [234 225 234 243 252 243]
729 [288 279 288 297 306 297]
730 [216 207 216 225 234 225]
731 ]
732 Use a calculator to confirm that 144 is the sum of all the elements in
734 the first 3x3 box and 135 is the sum of all the elements in the second
735 3x3 box.
736 Counting neighbours
738 We are almost there!
740 Adding up all the elements in a 3x3 box is not quite what we want. We
742 don't want to count the center box. Fortunately, this is an easy fix:
743 pdl> p $n->sumover->sumover - $a
745 [
746 [144 134 142 150 158 148]
747 [ 66 56 64 72 80 70]
748 [114 104 112 120 128 118]
749 [162 152 160 168 176 166]
750 [210 200 208 216 224 214]
751 [258 248 256 264 272 262]
752 [180 170 178 186 194 184]
753 ]
754 When applied to Conway's Game of Life, this will tell us how many
756 living neighbours each cell has:
757 pdl> $a = zeroes(10,10)
759 pdl> $a(1:3,1:3) .= pdl ( [1,1,1],
760 ..( > [0,0,1],
761 ..( > [0,1,0] )
762 pdl> p $a
763 [
764 [0 0 0 0 0 0 0 0 0 0]
765 [0 1 1 1 0 0 0 0 0 0]
766 [0 0 0 1 0 0 0 0 0 0]
767 [0 0 1 0 0 0 0 0 0 0]
768 [0 0 0 0 0 0 0 0 0 0]
769 [0 0 0 0 0 0 0 0 0 0]
770 [0 0 0 0 0 0 0 0 0 0]
771 [0 0 0 0 0 0 0 0 0 0]
772 [0 0 0 0 0 0 0 0 0 0]
773 [0 0 0 0 0 0 0 0 0 0]
774 ]
775 pdl> $n = $a->range(ndcoords($a)-1,3,"periodic")->reorder(2,3,0,1)
776 pdl> $n = $n->sumover->sumover - $a
777 pdl> p $n
778 [
779 [1 2 3 2 1 0 0 0 0 0]
780 [1 1 3 2 2 0 0 0 0 0]
781 [1 3 5 3 2 0 0 0 0 0]
782 [0 1 1 2 1 0 0 0 0 0]
783 [0 1 1 1 0 0 0 0 0 0]
784 [0 0 0 0 0 0 0 0 0 0]
785 [0 0 0 0 0 0 0 0 0 0]
786 [0 0 0 0 0 0 0 0 0 0]
787 [0 0 0 0 0 0 0 0 0 0]
788 [0 0 0 0 0 0 0 0 0 0]
789 ]
790 For example, this tells us that cell (0,0) has 1 living neighbour,
792 while cell (2,2) has 5 living neighbours.
793 Calculating the next generation
795 At this point, the variable $n has the number of living neighbours for
797 every cell. Now we apply the rules of the game of life to calculate the
798 next generation.
799 If an empty cell has exactly three neighbours, a living cell is
801 generated.
802 Get a list of cells that have exactly three neighbours:
803 pdl> p ($n == 3)
805 [
806 [0 0 1 0 0 0 0 0 0 0]
807 [0 0 1 0 0 0 0 0 0 0]
808 [0 1 0 1 0 0 0 0 0 0]
809 [0 0 0 0 0 0 0 0 0 0]
810 [0 0 0 0 0 0 0 0 0 0]
811 [0 0 0 0 0 0 0 0 0 0]
812 [0 0 0 0 0 0 0 0 0 0]
813 [0 0 0 0 0 0 0 0 0 0]
814 [0 0 0 0 0 0 0 0 0 0]
815 [0 0 0 0 0 0 0 0 0 0]
816 ]
817 Get a list of empty cells that have exactly three neighbours:
819 pdl> p ($n == 3) * !$a
821 If a living cell has less than 2 or more than 3 neighbours, it dies.
823 Get a list of cells that have exactly 2 or 3 neighbours:
824 pdl> p (($n == 2) + ($n == 3))
826 [
827 [0 1 1 1 0 0 0 0 0 0]
828 [0 0 1 1 1 0 0 0 0 0]
829 [0 1 0 1 1 0 0 0 0 0]
830 [0 0 0 1 0 0 0 0 0 0]
831 [0 0 0 0 0 0 0 0 0 0]
832 [0 0 0 0 0 0 0 0 0 0]
833 [0 0 0 0 0 0 0 0 0 0]
834 [0 0 0 0 0 0 0 0 0 0]
835 [0 0 0 0 0 0 0 0 0 0]
836 [0 0 0 0 0 0 0 0 0 0]
837 ]
838 Get a list of living cells that have exactly 2 or 3 neighbours:
840 pdl> p (($n == 2) + ($n == 3)) * $a
842 Putting it all together, the next generation is:
844 pdl> $a = ((($n == 2) + ($n == 3)) * $a) + (($n == 3) * !$a)
846 pdl> p $a
847 [
848 [0 0 1 0 0 0 0 0 0 0]
849 [0 0 1 1 0 0 0 0 0 0]
850 [0 1 0 1 0 0 0 0 0 0]
851 [0 0 0 0 0 0 0 0 0 0]
852 [0 0 0 0 0 0 0 0 0 0]
853 [0 0 0 0 0 0 0 0 0 0]
854 [0 0 0 0 0 0 0 0 0 0]
855 [0 0 0 0 0 0 0 0 0 0]
856 [0 0 0 0 0 0 0 0 0 0]
857 [0 0 0 0 0 0 0 0 0 0]
858 ]
859 Bonus feature: Graphics!
861 If you have PDL::Graphics::TriD installed, you can make a graphical
862 version of the program by just changing three lines:
863 #!/usr/bin/perl
865 use PDL;
866 use PDL::NiceSlice;
867 use PDL::Graphics::TriD;
868 my $a = zeroes(20,20);
870 # Put in a simple glider.
872 $a(1:3,1:3) .= pdl ( [1,1,1],
873 [0,0,1],
874 [0,1,0] );
875 my $n;
877 for (my $i = 0; $i < 100; $i++) {
878 # Calculate the number of neighbours per cell.
879 $n = $a->range(ndcoords($a)-1,3,"periodic")->reorder(2,3,0,1);
880 $n = $n->sumover->sumover - $a;
881 # Calculate the next generation.
883 $a = ((($n == 2) + ($n == 3))* $a) + (($n==3) * !$a);
884 # Display.
886 nokeeptwiddling3d();
887 imagrgb [$a];
888 }
889 But if we really want to see something interesting, we should make a
891 few more changes:
892 1) Start with a random collection of 1's and 0's.
894 2) Make the grid larger.
896 3) Add a small timeout so we can see the game evolve better.
898 4) Use a while loop so that the program can run as long as it needs to.
900 #!/usr/bin/perl
902 use PDL;
903 use PDL::NiceSlice;
904 use PDL::Graphics::TriD;
905 use Time::HiRes qw(usleep);
906 my $a = random(100,100);
908 $a = ($a < 0.5);
909 my $n;
911 while (1) {
912 # Calculate the number of neighbours per cell.
913 $n = $a->range(ndcoords($a)-1,3,"periodic")->reorder(2,3,0,1);
914 $n = $n->sumover->sumover - $a;
915 # Calculate the next generation.
917 $a = ((($n == 2) + ($n == 3))* $a) + (($n==3) * !$a);
918 # Display.
920 nokeeptwiddling3d();
921 imagrgb [$a];
922 # Sleep for 0.1 seconds.
924 usleep(100000);
925 }
926
928 The general strategy is: Move the dimensions you want to operate on to
929 the start of your piddle's dimension list. Then let PDL thread over the
930 higher dimensions.
931 Threading is a powerful tool that helps eliminate for-loops and can
933 make your code more concise. Hopefully this tutorial has shown why it
934 is worth getting to grips with threading in PDL.
935
937 Copyright 2010 Matthew Kenworthy (kenworthy@strw.leidenuniv.nl) and
938 Daniel Carrera (dcarrera@gmail.com). You can distribute and/or modify
939 this document under the same terms as the current Perl license.
940 See: http://dev.perl.org/licenses/
942perl v5.28.1 2018-05-05 THREADING(1)