

General Strategy:  Move the dimensions you want to operate on to the start of your piddle’s dimension list. Then let PDL thread over the higher dimensions. 
Okay, enough theory. Let’s do something a bit more interesting: We’ll write <B>Conway’s Game of LifeB> in PDL and see how powerful PDL can be!The <B>Game of LifeB> is a simulation run on a big two dimensional grid. Each cell in the grid can either be alive or dead (represented by 1 or 0). The next generation of cells in the grid is calculated with simple rules according to the number of living cells in it’s immediate neighbourhood:
1) If an empty cell has exactly three neighbours, a living cell is generated.
2) If a living cell has less than two neighbours, it dies of overfeeding.
3) If a living cell has 4 or more neighbours, it dies from starvation.
Only the first generation of cells is determined by the programmer. After that, the simulation runs completely according to these rules. To calculate the next generation, you need to look at each cell in the 2D field (requiring two loops), calculate the number of live cells adjacent to this cell (requiring another two loops) and then fill the next generation grid.
Here’s a classic way of writing this program in Perl. We only use PDL for addressing individual cells.
#!/usr/local/bin/perl w use PDL; use PDL::NiceSlice; # Make a board for the game of life. my $nx = 20; my $ny = 20; # Current generation. my $a = zeroes($nx, $ny); # Next generation. my $n = zeroes($nx, $ny); # Put in a simple glider. $a(1:3,1:3) .= pdl ( [1,1,1], [0,0,1], [0,1,0] ); for (my $i = 0; $i < 100; $i++) { $n = zeroes($nx, $ny); $new_a = $a>copy; for ($x = 0; $x < $nx; $x++) { for ($y = 0; $y < $ny; $y++) { # For each cell, look at the surrounding neighbours. for ($dx = 1; $dx <= 1; $dx++) { for ($dy = 1; $dy <= 1; $dy++) { $px = $x + $dx; $py = $y + $dy; # Wrap around at the edges. if ($px < 0) {$px = $nx1}; if ($py < 0) {$py = $ny1}; if ($px >= $nx) {$px = 0}; if ($py >= $ny) {$py = 0}; $n($x,$y) .= $n($x,$y) + $a($px,$py); } } # Do not count the central cell itself. $n($x,$y) = $a($x,$y); # Work out if cell lives or dies: # Dead cell lives if n = 3 # Live cell dies if n is not 2 or 3 if ($a($x,$y) == 1) { if ($n($x,$y) < 2) {$new_a($x,$y) .= 0}; if ($n($x,$y) > 3) {$new_a($x,$y) .= 0}; } else { if ($n($x,$y) == 3) {$new_a($x,$y) .= 1} } } } print $a; $a = $new_a; }If you run this, you will see a small glider crawl diagonally across the grid of zeroes. On my machine, it prints out a couple of generations per second.
And here’s the threaded version in PDL. Just four lines of PDL code, and one of those is printing out the latest generation!
#!/usr/local/bin/perl w use PDL; use PDL::NiceSlice; my $a = zeroes(20,20); # Put in a simple glider. $a(1:3,1:3) .= pdl ( [1,1,1], [0,0,1], [0,1,0] ); my $n; for (my $i = 0; $i < 100; $i++) { # Calculate the number of neighbours per cell. $n = $a>range(ndcoords($a)1,3,"periodic")>reorder(2,3,0,1); $n = $n>sumover>sumover  $a; # Calculate the next generation. $a = ((($n == 2) + ($n == 3))* $a) + (($n==3) * !$a); print $a; }The threaded PDL version is much faster:
Classical => 32.79 seconds. Threaded => 0.41 seconds.
How does the threaded version work?There are many PDL functions designed to help you carry out PDL threading. In this example, the key functions are:
Method: range
At the simplest level, the range method is a different way to select a portion of a piddle. Instead of using the $a(2,3) notation, we use another piddle.
pdl> $a = sequence(6,7) pdl> p $a [ [ 0 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] ] pdl> p $a>range( pdl [1,2] ) 13 pdl> p $a(1,2) [ [13] ]At this point, the range method looks very similar to a regular PDL slice. But the range method is more general. For example, you can select several components at once:
pdl> $index = pdl [ [1,2],[2,3],[3,4],[4,5] ] pdl> p $a>range( $index ) [13 20 27 34]Additionally, range takes a second parameter which determines the size of the chunk to return:
pdl> $size = 3 pdl> p $a>range( pdl([1,2]) , $size ) [ [13 14 15] [19 20 21] [25 26 27] ]We can use this to select one or more 3x3 boxes.
Finally, range can take a third parameter called the boundary condition. It tells PDL what to do if the size box you request goes beyond the edge of the piddle. We won’t go over all the options. We’ll just say that the option periodic means that the piddle wraps around. For example:
pdl> p $a [ [ 0 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] ] pdl> $size = 3 pdl> p $a>range( pdl([4,2]) , $size , "periodic" ) [ [16 17 12] [22 23 18] [28 29 24] ] pdl> p $a>range( pdl([5,2]) , $size , "periodic" ) [ [17 12 13] [23 18 19] [29 24 25] ]Notice how the box wraps around the boundary of the piddle.
Method: ndcoords
The ndcoords method is a convenience method that returns an enumerated list of coordinates suitable for use with the range method.
pdl> p $piddle = sequence(3,3) [ [0 1 2] [3 4 5] [6 7 8] ] pdl> p ndcoords($piddle) [ [ [0 0] [1 0] [2 0] ] [ [0 1] [1 1] [2 1] ] [ [0 2] [1 2] [2 2] ] ]This can be a little hard to read. Basically it’s saying that the coordinates for every element in $piddle is given by:
(0,0) (1,0) (2,0) (1,0) (1,1) (2,1) (2,0) (2,1) (2,2)Combining range and ndcoords
What really matters is that ndcoords is designed to work together with range, with no $size parameter, you get the same piddle back.
pdl> p $piddle [ [0 1 2] [3 4 5] [6 7 8] ] pdl> p $piddle>range( ndcoords($piddle) ) [ [0 1 2] [3 4 5] [6 7 8] ]Why would this be useful? Because now we can ask for a series of boxes for the entire piddle. For example, 2x2 boxes:
pdl> p $piddle>range( ndcoords($piddle) , 2 , "periodic" )The output of this function is difficult to read because the boxes along the last two dimension. We can make the result more readable by rearranging the dimensions:
pdl> p $piddle>range( ndcoords($piddle) , 2 , "periodic" )>reorder(2,3,0,1) [ [ [ [0 1] [3 4] ] [ [1 2] [4 5] ] ... ]Here you can see more clearly that
[0 1] [3 4]Is the 2x2 box starting with the (0,0) element of $piddle.
We are not done yet. For the game of life, we want 3x3 boxes from $a:
pdl> p $a [ [ 0 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] ] pdl> p $a>range( ndcoords($a) , 3 , "periodic" )>reorder(2,3,0,1) [ [ [ [ 0 1 2] [ 6 7 8] [12 13 14] ] ... ]We can confirm that this is the 3x3 box starting with the (0,0) element of $a. But there is one problem. We actually want the 3x3 box to be <B>centeredB> on (0,0). That’s not a problem. Just subtract 1 from all the coordinates in ndcoords($a). Remember that the periodic option takes care of making everything wrap around.
pdl> p $a>range( ndcoords($a)  1 , 3 , "periodic" )>reorder(2,3,0,1) [ [ [ [41 36 37] [ 5 0 1] [11 6 7] ] [ [36 37 38] [ 0 1 2] [ 6 7 8] ] ...Now we see a 3x3 box with the (0,0) element in the centre of the box.
Method: sumover
The sumover method adds along only the first dimension. If we apply it twice, we will be adding all the elements of each 3x3 box.
pdl> $n = $a>range(ndcoords($a)1,3,"periodic")>reorder(2,3,0,1) pdl> p $n [ [ [ [41 36 37] [ 5 0 1] [11 6 7] ] [ [36 37 38] [ 0 1 2] [ 6 7 8] ] ... pdl> p $n>sumover>sumover [ [144 135 144 153 162 153] [ 72 63 72 81 90 81] [126 117 126 135 144 135] [180 171 180 189 198 189] [234 225 234 243 252 243] [288 279 288 297 306 297] [216 207 216 225 234 225] ]Use a calculator to confirm that 144 is the sum of all the elements in the first 3x3 box and 135 is the sum of all the elements in the second 3x3 box.
Counting neighbours
We are almost there!
Adding up all the elements in a 3x3 box is not <B>quiteB> what we want. We don’t want to count the center box. Fortunately, this is an easy fix:
pdl> p $n>sumover>sumover  $a [ [144 134 142 150 158 148] [ 66 56 64 72 80 70] [114 104 112 120 128 118] [162 152 160 168 176 166] [210 200 208 216 224 214] [258 248 256 264 272 262] [180 170 178 186 194 184] ]When applied to Conway’s Game of Life, this will tell us how many living neighbours each cell has:
pdl> $a = zeroes(10,10) pdl> $a(1:3,1:3) .= pdl ( [1,1,1], ..( > [0,0,1], ..( > [0,1,0] ) pdl> p $a [ [0 0 0 0 0 0 0 0 0 0] [0 1 1 1 0 0 0 0 0 0] [0 0 0 1 0 0 0 0 0 0] [0 0 1 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] ] pdl> $n = $a>range(ndcoords($a)1,3,"periodic")>reorder(2,3,0,1) pdl> $n = $n>sumover>sumover  $a pdl> p $n [ [1 2 3 2 1 0 0 0 0 0] [1 1 3 2 2 0 0 0 0 0] [1 3 5 3 2 0 0 0 0 0] [0 1 1 2 1 0 0 0 0 0] [0 1 1 1 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] ]For example, this tells us that cell (0,0) has 1 living neighbour, while cell (2,2) has 5 living neighbours.
Calculating the next generation
At this point, the variable $n has the number of living neighbours for every cell. Now we apply the rules of the game of life to calculate the next generation.
Putting it all together, the next generation is:
If an empty cell has exactly three neighbours, a living cell is generated. Get a list of cells that have exactly three neighbours:
pdl> p ($n == 3) [ [0 0 1 0 0 0 0 0 0 0] [0 0 1 0 0 0 0 0 0 0] [0 1 0 1 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] ]Get a list of <B>emptyB> cells that have exactly three neighbours:
pdl> p ($n == 3) * !$aIf a living cell has less than 2 or more than 3 neighbours, it dies. Get a list of cells that have exactly 2 or 3 neighbours:
pdl> p (($n == 2) + ($n == 3)) [ [0 1 1 1 0 0 0 0 0 0] [0 0 1 1 1 0 0 0 0 0] [0 1 0 1 1 0 0 0 0 0] [0 0 0 1 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] ]Get a list of <B>livingB> cells that have exactly 2 or 3 neighbours:
pdl> p (($n == 2) + ($n == 3)) * $a
pdl> $a = ((($n == 2) + ($n == 3)) * $a) + (($n == 3) * !$a) pdl> p $a [ [0 0 1 0 0 0 0 0 0 0] [0 0 1 1 0 0 0 0 0 0] [0 1 0 1 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] ]
If you have PDL::Graphics::TriD installed, you can make a graphical version of the program by just changing three lines:
#!/usr/local/bin/perl use PDL; use PDL::NiceSlice; use PDL::Graphics::TriD; my $a = zeroes(20,20); # Put in a simple glider. $a(1:3,1:3) .= pdl ( [1,1,1], [0,0,1], [0,1,0] ); my $n; for (my $i = 0; $i < 100; $i++) { # Calculate the number of neighbours per cell. $n = $a>range(ndcoords($a)1,3,"periodic")>reorder(2,3,0,1); $n = $n>sumover>sumover  $a; # Calculate the next generation. $a = ((($n == 2) + ($n == 3))* $a) + (($n==3) * !$a); # Display. nokeeptwiddling3d(); imagrgb [$a]; }But if we really want to see something interesting, we should make a few more changes:
1) Start with a random collection of 1’s and 0’s.
2) Make the grid larger.
3) Add a small timeout so we can see the game evolve better.
4) Use a while loop so that the program can run as long as it needs to.
#!/usr/local/bin/perl use PDL; use PDL::NiceSlice; use PDL::Graphics::TriD; use Time::HiRes qw(usleep); my $a = random(100,100); $a = ($a < 0.5); my $n; while (1) { # Calculate the number of neighbours per cell. $n = $a>range(ndcoords($a)1,3,"periodic")>reorder(2,3,0,1); $n = $n>sumover>sumover  $a; # Calculate the next generation. $a = ((($n == 2) + ($n == 3))* $a) + (($n==3) * !$a); # Display. nokeeptwiddling3d(); imagrgb [$a]; # Sleep for 0.1 seconds. usleep(100000); }
The general strategy is: Move the dimensions you want to operate on to the start of your piddle’s dimension list. Then let PDL thread over the higher dimensions.Threading is a powerful tool that helps eliminate forloops and can make your code more concise. Hopefully this tutorial has shown why it is worth getting to grips with threading in PDL.
Copyright 2010 Matthew Kenworthy (kenworthy@strw.leidenuniv.nl) and Daniel Carrera (dcarrera@gmail.com). You can distribute and/or modify this document under the same terms as the current Perl license.
perl v5.20.3  THREADING (1)  20160403 
Visit the GSP FreeBSD Man Page Interface.
Output converted with manServer 1.07.