.\" Automatically generated by Pod::Man 4.14 (Pod::Simple 3.40)
.\"
.\" Standard preamble:
.\" ========================================================================
.de Sp \" Vertical space (when we can't use .PP)
.if t .sp .5v
.if n .sp
..
.de Vb \" Begin verbatim text
.ft CW
.nf
.ne \\$1
..
.de Ve \" End verbatim text
.ft R
.fi
..
.\" Set up some character translations and predefined strings.  \*(-- will
.\" give an unbreakable dash, \*(PI will give pi, \*(L" will give a left
.\" double quote, and \*(R" will give a right double quote.  \*(C+ will
.\" give a nicer C++.  Capital omega is used to do unbreakable dashes and
.\" therefore won't be available.  \*(C` and \*(C' expand to `' in nroff,
.\" nothing in troff, for use with C<>.
.tr \(*W-
.ds C+ C\v'-.1v'\h'-1p'\s-2+\h'-1p'+\s0\v'.1v'\h'-1p'
.ie n \{\
.    ds -- \(*W-
.    ds PI pi
.    if (\n(.H=4u)&(1m=24u) .ds -- \(*W\h'-12u'\(*W\h'-12u'-\" diablo 10 pitch
.    if (\n(.H=4u)&(1m=20u) .ds -- \(*W\h'-12u'\(*W\h'-8u'-\"  diablo 12 pitch
.    ds L" ""
.    ds R" ""
.    ds C` ""
.    ds C' ""
'br\}
.el\{\
.    ds -- \|\(em\|
.    ds PI \(*p
.    ds L" ``
.    ds R" ''
.    ds C`
.    ds C'
'br\}
.\"
.\" Escape single quotes in literal strings from groff's Unicode transform.
.ie \n(.g .ds Aq \(aq
.el       .ds Aq '
.\"
.\" If the F register is >0, we'll generate index entries on stderr for
.\" titles (.TH), headers (.SH), subsections (.SS), items (.Ip), and index
.\" entries marked with X<> in POD.  Of course, you'll have to process the
.\" output yourself in some meaningful fashion.
.\"
.\" Avoid warning from groff about undefined register 'F'.
.de IX
..
.nr rF 0
.if \n(.g .if rF .nr rF 1
.if (\n(rF:(\n(.g==0)) \{\
.    if \nF \{\
.        de IX
.        tm Index:\\$1\t\\n%\t"\\$2"
..
.        if !\nF==2 \{\
.            nr % 0
.            nr F 2
.        \}
.    \}
.\}
.rr rF
.\" ========================================================================
.\"
.IX Title "THREADING 1p"
.TH THREADING 1p "2022-01-30" "perl v5.32.1" "User Contributed Perl Documentation"
.\" For nroff, turn off justification.  Always turn off hyphenation; it makes
.\" way too many mistakes in technical documents.
.if n .ad l
.nh
.SH "NAME"
PDL::Threading \- Tutorial for PDL's Threading feature
.SH "INTRODUCTION"
.IX Header "INTRODUCTION"
One of the most powerful features of \s-1PDL\s0 is \fBthreading\fR, which can
produce very compact and very fast \s-1PDL\s0 code by avoiding multiple nested
for loops that C and \s-1BASIC\s0 users may be familiar with. The trouble is
that it can take some getting used to, and new users may not appreciate
the benefits of threading.
.PP
Other vector based languages, such as \s-1MATLAB,\s0 use a subset of threading
techniques, but \s-1PDL\s0 shines by completely generalizing them for all sorts
of vector-based applications.
.SH "TERMINOLOGY: NDARRAY"
.IX Header "TERMINOLOGY: NDARRAY"
\&\s-1MATLAB\s0 typically refers to vectors, matrices, and arrays. Perl already
has arrays, and the terms \*(L"vector\*(R" and \*(L"matrix\*(R" typically refer to one\-
and two-dimensional collections of data. Having no good term to describe
their object, \s-1PDL\s0 developers coined the term "\fIndarray\fR" to give a name to
their data type.
.PP
An \fIndarray\fR consists of a series of numbers organized as an N\-dimensional
data set. ndarrays provide efficient storage and fast computation of large
N\-dimensional matrices. They are highly optimized for numerical work.
.SH "THINKING IN TERMS OF THREADING"
.IX Header "THINKING IN TERMS OF THREADING"
If you have used \s-1PDL\s0 for a little while already, you may have been using
threading without realising it. Start the \s-1PDL\s0 shell (type \f(CW\*(C`perldl\*(C'\fR or
\&\f(CW\*(C`pdl2\*(C'\fR on a terminal). Most examples in this tutorial use the \s-1PDL\s0 shell.
Make sure that PDL::NiceSlice and PDL::AutoLoader
are enabled. For example:
.PP
.Vb 6
\&  % pdl2
\&  perlDL shell v1.352
\&  ...
\&  ReadLines, NiceSlice, MultiLines  enabled
\& ...
\&  Note: AutoLoader not enabled (\*(Aquse PDL::AutoLoader\*(Aq recommended)
\&
\&  pdl>
.Ve
.PP
In this example, NiceSlice was automatically enabled, but AutoLoader was not.
To enable it, type \f(CW\*(C`use PDL::AutoLoader\*(C'\fR.
.PP
Let's start with a two-dimensional \fIndarray\fR:
.PP
.Vb 10
\&  pdl> $x = sequence(11,9)
\&  pdl> p $x
\&  [
\&    [ 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 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]
\&  ]
.Ve
.PP
The \f(CW\*(C`info\*(C'\fR method gives you basic information about an \fIndarray\fR:
.PP
.Vb 2
\&  pdl> p $x\->info
\&  PDL: Double D [11,9]
.Ve
.PP
This tells us that \f(CW$x\fR is an 11 x 9 \fIndarray\fR composed of double
precision numbers. If we wanted to add 3 to all elements in an \f(CW\*(C`n x m\*(C'\fR
ndarray, a traditional language would use two nested for-loops:
.PP
.Vb 6
\&  # Pseudo\-code. Traditional way to add 3 to an array.
\&  for (i=0; i < n; i++) {
\&      for (j=0; j < m; j++) {
\&          a(i,j) = a(i,j) + 3
\&      }
\&  }
.Ve
.PP
\&\fBNote\fR: Notice that indices start at 0, as in Perl, C and Java (and unlike
\&\s-1MATLAB\s0 and \s-1IDL\s0).
.PP
But with \s-1PDL,\s0 we can just write:
.PP
.Vb 10
\&  pdl> $y = $x + 3
\&  pdl> p $y
\&  [
\&    [  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]
\&  ]
.Ve
.PP
This is the simplest example of threading, and it is something that
all numerical software tools do. The \f(CW\*(C`+ 3\*(C'\fR operation was automatically
applied along two dimensions. Now suppose you want to to subtract a
line from every row in \f(CW$x\fR:
.PP
.Vb 10
\&  pdl> $line = sequence(11)
\&  pdl> p $line
\&  [0 1 2 3 4 5 6 7 8 9 10]
\&  pdl> $c = $x \- $line
\&  pdl> p $c
\&  [
\&   [ 0  0  0  0  0  0  0  0  0  0  0]
\&   [11 11 11 11 11 11 11 11 11 11 11]
\&   [22 22 22 22 22 22 22 22 22 22 22]
\&   [33 33 33 33 33 33 33 33 33 33 33]
\&   [44 44 44 44 44 44 44 44 44 44 44]
\&   [55 55 55 55 55 55 55 55 55 55 55]
\&   [66 66 66 66 66 66 66 66 66 66 66]
\&   [77 77 77 77 77 77 77 77 77 77 77]
\&   [88 88 88 88 88 88 88 88 88 88 88]
\&  ]
.Ve
.PP
Two things to note here: First, the value of \f(CW$x\fR is still the same. Try
\&\f(CW\*(C`p $x\*(C'\fR to check. Second, \s-1PDL\s0 automatically subtracted \f(CW$line\fR from each
row in \f(CW$x\fR. Why did it do that? Let's look at the dimensions of \f(CW$x\fR,
\&\f(CW$line\fR and \f(CW$c\fR:
.PP
.Vb 3
\&  pdl> p $line\->info  =>  PDL: Double D [11]
\&  pdl> p $x\->info     =>  PDL: Double D [11,9]
\&  pdl> p $c\->info     =>  PDL: Double D [11,9]
.Ve
.PP
So, both \f(CW$x\fR and \f(CW$line\fR have the same number of elements in the 0th
dimension! What \s-1PDL\s0 then did was thread over the higher dimensions in \f(CW$x\fR
and repeated the same operation 9 times to all the rows on \f(CW$x\fR. This is
\&\s-1PDL\s0 threading in action.
.PP
What if you want to subtract \f(CW$line\fR from the first line in \f(CW$x\fR only?
You can do that by specifying the line explicitly:
.PP
.Vb 10
\&  pdl> $x(:,0) \-= $line
\&  pdl> p $x
\&  [
\&   [ 0  0  0  0  0  0  0  0  0  0  0]
\&   [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]
\&  ]
.Ve
.PP
See PDL::Indexing and PDL::NiceSlice to
learn more about specifying subsets from ndarrays.
.PP
The true power of threading comes when you realise that the ndarray can
have any number of dimensions! Let's make a 4 dimensional ndarray:
.PP
.Vb 2
\&  pdl> $ndarray_4D = sequence(11,3,7,2)
\&  pdl> $c = $ndarray_4D \- $line
.Ve
.PP
Now \f(CW$c\fR is an ndarray of the same dimension as \f(CW$ndarray_4D\fR.
.PP
.Vb 2
\&  pdl> p $ndarray_4D\->info  =>  PDL: Double D [11,3,7,2]
\&  pdl> p $c\->info          =>  PDL: Double D [11,3,7,2]
.Ve
.PP
This time \s-1PDL\s0 has threaded over three higher dimensions automatically,
subtracting \f(CW$line\fR all the way.
.PP
But, maybe you don't want to subtract from the rows (dimension 0), but from
the columns (dimension 1). How do I subtract a column of numbers from each
column in \f(CW$x\fR?
.PP
.Vb 3
\&  pdl> $cols = sequence(9)
\&  pdl> p $x\->info      =>  PDL: Double D [11,9]
\&  pdl> p $cols\->info   =>  PDL: Double D [9]
.Ve
.PP
Naturally, we can't just type \f(CW\*(C`$x \- $cols\*(C'\fR. The dimensions don't match:
.PP
.Vb 3
\&  pdl> p $x \- $cols
\&  PDL: PDL::Ops::minus(a,b,c): Parameter \*(Aqb\*(Aq
\&  PDL: Mismatched implicit thread dimension 0: should be 11, is 9
.Ve
.PP
How do we tell \s-1PDL\s0 that we want to subtract from  dimension 1 instead?
.SH "MANIPULATING DIMENSIONS"
.IX Header "MANIPULATING DIMENSIONS"
There are many \s-1PDL\s0 functions that let you rearrange the dimensions of \s-1PDL\s0
arrays. They are mostly covered in PDL::Slices. The three
most common ones are:
.PP
.Vb 3
\& xchg
\& mv
\& reorder
.Ve
.ie n .SS "Method: ""xchg"""
.el .SS "Method: \f(CWxchg\fP"
.IX Subsection "Method: xchg"
The \f(CW\*(C`xchg\*(C'\fR method "\fBexchanges\fR" two dimensions in an ndarray:
.PP
.Vb 2
\&  pdl> $x = sequence(6,7,8,9)
\&  pdl> $x_xchg = $x\->xchg(0,3)
\&  
\&  pdl> p $x\->info       =>  PDL: Double D [6,7,8,9]
\&  pdl> p $x_xchg\->info  =>  PDL: Double D [9,7,8,6]
\&                                           |     |
\&                                           V     V
\&                                       (dim 0) (dim 3)
.Ve
.PP
Notice that dimensions 0 and 3 were exchanged without affecting the other
dimensions. Notice also that \f(CW\*(C`xchg\*(C'\fR does not alter \f(CW$x\fR. The original
variable \f(CW$x\fR remains untouched.
.ie n .SS "Method: ""mv"""
.el .SS "Method: \f(CWmv\fP"
.IX Subsection "Method: mv"
The \f(CW\*(C`mv\*(C'\fR method "\fBmoves\fR" one dimension, in an ndarray, shifting other
dimensions as necessary.
.PP
.Vb 8
\&  pdl> $x = sequence(6,7,8,9)         (dim 0)
\&  pdl> $x_mv = $x\->mv(0,3)               |
\&  pdl>                                   V _\|_\|_\|_\|_
\&  pdl> p $x\->info     =>  PDL: Double D [6,7,8,9]
\&  pdl> p $x_mv\->info  =>  PDL: Double D [7,8,9,6]
\&                                          \-\-\-\-\- |
\&                                                V
\&                                              (dim 3)
.Ve
.PP
Notice that when dimension 0 was moved to position 3, all the other dimensions
had to be shifted as well. Notice also that \f(CW\*(C`mv\*(C'\fR does not alter \f(CW$x\fR. The
original variable \f(CW$x\fR remains untouched.
.ie n .SS "Method: ""reorder"""
.el .SS "Method: \f(CWreorder\fP"
.IX Subsection "Method: reorder"
The \f(CW\*(C`reorder\*(C'\fR method is a generalization of the \f(CW\*(C`xchg\*(C'\fR and \f(CW\*(C`mv\*(C'\fR methods.
It "\fBreorders\fR" the dimensions in any way you specify:
.PP
.Vb 8
\&  pdl> $x = sequence(6,7,8,9)
\&  pdl> $x_reorder = $x\->reorder(3,0,2,1)
\&  pdl>
\&  pdl> p $x\->info          =>  PDL: Double D [6,7,8,9]
\&  pdl> p $x_reorder\->info  =>  PDL: Double D [9,6,8,7]
\&                                              | | | |
\&                                              V V v V
\&                                 dimensions:  0 1 2 3
.Ve
.PP
Notice what happened. When we wrote \f(CW\*(C`reorder(3,0,2,1)\*(C'\fR we instructed \s-1PDL\s0 to:
.PP
.Vb 4
\& * Put dimension 3 first.
\& * Put dimension 0 next.
\& * Put dimension 2 next.
\& * Put dimension 1 next.
.Ve
.PP
When you use the \f(CW\*(C`reorder\*(C'\fR method, all the dimensions are shuffled. Notice that
\&\f(CW\*(C`reorder\*(C'\fR does not alter \f(CW$x\fR. The original variable \f(CW$x\fR remains untouched.
.SH "GOTCHA: LINKING VS ASSIGNMENT"
.IX Header "GOTCHA: LINKING VS ASSIGNMENT"
.SS "Linking"
.IX Subsection "Linking"
By default, ndarrays are \fBlinked together\fR so that changes on one will go
back and affect the original \fBas well\fR.
.PP
.Vb 2
\&  pdl> $x = sequence(4,5)
\&  pdl> $x_xchg = $x\->xchg(1,0)
.Ve
.PP
Here, \f(CW$x_xchg\fR \fBis not a separate object\fR. It is merely a different way
of looking at \f(CW$x\fR. Any change in \f(CW$x_xchg\fR will appear in \f(CW$x\fR as well.
.PP
.Vb 10
\&  pdl> p $x
\&  [
\&   [ 0  1  2  3]
\&   [ 4  5  6  7]
\&   [ 8  9 10 11]
\&   [12 13 14 15]
\&   [16 17 18 19]
\&  ]
\&  pdl> $x_xchg += 3
\&  pdl> p $x
\&  [
\&   [ 3  4  5  6]
\&   [ 7  8  9 10]
\&   [11 12 13 14]
\&   [15 16 17 18]
\&   [19 20 21 22]
\&  ]
.Ve
.SS "Assignment"
.IX Subsection "Assignment"
Some times, linking is not the behaviour you want. If you want to make the
ndarrays independent, use the \f(CW\*(C`copy\*(C'\fR method:
.PP
.Vb 2
\&  pdl> $x = sequence(4,5)
\&  pdl> $x_xchg = $x\->copy\->xchg(1,0)
.Ve
.PP
Now \f(CW$x\fR and \f(CW$x_xchg\fR are completely separate objects:
.PP
.Vb 10
\&  pdl> p $x
\&  [
\&   [ 0  1  2  3]
\&   [ 4  5  6  7]
\&   [ 8  9 10 11]
\&   [12 13 14 15]
\&   [16 17 18 19]
\&  ]
\&  pdl> $x_xchg += 3
\&  pdl> p $x
\&  [
\&   [ 0  1  2  3]
\&   [ 4  5  6  7]
\&   [ 8  9 10 11]
\&   [12 13 14 15]
\&   [16 17 18 19]
\&  ]
\&  pdl> $x_xchg
\&  [
\&   [ 3  7 11 15 19]
\&   [ 4  8 12 16 20]
\&   [ 5  9 13 17 21]
\&   [ 6 10 14 18 22]
\&  ]
.Ve
.SH "PUTTING IT ALL TOGETHER"
.IX Header "PUTTING IT ALL TOGETHER"
Now we are ready to solve the problem that motivated this whole discussion:
.PP
.Vb 5
\&  pdl> $x = sequence(11,9)
\&  pdl> $cols = sequence(9)
\&  pdl>
\&  pdl> p $x\->info     =>  PDL: Double D [11,9]
\&  pdl> p $cols\->info  =>  PDL: Double D [9]
.Ve
.PP
How do we tell \s-1PDL\s0 to subtract \f(CW$cols\fR along dimension 1 instead of dimension 0?
The simplest way is to use the \f(CW\*(C`xchg\*(C'\fR method and rely on \s-1PDL\s0 linking:
.PP
.Vb 10
\&  pdl> p $x
\&  [
\&   [ 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 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]
\&  ]
\&  pdl> $x\->xchg(1,0) \-= $cols
\&  pdl> p $x
\&  [
\&   [ 0  1  2  3  4  5  6  7  8  9 10]
\&   [10 11 12 13 14 15 16 17 18 19 20]
\&   [20 21 22 23 24 25 26 27 28 29 30]
\&   [30 31 32 33 34 35 36 37 38 39 40]
\&   [40 41 42 43 44 45 46 47 48 49 50]
\&   [50 51 52 53 54 55 56 57 58 59 60]
\&   [60 61 62 63 64 65 66 67 68 69 70]
\&   [70 71 72 73 74 75 76 77 78 79 80]
\&   [80 81 82 83 84 85 86 87 88 89 90]
\&  ]
.Ve
.IP "General Strategy:" 5
.IX Item "General Strategy:"
Move the dimensions you want to operate on to the start of your ndarray's
dimension list. Then let \s-1PDL\s0 thread over the higher dimensions.
.SH "EXAMPLE: CONWAY'S GAME OF LIFE"
.IX Header "EXAMPLE: CONWAY'S GAME OF LIFE"
Okay, enough theory. Let's do something a bit more interesting: We'll write
\&\fBConway's Game of Life\fR in \s-1PDL\s0 and see how powerful \s-1PDL\s0 can be!
.PP
The \fBGame of Life\fR 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:
.PP
1) If an empty cell has exactly three neighbours, a living cell is generated.
.PP
2) If a living cell has less than two neighbours, it dies of overfeeding.
.PP
3) If a living cell has 4 or more neighbours, it dies from starvation.
.PP
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.
.SS "Classical implementation"
.IX Subsection "Classical implementation"
Here's a classic way of writing this program in Perl. We only use \s-1PDL\s0 for
addressing individual cells.
.PP
.Vb 3
\&  #!/usr/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 $a1 = zeroes($nx, $ny);
\&  
\&  # Next generation.
\&  my $n = zeroes($nx, $ny);
\&  
\&  # Put in a simple glider.
\&  $a1(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 = $a1\->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 = $nx\-1};
\&                     if ($py < 0) {$py = $ny\-1};
\&                     if ($px >= $nx) {$px = 0};
\&                     if ($py >= $ny) {$py = 0};
\&  
\&                    $n($x,$y)  .= $n($x,$y) + $a1($px,$py);
\&                }
\&            }
\&            # Do not count the central cell itself.
\&            $n($x,$y) \-= $a1($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 ($a1($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 $a1;
\&  
\&    $a1 = $new_a;
\&  }
.Ve
.PP
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.
.SS "Threaded \s-1PDL\s0 implementation"
.IX Subsection "Threaded PDL implementation"
And here's the threaded version in \s-1PDL.\s0 Just four lines of \s-1PDL\s0 code, and
one of those is printing out the latest generation!
.PP
.Vb 3
\&  #!/usr/bin/perl \-w
\&  use PDL;
\&  use PDL::NiceSlice;
\&  
\&  my $x = zeroes(20,20);
\&  
\&  # Put in a simple glider.
\&  $x(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 = $x\->range(ndcoords($x)\-1,3,"periodic")\->reorder(2,3,0,1);
\&    $n = $n\->sumover\->sumover \- $x;
\&    
\&    # Calculate the next generation.
\&    $x = ((($n == 2) + ($n == 3))* $x) + (($n==3) * !$x);
\&    
\&    print $x;
\&  }
.Ve
.PP
The threaded \s-1PDL\s0 version is much faster:
.PP
.Vb 2
\&  Classical => 32.79 seconds.
\&  Threaded  =>  0.41 seconds.
.Ve
.SS "Explanation"
.IX Subsection "Explanation"
How does the threaded version work?
.PP
There are many \s-1PDL\s0 functions designed to help you carry out \s-1PDL\s0 threading.
In this example, the key functions are:
.PP
\fIMethod: \f(CI\*(C`range\*(C'\fI\fR
.IX Subsection "Method: range"
.PP
At the simplest level, the \f(CW\*(C`range\*(C'\fR method is a different way to select a
portion of an ndarray. Instead of using the \f(CW\*(C`$x(2,3)\*(C'\fR notation, we use
another ndarray.
.PP
.Vb 10
\&  pdl> $x = sequence(6,7)
\&  pdl> p $x
\&  [
\&   [ 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 $x\->range( pdl [1,2] )
\&  13
\&  pdl> p $x(1,2)
\&  [
\&   [13]
\&  ]
.Ve
.PP
At this point, the \f(CW\*(C`range\*(C'\fR method looks very similar to a regular \s-1PDL\s0 slice.
But the \f(CW\*(C`range\*(C'\fR method is more general. For example, you can select several
components at once:
.PP
.Vb 3
\&  pdl> $index = pdl [ [1,2],[2,3],[3,4],[4,5] ]
\&  pdl> p $x\->range( $index )
\&  [13 20 27 34]
.Ve
.PP
Additionally, \f(CW\*(C`range\*(C'\fR takes a second parameter which determines the size
of the chunk to return:
.PP
.Vb 7
\&  pdl> $size = 3
\&  pdl> p $x\->range( pdl([1,2]) , $size )
\&  [
\&   [13 14 15]
\&   [19 20 21]
\&   [25 26 27]
\&  ]
.Ve
.PP
We can use this to select one or more 3x3 boxes.
.PP
Finally, \f(CW\*(C`range\*(C'\fR can take a third parameter called the \*(L"boundary\*(R" condition.
It tells \s-1PDL\s0 what to do if the size box you request goes beyond the edge of
the ndarray. We won't go over all the options. We'll just say that the option
\&\f(CW\*(C`periodic\*(C'\fR means that the ndarray \*(L"wraps around\*(R". For example:
.PP
.Vb 10
\&  pdl> p $x
\&  [
\&   [ 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 $x\->range( pdl([4,2]) , $size , "periodic" )
\&  [
\&   [16 17 12]
\&   [22 23 18]
\&   [28 29 24]
\&  ]
\&  pdl> p $x\->range( pdl([5,2]) , $size , "periodic" )
\&  [
\&   [17 12 13]
\&   [23 18 19]
\&   [29 24 25]
\&  ]
.Ve
.PP
Notice how the box wraps around the boundary of the ndarray.
.PP
\fIMethod: \f(CI\*(C`ndcoords\*(C'\fI\fR
.IX Subsection "Method: ndcoords"
.PP
The \f(CW\*(C`ndcoords\*(C'\fR method is a convenience method that returns an enumerated
list of coordinates suitable for use with the \f(CW\*(C`range\*(C'\fR method.
.PP
.Vb 10
\&  pdl> p $ndarray = sequence(3,3)
\&  [
\&   [0 1 2]
\&   [3 4 5]
\&   [6 7 8]
\&  ]
\&  pdl> p ndcoords($ndarray)
\&  [
\&   [
\&    [0 0]
\&    [1 0]
\&    [2 0]
\&   ]
\&   [
\&    [0 1]
\&    [1 1]
\&    [2 1]
\&   ]
\&   [
\&    [0 2]
\&    [1 2]
\&    [2 2]
\&   ]
\&  ]
.Ve
.PP
This can be a little hard to read. Basically it's saying that the coordinates
for every element in \f(CW$ndarray\fR is given by:
.PP
.Vb 3
\&   (0,0)     (1,0)     (2,0)
\&   (1,0)     (1,1)     (2,1)
\&   (2,0)     (2,1)     (2,2)
.Ve
.PP
\fICombining \f(CI\*(C`range\*(C'\fI and \f(CI\*(C`ndcoords\*(C'\fI\fR
.IX Subsection "Combining range and ndcoords"
.PP
What really matters is that \f(CW\*(C`ndcoords\*(C'\fR is designed to work together with
\&\f(CW\*(C`range\*(C'\fR, with no \f(CW$size\fR parameter, you get the same ndarray back.
.PP
.Vb 12
\&  pdl> p $ndarray
\&  [
\&   [0 1 2]
\&   [3 4 5]
\&   [6 7 8]
\&  ]
\&  pdl> p $ndarray\->range( ndcoords($ndarray) )
\&  [
\&   [0 1 2]
\&   [3 4 5]
\&   [6 7 8]
\&  ]
.Ve
.PP
Why would this be useful? Because now we can ask for a series of \*(L"boxes\*(R" for
the entire ndarray. For example, 2x2 boxes:
.PP
.Vb 1
\&  pdl> p $ndarray\->range( ndcoords($ndarray) , 2 , "periodic" )
.Ve
.PP
The output of this function is difficult to read because the \*(L"boxes\*(R" along
the last two dimension. We can make the result more readable by rearranging
the dimensions:
.PP
.Vb 10
\&  pdl> p $ndarray\->range( ndcoords($ndarray) , 2 , "periodic" )\->reorder(2,3,0,1)
\&  [
\&   [
\&    [
\&     [0 1]
\&     [3 4]
\&    ]
\&    [
\&     [1 2]
\&     [4 5]
\&    ]
\&    ...
\&  ]
.Ve
.PP
Here you can see more clearly that
.PP
.Vb 2
\&  [0 1]
\&  [3 4]
.Ve
.PP
Is the 2x2 box starting with the (0,0) element of \f(CW$ndarray\fR.
.PP
We are not done yet. For the game of life, we want 3x3 boxes from \f(CW$x\fR:
.PP
.Vb 10
\&  pdl> p $x
\&  [
\&   [ 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 $x\->range( ndcoords($x) , 3 , "periodic" )\->reorder(2,3,0,1)
\&  [
\&   [
\&    [
\&     [ 0  1  2]
\&     [ 6  7  8]
\&     [12 13 14]
\&    ]
\&    ...
\&  ]
.Ve
.PP
We can confirm that this is the 3x3 box starting with the (0,0) element of \f(CW$x\fR.
But there is one problem. We actually want the 3x3 box to be \fBcentered\fR on
(0,0). That's not a problem. Just subtract 1 from all the coordinates in
\&\f(CW\*(C`ndcoords($x)\*(C'\fR. Remember that the \*(L"periodic\*(R" option takes care of making
everything wrap around.
.PP
.Vb 10
\&  pdl> p $x\->range( ndcoords($x) \- 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]
\&    ]
\&    ...
.Ve
.PP
Now we see a 3x3 box with the (0,0) element in the centre of the box.
.PP
\fIMethod: \f(CI\*(C`sumover\*(C'\fI\fR
.IX Subsection "Method: sumover"
.PP
The \f(CW\*(C`sumover\*(C'\fR method adds along only the first dimension. If we apply it
twice, we will be adding all the elements of each 3x3 box.
.PP
.Vb 10
\&  pdl> $n = $x\->range(ndcoords($x)\-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]
\&  ]
.Ve
.PP
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.
.PP
\fICounting neighbours\fR
.IX Subsection "Counting neighbours"
.PP
We are almost there!
.PP
Adding up all the elements in a 3x3 box is not \fBquite\fR what we want. We
don't want to count the center box. Fortunately, this is an easy fix:
.PP
.Vb 10
\&  pdl> p $n\->sumover\->sumover \- $x
\&  [
\&   [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]
\&  ]
.Ve
.PP
When applied to Conway's Game of Life, this will tell us how many living
neighbours each cell has:
.PP
.Vb 10
\&  pdl> $x = zeroes(10,10)
\&  pdl> $x(1:3,1:3) .= pdl ( [1,1,1],
\&  ..(    >                  [0,0,1],
\&  ..(    >                  [0,1,0] )
\&  pdl> p $x
\&  [
\&   [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 = $x\->range(ndcoords($x)\-1,3,"periodic")\->reorder(2,3,0,1)
\&  pdl> $n = $n\->sumover\->sumover \- $x
\&  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]
\&  ]
.Ve
.PP
For example, this tells us that cell (0,0) has 1 living neighbour, while
cell (2,2) has 5 living neighbours.
.PP
\fICalculating the next generation\fR
.IX Subsection "Calculating the next generation"
.PP
At this point, the variable \f(CW$n\fR has the number of living neighbours for
every cell. Now we apply the rules of the game of life to calculate the next
generation.
.IP "If an empty cell has exactly three neighbours, a living cell is generated." 5
.IX Item "If an empty cell has exactly three neighbours, a living cell is generated."
Get a list of cells that have exactly three neighbours:
.Sp
.Vb 10
\&  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]
\&  ]
.Ve
.Sp
Get a list of \fBempty\fR cells that have exactly three neighbours:
.Sp
.Vb 1
\&  pdl> p ($n == 3) * !$x
.Ve
.IP "If a living cell has less than 2 or more than 3 neighbours, it dies." 5
.IX Item "If 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:
.Sp
.Vb 10
\&  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]
\&  ]
.Ve
.Sp
Get a list of \fBliving\fR cells that have exactly 2 or 3 neighbours:
.Sp
.Vb 1
\&  pdl> p (($n == 2) + ($n == 3)) * $x
.Ve
.PP
Putting it all together, the next generation is:
.PP
.Vb 10
\&  pdl> $x = ((($n == 2) + ($n == 3)) * $x) + (($n == 3) * !$x)
\&  pdl> p $x
\&  [
\&   [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]
\&  ]
.Ve
.SS "Bonus feature: Graphics!"
.IX Subsection "Bonus feature: Graphics!"
If you have PDL::Graphics::TriD installed, you can
make a graphical version of the program by just changing three lines:
.PP
.Vb 4
\&  #!/usr/bin/perl
\&  use PDL;
\&  use PDL::NiceSlice;
\&  use PDL::Graphics::TriD;
\&  
\&  my $x = zeroes(20,20);
\&  
\&  # Put in a simple glider.
\&  $x(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 = $x\->range(ndcoords($x)\-1,3,"periodic")\->reorder(2,3,0,1);
\&      $n = $n\->sumover\->sumover \- $x;
\&      
\&      # Calculate the next generation.
\&      $x = ((($n == 2) + ($n == 3))* $x) + (($n==3) * !$x);
\&    
\&      # Display.
\&      nokeeptwiddling3d();
\&      imagrgb [$x];
\&  }
.Ve
.PP
But if we really want to see something interesting, we should make a few more
changes:
.PP
1) Start with a random collection of 1's and 0's.
.PP
2) Make the grid larger.
.PP
3) Add a small timeout so we can see the game evolve better.
.PP
4) Use a while loop so that the program can run as long as it needs to.
.PP
.Vb 5
\&  #!/usr/bin/perl
\&  use PDL;
\&  use PDL::NiceSlice;
\&  use PDL::Graphics::TriD;
\&  use Time::HiRes qw(usleep);
\&  
\&  my $x = random(100,100);
\&  $x = ($x < 0.5);
\&  
\&  my $n;
\&  while (1) {
\&      # Calculate the number of neighbours per cell.
\&      $n = $x\->range(ndcoords($x)\-1,3,"periodic")\->reorder(2,3,0,1);
\&      $n = $n\->sumover\->sumover \- $x;
\&    
\&      # Calculate the next generation.
\&      $x = ((($n == 2) + ($n == 3))* $x) + (($n==3) * !$x);
\&      
\&      # Display.
\&      nokeeptwiddling3d();
\&      imagrgb [$x];
\&      
\&      # Sleep for 0.1 seconds.
\&      usleep(100000);
\&  }
.Ve
.SH "CONCLUSION: GENERAL STRATEGY"
.IX Header "CONCLUSION: GENERAL STRATEGY"
The general strategy is: \fIMove the dimensions you want to operate on to the
start of your ndarray's dimension list. Then let \s-1PDL\s0 thread over the higher
dimensions.\fR
.PP
Threading is a powerful tool that helps eliminate for-loops and can make your
code more concise. Hopefully this tutorial has shown why it is worth getting
to grips with threading in \s-1PDL.\s0
.SH "COPYRIGHT"
.IX Header "COPYRIGHT"
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.
.PP
See: http://dev.perl.org/licenses/