PDL::Indexing(3pm)

📦

version 2.106 · 5 versions

💖 PDL::Indexing — bestie's guide to indexing and slicing ndarrays 💖

(this is PDL::Indexing(3pm), from the User Contributed Perl Documentation, and it's literally your new fave tutorial)

✨ NAME ✨

PDL::Indexing — the slay intro to indexing and slicing ndarrays, no cap.

🌸 OVERVIEW 🌸

okay so this man page is your first-ever tutorial on the indexing and broadcasting glow-up of PDL. Think of it as the group chat where you finally figure out what everyone's been talking about.

Like every vectorized language out there, PDL automates the looping over multi-dimensional data structures ("ndarrays") using a variant of mathematical vector notation. That automatic looping is called "broadcasting" — same energy as NumPy and Julia, fr fr. And PDL also automatically runs broadcasting computation in parallel — go peek at PDL::ParallelCPU for that.

A huge chunk of PDL's flexibility and power lives in the indexing and broadcasting features of the Perl extension. Indexing = you get to access an ndarray's data in a super flexible way, periodt. Broadcasting = efficient vectorization of simple operations, which is iconic.

Here's the tea: the values of an ndarray are stored compactly as typed values in a single block of memory, not — like a normal Perl list-of-lists — as individual Perl scalars. Different vibe entirely.

In the sections coming up, lots of "methods" get called out — these are Perl operators that apply to ndarrays. From the perldl shell, you can find out more about each method by typing ? followed by the method name. iykyk.

🧠 Memory layout and indexing order

PDL follows a C-like memory layout; the fastest-moving index is generally thought of as a column-index into a row. Blame PDL's roots in image-processing for this one 💅. When a 2D ndarray is printed out, the 0-th dim is printed as rows:

# perldl with "do_print" set to true
pdl> sequence(3,2)
[
 [0 1 2]
 [3 4 5]
]
# in memory, the values are: 0 1 2 3 4 5
pdl> sequence(3,2)->at(0,1)
3

So this uses "column-major" indexing, though the layout is "row-major". Chaos? A little. It's different from the mathematical convention for matrices, where 0,1 would be (row,column), and give a value of 1. And that convention also ignores the zero-based indexing in C-like languages, versus the 1-based indexing in Fortran and maths. Not me explaining three indexing systems in one paragraph 😭

For interoperability between PDL and Fortran (and also OpenGL), an ndarray would need to be transposed with the transpose operation. Alternatively, functions in either of those often have a transpose flag for smooth interoperation. 🦋

📐 Dimension lists

A ndarray (PDL variable), in general, is an N-dimensional array where N can be:

0 (for a scalar — yes, a whole ndarray for one lil number),
1 (e.g. for a sound sample),
or higher values for images and more complex structures.

Every dimension of the ndarray has a positive integer size. The perl interpreter treats each ndarray as a special type of Perl scalar (a blessed Perl object, actually — but you don't have to know that to use them) that can be used anywhere you can put a normal scalar. So sneaky, so convenient. ✨

You can access the dimensions of an ndarray as a Perl list and otherwise determine the size of an ndarray with several methods. The important ones are:

nelem — the total number of elements in an ndarray
ndims — returns the number of dimensions in an ndarray
dims — returns the dimension list of an ndarray as a Perl list
dim — returns the size of a particular dimension of an ndarray

🔗 Indexing and Dataflow

PDL maintains a notion of "dataflow" between an ndarray and indexed subfields of that ndarray. When you produce an indexed subfield or single element of a parent ndarray, the child and parent remain attached until you manually disconnect them. That's the whole magic trick, bestie.

This lets you represent the same data different ways within your code — for example, you can consider an RGB image simultaneously as a collection of (R,G,B) values in a 3 x 1000 x 1000 image, and as three separate 1000 x 1000 color planes stored in different variables. Modifying any of the variables changes the underlying memory, and the changes are reflected in all representations of the data. One memory block, many main character moments. 👑

There are two important methods that let you control dataflow connections between a child and parent ndarray:

copy — forces an explicit copy of an ndarray
sever — breaks the dataflow connection between an ndarray and its parents (if any)

🌈 Broadcasting and Dimension Order

Most PDL operations act on the first few dimensions of their ndarray arguments. For example, sumover sums all elements along the first dimension in the list (dimension 0). If you feed in a three-dimensional ndarray, then the first dimension is considered the "active" dimension and the later dimensions are "broadcast" dimensions because they are simply looped over. Looping, but make it chic.

There are several ways to transpose or re-order the dimension list of an ndarray. Those techniques are very fast since they don't touch the underlying data, only change the way that PDL accesses the data. Iconic behaviour, honestly. The main dimension ordering functions are:

mv — moves a particular dimension somewhere else in the dimension list
xchg — exchanges two dimensions in the dimension list, leaving the rest alone
reorder — allows wholesale mixing of the dimensions
clump — clumps together two or more small dimensions into one larger one
squeeze — eliminates any dimensions of size 1

🔮 Physical and Dummy Dimensions

document Perl level broadcasting (see "broadcast_define" in PDL::Core)
broadcastids (see "Explicit broadcasting")
slice in PDL::Slices is crucial, useful in niches are lags in PDL::Slices and splitdim in PDL::Slices

🎀 Indexing and broadcasting with PDL 🎀

A lot of the flexibility and power of PDL relies on the indexing and looping features of the Perl extension. Indexing allows access to the data of an ndarray in a very flexible way. Broadcasting provides efficient implicit looping functionality (since the loops are implemented as optimized C code). Optimized C behind the scenes = the real MVP. 🔥

ndarrays are Perl objects that represent multidimensional arrays and operations on those. In contrast to simple Perl @x style lists the array data is compactly stored in a single block of memory thus taking up a lot less memory and enabling use of fast C code to implement operations (e.g. addition, etc) on ndarrays. Less memory, more speed, zero notes. 🫶

👶 ndarrays can have children

Central to many of the indexing capabilities of PDL are the relation of "parent" and "child" between ndarrays. Many of the indexing commands create a new ndarray from an existing ndarray. The new ndarray is the "child" and the old one is the "parent". The data of the new ndarray is defined by a transformation that specifies how to generate (compute) its data from the parent's data.

Another way to interpret the ndarrays created by our indexing commands is to view them as a kind of intelligent pointer that points back to some portion or all of its parent's data. Therefore, it is not surprising that the parent's data (or a portion of it) changes when manipulated through this "pointer". We will further illustrate the pointer/dataflow analogies in the context of some of the examples later on. Stay tuned bestie, it gets good. 💅

There are two different implementations of this "smart pointer" relationship:

the first one, which is a little slower but works for any transformation, is simply to do the transformation forwards and backwards as necessary.
the other is to consider the child ndarray a "virtual" ndarray, which only stores a pointer to the parent and access information so that routines which use the child ndarray actually directly access the data in the parent. If the virtual ndarray is given to a routine which cannot use it, PDL transparently physicalizes the virtual ndarray before letting the routine use it.

Currently (still true as of 2.086) all transformations which are "affine", i.e. the indices of the data item in the parent ndarray are determined by a linear transformation (+ constant) from the indices of the child ndarray result in virtual affine, or "vaffine" ndarrays. All other indexing routines (e.g. ->index(...)) result in physical ndarrays. All routines compiled by PP can accept affine ndarrays, except those routines that pass pointers to external library functions, and those will simply be "physicalised" (memory allocated, data copied in). There will then be dataflow back and forth, in a lazy-evaluation way. See PDL::Dataflow for more.

⚠️ Note that whether something is affine or not does not affect the semantics of what you do in any way: both perl $x->index(...) .= 5; $x->slice(...) .= 5; change the data in $x. The affinity does, however, have a significant impact on memory usage and performance. So like, same result, different budget. 💀

🔪 Slicing ndarrays

Probably the most important application of the concept of parent/child ndarrays is the representation of rectangular slices of a physical ndarray by a virtual ndarray. Having talked long enough about concepts let's get more specific. Suppose we are working with a 2D ndarray representing a 5x5 image (it's unusually small so that we can print it without filling several screens full of digits).

pdl> $im = sequence(5,5)
pdl> p $im

[
 [ 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]
]

pdl> help vars
PDL variables in package main::

Name         Type   Dimension       Flow  State          Mem
----------------------------------------------------------------
$im          Double D [5,5]                P            0.20Kb

The interactive perldl shells that come with PDL have a help vars command that provides information about ndarrays. Use help help for more.

Now suppose we want to create a 1-D ndarray that just references one line of the image, say line 2; or an ndarray that represents all even lines of the image (imagine we have to deal with even and odd frames of an interlaced image due to some peculiar behaviour of our frame grabber). As another frequent application of slices we might want to create an ndarray that represents a rectangular region of the image with top and bottom reversed. All these effects (and many more) can be easily achieved with the powerful slice function:

pdl> $line = $im->slice(':,(2)')
pdl> $even = $im->slice(':,1:-1:2')
pdl> $area = $im->slice('3:4,3:1')
pdl> help vars  # or just PDL->vars
PDL variables in package main::

Name         Type   Dimension       Flow  State          Mem
----------------------------------------------------------------
$even        Double D [5,2]                -C           0.00Kb
$im          Double D [5,5]                P            0.20Kb
$line        Double D [5]                  -C           0.00Kb
$area        Double D [2,3]                -C           0.00Kb

All three "child" ndarrays are children of $im or in the other (largely equivalent) interpretation pointers to data of $im. Operations on those virtual ndarrays access only those portions of the data as specified by the argument to slice. So we can just print line 2:

pdl> p $line
[10 11 12 13 14]

Also note the difference in the "Flow State" of $area above and below:

pdl> p $area
pdl> help $area
This variable is Double D [2,3]                VC           0.00Kb

Watch this next bit, bestie, because it proves that $im and $line are actually behaving like pointer-ish objects — or, in dataflow speak, the changes you make to $line's data get propagated all the way back up to $im:

pdl> $im++
pdl> p $line
[11 12 13 14 15]
pdl> $line += 2
pdl> p $im

[
 [ 1  2  3  4  5]
 [ 6  7  8  9 10]
 [13 14 15 16 17]
 [16 17 18 19 20]
 [21 22 23 24 25]
]

Okay so notice how assignment operations on a child virtual ndarray change the parent physical ndarray and vice versa — BUT the basic = assignment does not do that, you need .= for that effect. (Reasons below, keep reading, it's giving ✨plot twist✨.) Those virtual child ndarrays are basically "live links" to the "original" parent ndarray. Like I said earlier, you can think of them sort of like a C-pointer. But unlike a C-pointer they're carrying so much more information. First: they spell out the structure of the data they represent (aka the dimensionality of the new ndarray). Second: they spell out how to build that structure out of the parent's data (how that actually works is buried deep in PDL's internals and, ngl, you don't need to know (unless you wanna hack the core someday or become a full PDL guru — for a definition of this strange creature see PDL::Internals)).

So the previous examples showed you typical usage of the slice function. And since slicing is that important, here's the full breakdown of the syntax for the string argument you pass to slice:

$vpdl = $x->slice('ind0,ind1...')

where ind0 says what to do with index No 0 of the ndarray $x, etc. Every element of that comma-separated list can take one of these forms:

':' — use the whole dimension. 💯
'n' — use only index n. That index's dimension in the resulting virtual ndarray is 1. An example with these first two index formats:

pdl> $column = $im->slice('2,:')
pdl> $row = $im->slice(':,0')
pdl> p $column

[
 [ 3]
 [ 8]
 [15]
 [18]
 [23]
]

pdl> p $row

[
 [1 2 3 4 5]
]

pdl> help $column
This variable is Double D [1,5]                VC           0.00Kb

pdl> help $row
This variable is Double D [5,1]                VC           0.00Kb

'(n)' — use only index n. This dimension gets removed from the resulting ndarray (relying on the fact that a dimension of size 1 can always be removed). The difference between this and the previous case matters a lot in assignments, where the left and right hand sides have to have appropriate dimensions.

pdl> $line = $im->slice(':,(0)')
pdl> help $line
This variable is Double D [5]                  -C           0.00Kb

pdl> p $line
[1 2 3 4 5]

Spot the difference from the example before? 🧐

'n1:n2' or 'n1:n2:n3' — take the range of indices from n1 to n2, or (second form) take the range from n1 to n2 with step n3. An example of this format is that earlier definition of the sub-image made of even lines.

pdl> $even = $im->slice(':,1:-1:2')

This example also shows that negative indices work just like they do for normal Perl style arrays — counting backwards from the end of the dimension. If n2 is smaller than n1 (in the example -1 is equivalent to index 4) then the elements in the virtual ndarray are effectively reverted with respect to the parent. 🔄

'*[n]' — add a dummy dimension. Its size is 1 by default, or equal to n if you give the optional numerical argument.

Okay so this one is genuinely a bit strange on first sight. What even is a dummy dimension? A dummy dimension inserts a dimension where there wasn't one before. How?? Well, when the new dimension has size 1 it's easy to explain through the way you can identify a vector (with m elements) with a (1,m) or (m,1) matrix. Same obviously holds for higher dimensional objects. The more interesting case is a dummy dimension of size greater than one (e.g. slice('*5,:')). This works the same way as a call to the dummy function creating a new dummy dimension. So read on and check its explanation below. 👀

Generalised diagonals

With diagonal in PDL::Slices you can make generalised diagonals. Best way to explain is probably an example — here we make an ndarray that refers to the elements along the space diagonal of its parent ndarray (a cube):

$cube = zeroes(5,5,5);
$sdiag = $cube->diagonal(1..3);

That command creates a virtual ndarray representing the diagonal along the parent's dimension no. 0, 1 and 2, and makes its dimension 0 (the only dimension) out of it. You use the extended syntax if the dimension sizes of the parent dimensions you want to build the diagonal from have different sizes, or if you want to reverse the sequence of elements in the diagonal, e.g.

$rect = zeroes(12,3,5,6,2);
$vpdl = $rect->slice('2:7,0:1,(4),5:4'); # now 6,2,2,2
$vpdl = $vpdl->diagonal(1..3); # now 6,2

So the elements of $vpdl will then relate to those of its parent in a way we can express as:

vpdl(i,j) = rect(i+2,j,4,5-j,j)       0<=i<6, 0<=j<2

There are different kinds of assignments in PDL

The examples above already showed that virtual ndarrays can be used to operate on or access portions of a parent ndarray's data. They can also be used as lvalues in assignments (as that ++ in some of the examples already demonstrated 💅). For explicit assignments to the data represented by a virtual ndarray you have to use the overloaded .= operator (which in this context we call propagated assignment). So why can't you use the normal = operator?

Okay so, you definitely can still use the = operator, it just wouldn't do what you want. In recent Perls, = can be overloaded just like other assignment operators, but for back-compatibility reasons PDL uses = for a reference copy, so the variable representing the virtual ndarray (a reference to a blessed thingy) would, after the assignment, just contain a reference to another blessed thingy which would behave toward future assignments as a "physical" copy of the original rvalue.

E.g.

pdl> $line = $im->slice(':,(2)')
pdl> $line = zeroes(5);
pdl> $line++;
pdl> p $im

[
 [ 1  2  3  4  5]
 [ 6  7  8  9 10]
 [13 14 15 16 17]
 [16 17 18 19 20]
 [21 22 23 24 25]
]

pdl> p $line
[1 1 1 1 1]

But using .=:

pdl> $line = $im->slice(':,(2)')
pdl> $line .= zeroes(5)
pdl> $line++
pdl> p $im

[
 [ 1  2  3  4  5]
 [ 6  7  8  9 10]
 [ 1  1  1  1  1]
 [16 17 18 19 20]
 [21 22 23 24 25]
]

pdl> print $line
[1 1 1 1 1]

Also, you can substitute

pdl> $line .= 0;

for the assignment above (the zero gets converted to a scalar ndarray, with no dimensions, so it can be assigned to any ndarray).

A cute feature in recent perl versions is lvalue subroutines (i.e., versions 5.6.x and higher, including all perls currently supported by PDL). That lets you use the slicing syntax on both sides of the assignment:

pdl> $im->slice(':,(2)') .= xvals(float, 5)

⚠️ Heads up bestie: there can be a problem with assignments like this when the lvalue and rvalue ndarrays refer to overlapping portions of data in the parent ndarray:

# revert the elements of the first line of $x
$x->slice(':,(1)') .= $x->slice('-1:0,(1)');

Currently, the parent data on the right side of the assignment is not copied before the (internal) assignment loop proceeds. So the outcome of this assignment will depend on the sequence in which elements get assigned and will almost certainly not do what you wanted. The semantics are currently undefined for now and liable to change anytime. 😭 To get the behaviour you actually want, use

$x->slice(':,(1)') .= $x->slice('-1:0,(1)')->copy;

which makes a physical copy of the slice, or

$x->slice(':,(1)') .= $x->slice('-1:0,(1)')->sever;

which returns the same slice but first severs the connection of the slice to its parent. ✂️✨

Other functions that manipulate dimensions

Having gone this deep on the slice function, it should be noted that it's not the only PDL indexing function. There are additional indexing functions that are also useful (especially in the context of broadcasting, which we'll talk about later). Here's a list and some examples of how to use them.

dummy — inserts a dummy dimension of the size you specify (default 1) at the chosen location. Can't wait to hear how that's achieved? Well, all elements with index (X,x,Y) (0<=x<size_of_dummy_dim) just map to the element with index (X,Y) of the parent ndarray (where X and Y refer to the group of indices before and after the location where the dummy dimension was inserted).

This example calculates the x coordinate of the centroid of an image — later we'll learn that we didn't actually need the dummy dimension thanks to the magic of implicit broadcasting, but using dummy dimensions the code would also work in a broadcast-less world; though once you've worked with PDL broadcasting you wouldn't want to live without them again. 🫶

# centroid
($xd,$yd) = $im->dims;
$xc = sum($im*xvals($xd)->dummy(1,$yd))/sum($im);

Let's explain how that works in a little more detail. First, the product:

$xvs = xvals($xd);
print $xvs->dummy(1,$yd);      # repeat the line $yd times
$prod = $im*xvs->dummy(1,$yd); # form the pixel-wise product with
                               # the repeated line of x-values

The rest is then summing the results of the pixel-wise product together and normalizing with the sum of all pixel values in the original image, thereby calculating the x-coordinate of the "center of mass" of the image (interpreting pixel values as local mass) — which is known as the centroid of an image. 💖

Okay so this next bit is a (memory-consumption-wise) super cheap glow-up from grey-scale to RGB — every pixel now holds a whole triple of values instead of just one scalar. And the way that works out is genuinely iconic: all three values in that triple are literally the same for a grey image, so our trick maps all three members of the triple onto the same source element. No cap.

# a cheap grey-scale to RGB conversion
$rgb = $grey->dummy(0,3)

Tragically though, bestie, this trick will not turn your old B/W photos into color ones the way you'd want. :( It's giving scam, honestly.

Heads up: the memory usage of ndarrays with dummy dimensions is especially sensitive to the internal representation. If the ndarray can be represented as a virtual affine ("vaffine") ndarray, only the control structures get stored — slay, basically free. But if $y in

$x = zeroes(10000);
$y = $x->dummy(1,10000);

gets made physical by some routine, you will find your program's memory usage has suddenly grown by 100Mb. Not a vibe. 💀

✨ "diagonal" ✨

This one replaces two dimensions (which have to be the same size, no exceptions, periodt) with one dimension that references all the elements along the "diagonal" along those two dimensions. Here are two examples that should feel very familiar to anyone who's ever touched linear algebra 🎀

First up, making a unity matrix:

# unity matrix
$e = zeroes(float, 3, 3); # make everything zero
$e->diagonal(0,1) .= 1; # set the elements along the diagonal to 1
print $e;

Or the other diagonal:

$e->slice(':-1:0')->diagonal(0,1) .= 2;
print $e;

(Did you clock how we used the slice function to reverse the sequence of lines before setting the diagonal of the new child, thereby setting the cross diagonal of the parent? 🦋 Look at us go.)

Or, if you want a mapping from the space of diagonal matrices to the field over which the matrices are defined — that's the trace of a matrix, babes:

# trace of a matrix
$trace = sum($mat->diagonal(0,1));  # sum all the diagonal elements

✨ "xchg" and "mv" ✨

xchg exchanges — or "transposes" — the two specified dimensions. Straightforward moment:

# transpose a matrix (without explicitly reshuffling data and
# making a copy)
$prod = $x x $x->xchg(0,1);

$prod should now be pretty close to the unity matrix if $x is an orthogonal matrix. And often xchg is used in the context of broadcasting, but more on that later 👀

mv works in a similar fashion. It moves a dimension (specified by its number in the parent) to a new position in the new child ndarray:

$y = $x->mv(4,0);  # make the 5th dimension of $x the first in the
                   # new child $y

So what's the difference between xchg and mv? xchg only changes the position of two dimensions with each other, whereas mv inserts the first dimension into the place of the second, moving the other dimensions around accordingly. 💅

✨ "clump" ✨

clump collapses several dimensions into one. Its only argument says how many dimensions of the source ndarray should be collapsed (starting from the first). Our (admittedly unrealistic) example is a 3D ndarray holding data from a stack of image files you just read in. Except — plot twist — the data from each image actually represents a 1D time series, and it only got arranged that way because it was digitized with a frame grabber. So to get it back as an array of time sequences, you say:

pdl> $seqs = $stack->clump(2)
pdl> help vars
PDL variables in package main::

Name         Type   Dimension       Flow  State          Mem
----------------------------------------------------------------
$seqs        Double D [8000,50]            -C           0.00Kb
$stack       Double D [100,80,50]          P            3.05Mb

Unrealistic as it may seem, our confocal microscope software genuinely writes data (sometimes) this way 💀. But more often you use clump to achieve a certain effect when using implicit or explicit broadcasting.

🌈 Calls to indexing functions can be chained 🌈

As you might've noticed in some of the examples above, calls to the indexing functions chain together beautifully, since all of these functions return a newly created child object. Main character energy, honestly.

But when you're doing extensive index manipulations in a chain, bestie, please keep track of what you're doing, like:

$x->xchg(0,1)->mv(0,4)

This moves the dimension 1 of $x to position 4 — because when the second command is executed, the original dimension 1 has already been moved to position 0 of the new child that calls the mv function. I think you get the idea (in spite of my convoluted explanations 🙈).

💜 Propagated assignments (`.=`) and dummy dimensions 💜

Okay here's a subtlety related to indexing: assignment to ndarrays containing dummy dimensions of size greater than 1. These assignments (using .=) are forbidden, because several elements of the lvalue ndarray point to the same element of the parent. Which means the value of those parent elements is potentially ambiguous — it'd depend on the sequence in which the implementation makes the assignments to elements. Chaos. Not the cute kind. 😭

So an assignment like this:

$x = pdl [1,2,3];
$y = $x->dummy(1,4);
$y .= yvals(3,4);

can produce unexpected results, and the results are explicitly undefined, and will (as of 2.086) throw an error. So like. Don't. 🚫

From the point of view of dataflow, introducing greater-size-than-one dummy dimensions is regarded as an irreversible transformation (similar to the terminology in thermodynamics 🌡️) which precludes backward propagation of assignment to a parent — a parent you had explicitly requested using the .= assignment. A similar problem to watch out for shows up in the context of broadcasting, where sometimes dummy dimensions get created implicitly during the broadcast loop (see below 👇).

🧚 Reasons for the parent/child (or "pointer") concept 🧚

It is very memory efficient, with only metadata existing in the child. Skinny legend. ✨
It provides a very flexible and powerful way of accessing portions of ndarray data. 👑

🎀 How to make things physical again 🎀

A "virtual affine" (or "vaffine") ndarray will be automatically "physicalised" — meaning memory gets allocated and its data copied in — if that is requested with the make_physical in PDL::Core method, if it is severed, or if it is passed to an operation as a parameter that is marked as [phys]. See PDL::PP for more.

It will still have dataflow, except in the sever case. 🫶

📣 Broadcasting 📣

So in the previous bit on indexing we already dropped the term here and there, but now it's really time to talk explicitly about "broadcasting" with ndarrays. Within the framework of PDL it could probably be loosely defined as an implicit looping facility. Implicit because you don't specify anything like enclosing for-loops — instead the loops are automatically (or 'magically' ✨) generated by PDL based on the dimensions of the ndarrays involved.

Which should give you a first idea of why the index/dimension manipulating functions from the previous paragraphs are especially important and useful in the context of broadcasting. The other ingredient for broadcasting (apart from the ndarrays involved) is a function that is broadcasting aware — generally, these are PDL::PP compiled functions — and that the ndarrays are "broadcast" over.

So much for the terminology. Now let's try to shed some light on what it all actually means. 🔥

🌸 Implicit broadcasting — a first example 🌸

There are two slightly different variants of broadcasting. We're starting with what we call "implicit broadcasting". Let's pick a practical example involving looping a function over many elements of an ndarray. Suppose we have an RGB image we want to convert to grey-scale. The RGB image is represented by a 3-dim ndarray im(3,x,y), where the first dimension holds the three color components of each pixel, and x and y are the width and height of the image, respectively.

Next we need to specify how to convert a color-triple at a given pixel into a grey-value (and to be a realistic example, it should represent the relative intensity with which our color-insensitive eye cells would detect that color — that's what we'd call a natural conversion from color to grey-scale). An approximation that works quite well is to compute the grey intensity from each RGB triplet (r,g,b) as a weighted sum:

grey-value = 77/256*r + 150/256*g + 29/256*b =
    inner([77,150,29]/256, [r,g,b])

where that last form indicates we can write this as an inner product of the 3-vector comprising the weights for red, green and blue components with the 3-vector containing the color components. 💅

Traditionally, we might have written a function like the following to process the whole image:

my @dims = $im->dims;
# here normally check that first dim has correct size (3), etc
$grey = zeroes(@dims[1,2]); # make the ndarray for the resulting grey image
$w = pdl [77,150,29] / 256; # the vector of weights
for ($j=0; $j<dims[2]; $j++) {
   for ($i=0; $i<dims[1]; $i++) {
       # compute the pixel value
       $tmp = inner($w,$im->slice(':,(i),(j)'));
       set($grey,$i,$j,$tmp); # and set it in the grey-scale image
   }
}

Now we write the same thing using broadcasting (noting that inner is a broadcasting aware function defined in the PDL::Primitive package):

$grey = inner($im, pdl([77,150,29]/256));

We have ended up with a one-liner 😭 that automatically creates the ndarray $grey with the right number and size of dimensions and performs the loops automatically (these loops are implemented as fast C code in the internals of PDL). Ate and left no crumbs, fr fr.

Well… we still owe you an explanation of how this 'magic' is achieved. 🪄

✨ How does the example work? ✨

okay so bestie, first thing you gotta know: every function that's broadcasting-aware (and these are, without exception, functions compiled from concise descriptions by PDL::PP — we call them PP-functions from here on, periodt) expects a defined minimum number of dimensions from each of its ndarray arguments. we call those the core dimensions 💖.

so like, the inner function expects two one-dimensional (input) parameters, and from those it calculates a zero-dimensional (output) parameter. we write that symbolically as inner((n),(n),[o]()) and call it inner's signature, where n is the size of that dimension. and when n is equal in the first and second parameter, that means those dimensions HAVE to be the same size in any call — no cap, that's a hard requirement.

another example, bestie: take the outer product. it takes two 1D vectors and generates a 2D matrix, written symbolically as outer((n),(m),[o](n,m)). the [o] in both examples means that argument (here, the third one) is an output argument 🌈. and in that outer example, the dimensions of the first and second argument don't have to agree — but you can see how they determine the size of the two dimensions of the output ndarray. iconic.

now HERE is the moment when broadcasting finally enters the chat 💅. if you call PP-functions with ndarrays that have more than the required core dimensions, then the first dimensions of the ndarray arguments get used as the core dimensions, and all the extra additional dimensions get broadcast over. let's demonstrate this with our example from above:

$grey = inner($im,$w); # w is the weight vector from above

in this case $w is 1D so it only supplied the core dimension. $im is 3D, more specifically (3,x,y). the first dimension (size 3) is the required core dimension that matches (as required by inner) the first (and only) dimension of $w. the second dimension is the first broadcast dimension (size x) and the third is here the second broadcast dimension (size y). the output ndarray is automatically created (as requested by setting $grey to null before invocation). the output dimensions come from appending the loop dimensions (here (x,y)) to the core output dimensions (here 0D), which gives you the final dimensions of the auto-created ndarray (here 0D+2D=2D, yielding a 2D output of size (x,y)).

so that command up there calls the core functionality that computes the inner product of two 1D vectors x*y times with $w and all 1D slices of the form (':,(i),(j)') of $im, and it sets the respective elements of the output ndarray $grey(i,j) to the result of each computation. we could write that symbolically as:

$grey(0,0) = f($w,$im(:,(0),(0)))
$grey(1,0) = f($w,$im(:,(1),(0)))
    .
    .
    .
$grey(x-2,y-1) = f($w,$im(:,(x-2),(y-1)))
$grey(x-1,y-1) = f($w,$im(:,(x-1),(y-1)))

but bestie, this is all done automatically by PDL — you're not writing ANY explicit Perl loops. not a single one. we see the command really does create an output ndarray with the right dimensions and sets the elements to the result of the computation for each pixel of the input image. slay.

when even MORE ndarrays and extra dimensions get involved, things get a bit more complicated. so first we'll give you the general rules for how broadcast dimensions depend on the dimensions of input ndarrays — this lets you figure out the dimensionality of an auto-created output ndarray (for any given set of input ndarrays and core dimensions of the PP-function in question). ngl, the general rules are probably gonna look a little confusing on first sight, so we'll illustrate the usage with a bunch of further examples (which will hopefully also show you that there really are SO many practical situations where broadcasting is extremely handy).

🎀 A call for coding discipline 🎀

before we point out the other technical details of broadcasting, please note this call for programming discipline when using broadcasting:

In order to preserve human readability, PLEASE comment any nontrivial expression in your code involving broadcasting. Most importantly, for any subroutine, include information at the beginning about what you expect the dimensions to represent (or ranges of dimensions).

as a warning 👀, look at this undocumented function and try to guess what might be going on:

sub lookup {
  my ($im,$palette) = @_;
  $palette->xchg(0,1)->index($im->indx->dummy(0,$palette->dim(0)));
}

would you agree that it might be difficult to figure out expected dimensions, purpose of the routine, etc? 😭 (if you wanna find out what this little piece of code does, see below)

👑 How to figure out the loop dimensions 👑

there are a couple of rules that let you figure out the number and size of loop dimensions (and whether the size of your input ndarrays complies with the broadcasting rules). dimensions of any ndarray argument get broken down into two groups here: Core dimensions (as defined by the PP-function — see Appendix A for a list of PDL primitives) and extra dimensions, which is all the remaining dimensions of that ndarray. for example, calling a function func with the signature func((n,m),[o](n)) with an ndarray $x(2,4,7,1,3) as f($x,($o = null)) results in the semantic splitting of x's dimensions into: core dimensions (2,4) and extra dimensions (7,1,3).

R0 — Core, or "active", dimensions are identified with the first N dimensions of the respective ndarray argument (and are required). Any further dimensions are extra dimensions and used to determine the loop dimensions.
R1 — The number of (implicit) loop dimensions is equal to the maximal number of extra dimensions taken over the set of ndarray arguments.
R2 — The size of each of the loop dimensions is derived from the size of the respective dimensions of the ndarray arguments. The size of a loop dimension is given by the maximal size found in any of the ndarrays having this extra dimension.
R3 — If an ndarray doesn't have a loop dimension, in the broadcast loop this ndarray is "promoted" as if having a dummy dimension of size equal to the size of that loop dimension.
R4 — For all ndarrays that have a given extra dimension the size must be equal to the size of the loop dimension (as determined by the previous rule) or 1; otherwise you raise a runtime exception 🚨. If the size of the extra dimension in an ndarray is 1 it is implicitly treated as a dummy dimension of size equal to ("inflated" to) that loop dim size when performing the broadcast loop.
R5 — If output auto-creation is used (by setting the relevant ndarray to PDL->null before invocation) the number of dimensions of the created ndarray is equal to the sum of the number of core output dimensions + number of loop dimensions. The size of the core output dimensions is derived from the relevant dimension of input ndarrays (as specified in the function definition) and the sizes of the other dimensions are equal to the size of the loop dimension it is derived from. The automatically created ndarray will be physical (unless dataflow is in operation).

in this context, note that you can run into the problem with assignment to ndarrays containing greater-than-one dummy dimensions (see above). even though your output ndarray(s) didn't contain any dummy dimensions in the first place, they may end up with implicitly created dummy dimensions according to R3. lowkey wild, i know.

as an example, suppose we have a (here unspecified) PP-function with the signature:

func((m,n),(m,n,o),(m),[o](m,o))

and you call it with 3 ndarrays $x(5,3,10,11), $y(5,3,2,10,1,12), and $z(5,1,11,12) as:

func($x,$y,$z,($d=null))

then the number of loop dimensions is 3 (by R0+R1 from $y and $z) with sizes (10,11,12) (by R2). the two output core dimensions are (5,2) (from the signature of func), resulting in a 5-dimensional output ndarray $c of size (5,2,10,11,12) (see R5), and (the automatically created) $d is derived from ($x,$y,$z) in a way that can be expressed in pdl pseudo-code as:

$d(:,:,i,j,k) .= func($x(:,:,i,j),$y(:,:,:,i,0,k),$z(:,0,j,k))
   with 0<=i<10, 0<=j<=11, 0<=k<12

if we analyze the color to grey-scale conversion again with these rules in mind, we notice another great advantage of implicit broadcasting 💖. we can call the conversion with an ndarray representing a pixel (im(3)), a line of rgb pixels (im(3,x)), a proper color image (im(3,x,y)) or a whole stack of RGB images (im(3,x,y,z)). as long as $im is of the form (3,...), the automatically created output ndarray will contain the right number of dimensions and contain the intensity data as we expect, since the loops have been implicitly performed thanks to implicit broadcasting. you can easily convince yourself: calling with a color pixel $grey is 0D, with a line it turns out 1D grey(x), with an image we get grey(x,y), and finally we get a converted image stack grey(x,y,z). chef's kiss 👑

let's fill these general rules with some more life by going through a couple more examples. the reader may try to figure out equivalent formulations with explicit for-looping and compare the flexibility of those routines using implicit broadcasting to the explicit formulation. furthermore — especially when using several broadcast dimensions — it's a useful exercise to check the relative speed by doing some benchmark tests (which we still have to do, lol).

first in the row is a slightly reworked centroid example, now coded with broadcasting in mind.

# broadcast mult to calculate centroid coords, works for stacks as well
$xc = sumover(($im*xvals($im->dim(0))->clump(2)) /
      sumover($im->clump(2));

let's analyze what's going on step by step. first the product:

$prod = $im*xvals($im->dim(0))

this will actually work for $im being one, two, three, and higher dimensional. if $im is one-dimensional it's just an ordinary product (in the sense that every element of $im gets multiplied with the respective element of xvals(...)). if $im has more dimensions, further broadcasting is done by adding appropriate dummy dimensions to xvals(...) according to R3. more importantly, the two sumover operations show a first example of how to make use of the dimension-manipulating commands. a quick look at sumover's signature will remind you that it will only "gobble up" the first dimension of a given input ndarray. but what if we want to really compute the sum over all elements of the first two dimensions? well, nothing keeps us from passing a virtual ndarray into sumover — in this case it's formed by clumping the first two dimensions of the "parent ndarray" into one. from the point of view of the parent ndarray, the sum is now computed over the first two dimensions, just as we wanted, even though sumover has just done the job as specified by its signature. got it? 💅

Another lil finesse of writing the code like that: we intentionally used sumover($pdl->clump(2)) instead of sum($pdl) so that we can either pass just an image (x,y) or a stack of images (x,y,t) into this routine and get back either just one x-coordinate or a vector of x-coordinates (of size t). One function, two vibes, no extra work. ✨ iykyk ✨

Another set of common operations are what one could call "projection operations". These operations take an N-D ndarray as input and return a (N-1)-D "projected" ndarray. These operations are often performed with functions like sumover, prodover, minimum and maximum. 🎀

Using images as our example again 🦋 we might want to calculate the maximum pixel value for each line of an image or image stack. We know how to do that:

# maxima of lines (as function of line number and time)
maximum($stack,($ret=null));

But bestie, what if you want to calculate maxima per column, when implicit broadcasting always applies the core functionality to the first dimension and broadcasts over all the others? 😭 How do we make it so the core functionality hits the second dimension instead, and broadcasting goes over the rest? Can you guess it? Yes queen, you guessed it — we make a virtual ndarray that has the second dimension of the "parent ndarray" as its first dimension, using the "mv" command. 💅

# maxima of columns (as function of column number and time)
maximum($stack->mv(1,0),($ret=null));

and calculating all the sums of sub-slices over the third dimension is now almost too easy:

# sums of pixels in time (assuming time is the third dim)
sumover($stack->mv(2,0),($ret=null));

Finally, if you want to apply the operation to all elements (like max over all elements or sum over all elements) regardless of the dimensions of the ndarray in question, "clump" comes in handy. As an example, look at a definition of "sum" (summarised from Basic/Ufunc/ufunc.pd):

sub sum {
  PDL::Ufunc::sumover($name->clump(-1),($tmp=null));
  return $tmp; # return a 0D ndarray
}

We have already mentioned that all basic operations support broadcasting, and assignment is no exception. So here are a couple of broadcasted assignments 🌸

pdl> $im = zeroes(byte, 10,20)
pdl> $line = exp(-rvals(10)**2/9)
# broadcasted assignment
pdl> $im .= $line      # set every line of $im to $line
pdl> $im2 .= 5         # set every element of $im2 to 5

By now you probably see how it works and what it does, don't you? 👑

To finish the examples in this paragraph, here is a function to create an RGB image from what is called a palette image. The palette image consists of two parts: an image of indices into a color lookup table, and the color lookup table itself. [ describe how it works ] We are going to use a PP-function we haven't encountered yet in the previous examples. It is the aptly named index function, signature "((n),(),[o]())" (see Appendix A), with the core functionality that "index(pdl (0,2,4,5),2,($ret=null))" will return the element with index 2 of the first input ndarray. In this case, $ret will contain the value 4. So here is the example:

# a broadcasted index lookup to generate an RGB, or RGBA or YMCK image
# from a palette image (represented by a lookup table $palette and
# an color-index image $im)
# you can say just dummy(0) since the rules of broadcasting make it fit
pdl> index($palette->xchg(0,1),
              $im->indx->dummy(0,$palette->dim(0)),
              $res=null);

Let's go through it and explain the steps involved. Assuming we are dealing with an RGB lookup-table, $palette is of size (3,x). First we exchange the dimensions of the palette so that looping is done over the first dimension of $palette (of size 3, that represent r, g, and b components). Now looking at $im, we add a dummy dimension of size equal to the length of the number of components (in the case we are discussing here we could have just used the number 3 since we have 3 color components). We can use a dummy dimension since for red, green and blue color components we use the same index from the original image — e.g. assuming a certain pixel of $im had the value 4, then the lookup should produce the triple

[palette(0,4),palette(1,4),palette(2,4)]

for the new red, green and blue components of the output image. 🎨 Hopefully by now you have some sort of idea what the above piece of code is supposed to do (it is often actually quite complicated to describe in detail how a piece of broadcasting code works; just go ahead and experiment a bit to get a better feeling for it). 🙈

If you have read the broadcasting rules carefully, then you might have noticed that we didn't have to explicitly state the size of the dummy dimension that we created for $im; when we create it with size 1 (the default) the rules of broadcasting make it automatically fit to the desired size (by rule R4 — in our example the size would be 3, assuming a palette of size (3,x)). Since situations like this do occur often in practice, this is actually why rule R4 has been introduced (the part that makes dimensions of size 1 fit to the broadcast loop dim size). So we can just say:

pdl> index($palette->xchg(0,1),$im->indx->dummy(0),($res=null));

Again, you can convince yourself that this routine will create the right output if called with a pixel ($im is 0D), a line ($im is 1D), an image ($im is 2D), ..., an RGB lookup table (palette is (3,x)) and RGBA lookup table (palette is (4,x), see e.g. OpenGL). This flexibility is achieved by the rules of broadcasting, which are made to do the right thing in most situations. 🦋

To wrap it all up once again, the general idea is as follows. 💖 If you want to achieve looping over certain dimensions and have the core functionality applied to another specified set of dimensions, you use the dimension manipulating commands to create a (or several) virtual ndarray(s) so that from the point of view of the parent ndarray(s) you get what you want (always having the signature of the function in question and R1-R5 in mind!). Easy, isn't it? ✨

See "dice" in PDL::Slices, and "dice_axis" in PDL::Slices for more indexed lookups, and "indexND" in PDL::Slices, "whereND" in PDL::Slices, "whichND" in PDL::Slices for even more powerful constructs. 🔥

🎀 Output auto-creation and PP-function calling conventions 🎀

At this point we have to divert to some technical detail that has to do with the general calling conventions of PP-functions and the automatic creation of output arguments. Basically, there are two ways of invoking PDL routines, namely

$result = func($x,$y);

and

func($x,$y,$result);

If you are only using implicit broadcasting, then the output variable can be automatically created by PDL. 💅 You flag that to the PP-function by setting the output argument to a special kind of ndarray that is returned from a call to the function "PDL->null", which returns an essentially "empty" ndarray (for those interested in details, there is a flag in the C pdl structure for this). The dimensions of the created ndarray are determined by the rules of implicit broadcasting: the first dimensions are the core output dimensions to which the broadcasting dimensions are appended (which are in turn determined by the dimensions of the input ndarrays as described above). So you can say

func($x,$y,($result=PDL->null));

$result = func($x,$y)

which are exactly equivalent. 👑

🚨 Be warned that you can not use output auto-creation when using explicit broadcasting (for reasons explained in the following section on explicit broadcasting, the second variant of broadcasting).

In "tight" loops you probably want to avoid the implicit creation of a temporary ndarray in each step of the loop that comes along with the "functional" style, but rather say

# create output ndarray of appropriate size only at first invocation
$result = null;
for (0...$n) {
     func($x,$y,$result); # in all but the first invocation $result
     func2($y);           # is defined and has the right size to
                          # take the output provided $y's dims don't change
     twiddle($result,$x); # do something from $result to $x for iteration
}

The take-home message of this section once more: be aware of the limitation on output creation when using explicit broadcasting. 🫶

🌈 Explicit broadcasting 🌈

Having so far only talked about the first flavour of broadcasting, it is now about time to introduce the second variant. 🔥 Instead of shuffling around dimensions all the time and relying on the rules of implicit broadcasting to get it all right, you sometimes might want to specify in a more explicit way how to perform the broadcast loop. It is probably not too surprising that this variant of the game is called explicit broadcasting. Now, before we create the wrong impression: it is not either implicit or explicit; the two flavours do mix. 💜 But more about that later.

The two most used functions with explicit broadcasting are broadcast and unbroadcast. We start with an example that illustrates typical usage of the former:

# code to show that $mat += $line is different from
#                   $mat->broadcast(0) += $line
# explicit broadcasting to add a vector to each column of a matrix
pdl> $mat  = zeroes(4,3)
pdl> $line = pdl (3.1416,2,-2)
pdl> $mat->broadcast(0) += $line

In this example, "$mat->broadcast(0)" tells PDL that you want the second dimension of this ndarray to be broadcast over first, leading to a broadcast loop that can be expressed as

for (j=0; j<3; j++) {
   for (i=0; i<4; i++) {
       mat(i,j) += src(j);
   }
}

"broadcast" takes a list of numbers as arguments which explicitly specify which dimensions to broadcast over first. 🎀 With the introduction of explicit broadcasting, the dimensions of an ndarray are conceptually split into three different groups, the latter two of which we have already encountered: broadcast dimensions, core dimensions and extra dimensions.

Conceptually, it is best to think of those dimensions of an ndarray that have been specified in a call to "broadcast" as being taken away from the set of normal dimensions and put on a separate stack. So assuming we have an ndarray "x(4,7,2,8)", saying

$y = $x->broadcast(2,1)

creates a new virtual ndarray of dimension y(4,8) (which we call the remaining dims) that also has 2 broadcast dimensions of size (2,7). For the purposes of this doc we write that symbolically as y(4,8){2,7}. 💅

And here's the big plot twist compared to the earlier examples where only implicit broadcasting was in play: now the core dimensions are matched against the remaining dimensions, which are not necessarily the first dimensions of the ndarray. So let's lay out exactly how having broadcast dimensions rewrites the rules R1–R5 for broadcast loops (remember, R1–R5 apply to the special case where none of the ndarray arguments has any broadcast dimensions at all).

B0 — Core dimensions are matched against the first n remaining dimensions of the ndarray argument (plot twist: this is different from R1!). Any further remaining dimensions are extra dimensions, and those get used to determine the implicit loop dimensions.
B1a — The number of implicit loop dimensions equals the maximal number of extra dimensions taken over the whole set of ndarray arguments.
B1b — The number of explicit loop dimensions equals the maximal number of broadcast dimensions taken over the whole set of ndarray arguments.
B1c — The total number of loop dimensions = explicit loop dimensions + implicit loop dimensions. In the broadcast loop, explicit loop dimensions get broadcasted over first, then implicit loop dimensions follow.
B2 — The size of each loop dimension comes from the size of the respective dimensions of the ndarray arguments. It's the maximal size found in any ndarray that has this broadcast dimension (for explicit loop dimensions) or this extra dimension (for implicit loop dimensions).
B3 — This rule hits both explicit and implicit loop dimensions. For all ndarrays that have a given broadcast/extra dimension, the size must be equal to the size of the respective explicit/implicit loop dimension or 1 — otherwise you raise a runtime exception 💀. If the size of an ndarray's broadcast/extra dimension is one, it gets implicitly treated as a dummy dimension of size equal to the explicit/implicit loop dimension.
B4 — If an ndarray doesn't have a broadcast/extra dimension corresponding to an explicit/implicit loop dimension, then in the broadcast loop it's treated as if it had a dummy dimension of size equal to the size of that loop dimension.
B4a — All ndarrays that do have broadcast dimensions must have the same number of broadcast dimensions.
B5 — Output auto-creation cannot be used if any of the ndarray arguments has any broadcast dimensions. Otherwise R5 applies.

The same restrictions apply with regard to implicit dummy dimensions (created by application of B4) as already mentioned in the section on implicit broadcasting: if any of the output ndarrays has an (explicit or implicitly created) greater-than-one dummy dimension, a runtime exception will be raised. 🚨

Okay bestie, let's watch these rules do their thing in a generic case. Suppose we have a (here unspecified) PP-function with the signature:

func((m,n),(m),(),[o](m))

and you call it with 3 ndarrays a(5,3,10,11), b(3,5,10,1,12), c(10) and an output ndarray d(3,11,5,10,12) (which here can't be automatically created) like this:

func($x->broadcast(1,3),$y->broadcast(0,3),$c,$d->broadcast(0,1))

From the signature of func and the above call, the ndarrays split into the following groups of core, extra and broadcast dimensions (written in the form pdl(core dims){broadcast dims}[extra dims]):

a(5,10){3,11}[] b(5){3,1}[10,12] c(){}[10] d(5){3,11}[10,12]

With this to help us along (and ngl, it's genuinely helpful to write your arguments down like this whenever you start playing with broadcasting and want to keep track of what's going on), we further deduce that the number of explicit loop dimensions is 2 (by B1b, from $a and $b) with sizes (3,11) (by B2); 2 implicit loop dimensions (by B1a, from $b and $d) of size (10,12) (by B2); and the elements are computed from the input ndarrays in a way that can be expressed in pdl pseudo-code as

for (l=0; l<12; l++)
 for (k=0; k<10; k++)
  for (j=0; j<11; j++)         effect of treating it as dummy dim (index j)
   for (i=0; i<3; i++)                       |
      d(i,j,:,k,l) = func(a(:,i,:,j),b(i,:,k,0,l),c(k))

Ugh, this example was really not easy in terms of bookkeeping 😭. It serves mostly as an example of how to figure out what's going on when you encounter a complicated looking expression. But now it is really time to show that broadcasting is useful by giving some more of our so called "practical" examples. ✨

[ The following examples will need some additional explanations in the future. For the moment please try to live with the comments in the code fragments. ]

Example 1:

*** inverse of matrix represented by eigvecs and eigvals
** given a symmetrical matrix M = A^T x diag(lambda_i) x A
**    =>  inverse M^-1 = A^T x diag(1/lambda_i) x A
** first $tmp = diag(1/lambda_i)*A
** then  A^T * $tmp by broadcasted inner product
# index handling so that matrices print correct under pdl
$inv .= $evecs*0;  # just copy to get appropriately sized output
$tmp .= $evecs;    # initialise, no back-propagation
$tmp->broadcast(0) /= $evals;    #  broadcasted division
# and now a matrix multiplication in disguise
PDL::Primitive::inner($evecs->xchg(0,1)->broadcast(-1,1),
                      $tmp->broadcast(0,-1),
                      $inv->broadcast(0,1));
# alternative for matrix mult using implicit broadcasting,
# first xchg only for transpose
PDL::Primitive::inner($evecs->xchg(0,1)->dummy(1),
                      $tmp->xchg(0,1)->dummy(2),
                      ($inv=null));

Example 2:

# outer product by broadcasted multiplication
# stress that we need to do it with explicit call to my_biop1
# when using explicit broadcasting
$res = zeroes($x->dim(0), $y->dim(0));
my_biop1($x->broadcast(0,-1),$y->broadcast(-1,0),$res->(0,1),"*");
# similar thing by implicit broadcasting with auto-created ndarray
$res = $x->dummy(1) * $y->dummy(0);

Example 3:

# different use of broadcast and unbroadcast to shuffle a number of
# dimensions in one go without lots of calls to ->xchg and ->mv

# use broadcast/unbroadcast to shuffle dimensions around
# just try it out and compare the child ndarray with its parent
$trans = $x->broadcast(4,1,0,3,2)->unbroadcast;

Example 4:

# calculate a couple of bounding boxes
# $bb will hold BB as [xmin,xmax],[ymin,ymax],[zmin,zmax]
# we use again broadcast and unbroadcast to shuffle dimensions around
pdl> $bb = zeroes(double, 2, 3);
pdl> minimum($vertices->broadcast(0)->clump->unbroadcast(1), $bb->slice('(0),:'));
pdl> maximum($vertices->broadcast(0)->clump->unbroadcast(1), $bb->slice('(1),:'));

Example 5:

# calculate a self-rationed (i.e. self normalized) sequence of images
# uses explicit broadcasting and an implicitly broadcasted division
$stack = read_image_stack();
# calculate the average (per pixel average) of the first $n+1 images
$aver = zeroes([stack->dims]->[0,1]);  # make the output ndarray
sumover($stack->slice(":,:,0:$n")->broadcast(0,1),$aver);
$aver /= ($n+1);
$stack /= $aver;  # normalize the stack by doing a broadcasted division
# implicit versus explicit
# alternatively calculate $aver with implicit broadcasting and auto-creation
sumover($stack->slice(":,:,0:$n")->mv(2,0),($aver=null));
$aver /= ($n+1);
#

🌸 Implicit versus explicit broadcasting 🌸

In this paragraph we're going to spill on when explicit broadcasting is preferable over implicit broadcasting and vice versa. But honestly, that framing is lowkey not the best way of putting it, because you already know: the two flavours do mix. So it's more about how to get the best of both worlds and, anyway, in the best of Perl traditions: TIMTOWTDI ! 🫶

[ Sorry, this still has to be filled in in a later release; either refer to above examples or choose some new ones ]

Finally, this may be a good place to justify all the technical detail we have been going on about for a couple pages: why broadcasting?? 💭

Well:

Code that uses broadcasting should be (considerably) faster than code that uses explicit for-loops (or similar Perl constructs) to achieve the same functionality. Especially on supercomputers (with vector computing facilities / parallel processing), PDL broadcasting is implemented in a way that takes advantage of the additional facilities of these machines. 🔥
Furthermore, it's a conceptually simple construct (even though the technical details might get involved at times) and it can greatly reduce the syntactical complexity of PDL code — but keep the admonition for documentation in mind, bestie, fr fr.
Once you're comfortable with the broadcasting way of thinking (and coding), it shouldn't be too difficult to understand code that somebody else wrote (provided they gave you an idea what expected input dimensions are, etc.).

And as a general tip to increase the performance of your code: if you have to introduce a loop into your code, try to reformulate the problem so that you can use broadcasting to perform the loop. (As with anything, there are exceptions to this rule of thumb — but the authors of this document tend to think these are rare cases ;).)

✨ PDL::PP ✨

💖 An easy way to define functions that are aware of indexing and broadcasting 💖

(and the universe and everything 🌈)

PDL::PP is part of the PDL distribution. It is used to generate functions that are aware of indexing and broadcasting rules from very concise descriptions. It can be useful for you if you want to write your own functions, or if you want to interface functions from an external library so that they support indexing and broadcasting (and maybe dataflow as well, see PDL::Dataflow). For further details check PDL::PP. 🎀

📎 Appendix A 💅

okay bestie, this is the receipts section. Here's a little hand-picked pile of signatures from actual PDL primitives, so you can see how many dimensions these PP-compiled functions gobble up — and therefore figure out what's gonna get broadcasted over. Most of these live in the basic ones defined in Primitive.pd ✨

# functions in Primitive.pd
#
axisvalues     ([o](n))
inner      ((n),(n),[o]())
outer      ((n),(m),[o](n,m))
innerwt    ((n),(n),(n),[o]())
inner2         ((m),(m,n),(n),[o]())
inner2d      ((n,m),(n,m),[o]())
inner2t    ((j,n),(n,m),(m,k),[o](j,k))
wtstat     ((n),(n),(),[o],())

# functions in Ufunc.pd
#
sumover    ((n),[o]())
prodover   ((n),[o]())
minimum    (1D,[o])
maximum    (1D,[o])

index      (1D,0D,[o])    in Slices.pd
assgn      ((),())        in Ops.pd

# basic operations
binary operations ((),(),[o]())
unary operations  ((),[o]())

👑 AUTHOR & COPYRIGHT 👑

and that's a wrap on the whole page, bestie!! 🎀 This version came from perl v5.42.3, dated 2026-08-20, and it's PDL::Indexing(3pm) — go forth and broadcast responsibly, you absolute indexing icon 💅✨🦋

📦 where this page comes from

ubuntu stonking: pdl 1:2.106-1build1 · Artistic-1.0-Perl OR GPL-1.0-or-later
💖 girlified by deepseek-ai/DeepSeek-V4.1-Flash on 2026-11-03: a language model's rewrite of the original, under its terms