1THREADING(1)          User Contributed Perl Documentation         THREADING(1)
2
3
4

NAME

6       PDL::Threading - Tutorial for PDL's Threading feature
7

INTRODUCTION

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
15       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

TERMINOLOGY: PIDDLE

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
26       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

THINKING IN TERMS OF THREADING

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
38         % pdl2
39         perlDL shell v1.352
40         ...
41         ReadLines, NiceSlice, MultiLines  enabled
42        ...
43         Note: AutoLoader not enabled ('use PDL::AutoLoader' recommended)
44
45         pdl>
46
47       In this example, NiceSlice was automatically enabled, but AutoLoader
48       was not.  To enable it, type "use PDL::AutoLoader".
49
50       Let's start with a two-dimensional piddle:
51
52         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
66       The "info" method gives you basic information about a piddle:
67
68         pdl> p $a->info
69         PDL: Double D [11,9]
70
71       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
75         # 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
82       Note: Notice that indices start at 0, as in Perl, C and Java (and
83       unlike MATLAB and IDL).
84
85       But with PDL, we can just write:
86
87         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
101       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
106         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
123       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
128         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
132       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
137       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
140         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
154       See PDL::Indexing and PDL::NiceSlice to learn more about specifying
155       subsets from piddles.
156
157       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
160         pdl> $piddle_4D = sequence(11,3,7,2)
161         pdl> $c = $piddle_4D - $line
162
163       Now $c is a piddle of the same dimension as $piddle_4D.
164
165         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
168       This time PDL has threaded over three higher dimensions automatically,
169       subtracting $line all the way.
170
171       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
175         pdl> $cols = sequence(9)
176         pdl> p $a->info      =>  PDL: Double D [11,9]
177         pdl> p $cols->info   =>  PDL: Double D [9]
178
179       Naturally, we can't just type "$a - $cols". The dimensions don't match:
180
181         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
185       How do we tell PDL that we want to subtract from  dimension 1 instead?
186

MANIPULATING DIMENSIONS

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
192        xchg
193        mv
194        reorder
195
196   Method: "xchg"
197       The "xchg" method "exchanges" two dimensions in a piddle:
198
199         pdl> $a = sequence(6,7,8,9)
200         pdl> $a_xchg = $a->xchg(0,3)
201
202         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
208       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
212   Method: "mv"
213       The "mv" method "moves" one dimension, in a piddle, shifting other
214       dimensions as necessary.
215
216         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
225       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
229   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
233         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
242       Notice what happened. When we wrote "reorder(3,0,2,1)" we instructed
243       PDL to:
244
245        * Put dimension 3 first.
246        * Put dimension 0 next.
247        * Put dimension 2 next.
248        * Put dimension 1 next.
249
250       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

GOTCHA: LINKING VS ASSIGNMENT

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
259         pdl> $a = sequence(4,5)
260         pdl> $a_xchg = $a->xchg(1,0)
261
262       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
265         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
283   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
287         pdl> $a = sequence(4,5)
288         pdl> $a_xchg = $a->copy->xchg(1,0)
289
290       Now $a and $a_xchg are completely separate objects:
291
292         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

PUTTING IT ALL TOGETHER

318       Now we are ready to solve the problem that motivated this whole
319       discussion:
320
321         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
327       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
331         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
357       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

EXAMPLE: CONWAY'S GAME OF LIFE

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
366       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
372       1) If an empty cell has exactly three neighbours, a living cell is
373       generated.
374
375       2) If a living cell has less than two neighbours, it dies of
376       overfeeding.
377
378       3) If a living cell has 4 or more neighbours, it dies from starvation.
379
380       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
387   Classical implementation
388       Here's a classic way of writing this program in Perl. We only use PDL
389       for addressing individual cells.
390
391         #!/usr/bin/perl -w
392         use PDL;
393         use PDL::NiceSlice;
394
395         # Make a board for the game of life.
396         my $nx = 20;
397         my $ny = 20;
398
399         # Current generation.
400         my $a = zeroes($nx, $ny);
401
402         # Next generation.
403         my $n = zeroes($nx, $ny);
404
405         # Put in a simple glider.
406         $a(1:3,1:3) .= pdl ( [1,1,1],
407                              [0,0,1],
408                              [0,1,0] );
409
410         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
416                   # 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
422                            # 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
428                           $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
434                   # 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
446           print $a;
447
448           $a = $new_a;
449         }
450
451       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
455   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
459         #!/usr/bin/perl -w
460         use PDL;
461         use PDL::NiceSlice;
462
463         my $a = zeroes(20,20);
464
465         # Put in a simple glider.
466         $a(1:3,1:3) .= pdl ( [1,1,1],
467                              [0,0,1],
468                              [0,1,0] );
469
470         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
476           # Calculate the next generation.
477           $a = ((($n == 2) + ($n == 3))* $a) + (($n==3) * !$a);
478
479           print $a;
480         }
481
482       The threaded PDL version is much faster:
483
484         Classical => 32.79 seconds.
485         Threaded  =>  0.41 seconds.
486
487   Explanation
488       How does the threaded version work?
489
490       There are many PDL functions designed to help you carry out PDL
491       threading.  In this example, the key functions are:
492
493       Method: "range"
494
495       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
499         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
517       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
521         pdl> $index = pdl [ [1,2],[2,3],[3,4],[4,5] ]
522         pdl> p $a->range( $index )
523         [13 20 27 34]
524
525       Additionally, "range" takes a second parameter which determines the
526       size of the chunk to return:
527
528         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
536       We can use this to select one or more 3x3 boxes.
537
538       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
544         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
568       Notice how the box wraps around the boundary of the piddle.
569
570       Method: "ndcoords"
571
572       The "ndcoords" method is a convenience method that returns an
573       enumerated list of coordinates suitable for use with the "range"
574       method.
575
576         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
601       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
604          (0,0)     (1,0)     (2,0)
605          (1,0)     (1,1)     (2,1)
606          (2,0)     (2,1)     (2,2)
607
608       Combining "range" and "ndcoords"
609
610       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
613         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
626       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
629         pdl> p $piddle->range( ndcoords($piddle) , 2 , "periodic" )
630
631       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
635         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
649       Here you can see more clearly that
650
651         [0 1]
652         [3 4]
653
654       Is the 2x2 box starting with the (0,0) element of $piddle.
655
656       We are not done yet. For the game of life, we want 3x3 boxes from $a:
657
658         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
679       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
685         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
700       Now we see a 3x3 box with the (0,0) element in the centre of the box.
701
702       Method: "sumover"
703
704       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
707         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
733       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
737       Counting neighbours
738
739       We are almost there!
740
741       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
744         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
755       When applied to Conway's Game of Life, this will tell us how many
756       living neighbours each cell has:
757
758         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
791       For example, this tells us that cell (0,0) has 1 living neighbour,
792       while cell (2,2) has 5 living neighbours.
793
794       Calculating the next generation
795
796       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
800       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
804              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
818            Get a list of empty cells that have exactly three neighbours:
819
820              pdl> p ($n == 3) * !$a
821
822       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
825              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
839            Get a list of living cells that have exactly 2 or 3 neighbours:
840
841              pdl> p (($n == 2) + ($n == 3)) * $a
842
843       Putting it all together, the next generation is:
844
845         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
860   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
864         #!/usr/bin/perl
865         use PDL;
866         use PDL::NiceSlice;
867         use PDL::Graphics::TriD;
868
869         my $a = zeroes(20,20);
870
871         # Put in a simple glider.
872         $a(1:3,1:3) .= pdl ( [1,1,1],
873                              [0,0,1],
874                              [0,1,0] );
875
876         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
882             # Calculate the next generation.
883             $a = ((($n == 2) + ($n == 3))* $a) + (($n==3) * !$a);
884
885             # Display.
886             nokeeptwiddling3d();
887             imagrgb [$a];
888         }
889
890       But if we really want to see something interesting, we should make a
891       few more changes:
892
893       1) Start with a random collection of 1's and 0's.
894
895       2) Make the grid larger.
896
897       3) Add a small timeout so we can see the game evolve better.
898
899       4) Use a while loop so that the program can run as long as it needs to.
900
901         #!/usr/bin/perl
902         use PDL;
903         use PDL::NiceSlice;
904         use PDL::Graphics::TriD;
905         use Time::HiRes qw(usleep);
906
907         my $a = random(100,100);
908         $a = ($a < 0.5);
909
910         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
916             # Calculate the next generation.
917             $a = ((($n == 2) + ($n == 3))* $a) + (($n==3) * !$a);
918
919             # Display.
920             nokeeptwiddling3d();
921             imagrgb [$a];
922
923             # Sleep for 0.1 seconds.
924             usleep(100000);
925         }
926

CONCLUSION: GENERAL STRATEGY

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
932       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
941       See: http://dev.perl.org/licenses/
942
943
944
945perl v5.28.1                      2018-05-05                      THREADING(1)
Impressum