NAME
PDL::Slices -- Indexing, slicing, and dicing
SYNOPSIS
use PDL;
$x = ones(3,3);
$y = $x->slice('-1:0,(1)');
$c = $x->dummy(2);DESCRIPTION
This package provides many of the powerful PerlDL core index manipulation routines. These routines mostly allow two-way data flow, so you can modify your data in the most convenient representation. For example, you can make a 1000x1000 unit matrix with
$x = zeroes(1000,1000);
$x->diagonal(0,1) ++;which is quite efficient. See PDL::Indexing and PDL::Tips for more examples. As of 2.090, backward dataflow will be turned off if any input has inward-only dataflow, to avoid creating "loops". See PDL::Dataflow for more.
Slicing is so central to the PDL language that a special compile-time syntax has been introduced to handle it compactly; see PDL::NiceSlice for details.
PDL indexing and slicing functions usually include two-way data flow, so that you can separate the actions of reshaping your data structures and modifying the data themselves. Two special methods, "copy" and "sever", help you control the data flow connection between related variables.
$y = $x->slice("1:3"); # Slice maintains a link between $x and $y.
$y += 5; # $x is changed!If you want to force a physical copy and no data flow, you can copy or sever the slice expression:
$y = $x->slice("1:3")->copy;
$y += 5; # $x is not changed.
$y = $x->slice("1:3")->sever;
$y += 5; # $x is not changed.The difference between sever and copy is that sever acts on (and returns) its argument, while copy produces a disconnected copy. If you say
$y = $x->slice("1:3");
$c = $y->sever;then the variables $y and $c point to the same object but with ->copy they would not.
index, index1d, and index2d provide rudimentary index indirection.
$c = index($source,$ind);
$c = index1d($source,$ind);
$c = index2d($source2,$ind1,$ind2);use the $ind variables as indices to look up values in $source. The three routines broadcast slightly differently.
indexuses direct broadcasting for 1-D indexing across the 0 dim of$source. It can broadcast over source broadcast dims or index broadcast dims, but not (easily) both: If$sourcehas more than 1 dimension and$indhas more than 0 dimensions, they must agree in a broadcasting sense.index1duses a single active dim in$indto produce a list of indexed values in the 0 dim of the output - it is useful for collapsing$sourceby indexing with a single row of values along$source's 0 dimension. The output has the same number of dims as$source. The 0 dim of the output has size 1 if$indis a scalar, and the same size as the 0 dim of$indif it is not. If$indand$sourceboth have more than 1 dim, then all dims higher than 0 must agree in a broadcasting sense.index2dworks likeindexbut uses separate ndarrays for X and Y coordinates. For more general N-dimensional indexing, see the PDL::NiceSlice syntax or PDL::Slices (in particularslice,indexND, andrange).
These functions are two-way, i.e. after
$c = $x->index(pdl[0,5,8]);
$c .= pdl [0,2,4];the changes in $c will flow back to $x.
index provids simple broadcasting: multiple-dimensioned arrays are treated as collections of 1-D arrays, so that
$x = xvals(10,10)+10*yvals(10,10);
$y = $x->index(3);
$c = $x->index(9-xvals(10));puts a single column from $x into $y, and puts a single element from each column of $x into $c. If you want to extract multiple columns from an array in one operation, see "dice" or "indexND".
indexNDb
Backwards-compatibility alias for indexNDindexND
Find selected elements in an N-D ndarray, with optional boundary handling$out = $source->indexND( $index, [$method] )
$source = 10*xvals(10,10) + yvals(10,10);
$index = pdl([[2,3],[4,5]],[[6,7],[8,9]]);
print $source->indexND( $index );
[
[23 45]
[67 89]
]IndexND collapses $index by lookup into $source. The 0th dimension of $index is treated as coordinates in $source, and the return value has the same dimensions as the rest of $index. The returned elements are looked up from $source. Dataflow works -- propagated assignment flows back into $source.
IndexND and IndexNDb were originally separate routines but they are both now implemented as a call to "range", and have identical syntax to one another.
SEE ALSO:
"whichND" in PDL::Primitive returns N-D indices into a multidimensional PDL, suitable for feeding to this.
range
Extract selected chunks from a source ndarray, with boundary conditions
$out = $source->range($index,[$size,[$boundary]])Returns elements or rectangular slices of the original ndarray, indexed by the $index ndarray. $source is an N-dimensional ndarray, and $index is an ndarray whose first dimension has size up to N. Each row of $index is treated as coordinates of a single value or chunk from $source, specifying the location(s) to extract.
If you specify a single index location, then range is essentially an expensive slice, with controllable boundary conditions.
INPUTS
$index and $size can be ndarrays or array refs such as you would feed to zeroes and its ilk. If $index's 0th dimension has size higher than the number of dimensions in $source, then $source is treated as though it had trivial dummy dimensions of size 1, up to the required size to be indexed by $index -- so if your source array is 1-D and your index array is a list of 3-vectors, you get two dummy dimensions of size 1 on the end of your source array.
You can extract single elements or N-D rectangular ranges from $source, by setting $size. If $size is undef or zero, then you get a single sample for each row of $index. This behavior is similar to "indexNDb", which is in fact implemented as a call to "range".
If $size is positive then you get a range of values from $source at each location, and the output has extra dimensions allocated for them. $size can be a scalar, in which case it applies to all dimensions, or an N-vector, in which case each element is applied independently to the corresponding dimension in $source. See below for details.
$boundary is a number, string, or list ref indicating the type of boundary conditions to use when ranges reach the edge of $source. If you specify no boundary conditions the default is to forbid boundary violations on all axes. If you specify exactly one boundary condition, it applies to all axes. If you specify more (as elements of a list ref, or as a packed string, see below), then they apply to dimensions in the order in which they appear, and the last one applies to all subsequent dimensions. (This is less difficult than it sounds; see the examples below).
- 0 (synonyms: 'f','forbid') (default)
-
Ranges are not allowed to cross the boundary of the original PDL. Disallowed ranges throw an error. The errors are thrown at evaluation time, not at the time of the range call (this is the same behavior as "slice").
- 1 (synonyms: 't','truncate')
-
Values outside the original ndarray get BAD if you've got bad value support compiled into your PDL and set the badflag for the source PDL; or 0 if you haven't (you must set the badflag if you want BADs for out of bound values, otherwise you get 0). Reverse dataflow works OK for the portion of the child that is in-bounds. The out-of-bounds part of the child is reset to (BAD|0) during each dataflow operation, but execution continues.
- 2 (synonyms: 'e','x','extend')
-
Values that would be outside the original ndarray point instead to the nearest allowed value within the ndarray. See the CAVEAT below on mappings that are not single valued.
- 3 (synonyms: 'p','periodic')
-
Periodic boundary conditions apply: the numbers in $index are applied, strict-modulo the corresponding dimensions of $source. This is equivalent to duplicating the $source ndarray throughout N-D space. See the CAVEAT below about mappings that are not single valued.
- 4 (synonyms: 'm','mirror')
-
Mirror-reflection periodic boundary conditions apply. See the CAVEAT below about mappings that are not single valued.
The boundary condition identifiers all begin with unique characters, so you can feed in multiple boundary conditions as either a list ref or a packed string. (The packed string is marginally faster to run). For example, the four expressions [0,1], ['forbid','truncate'], ['f','t'], and 'ft' all specify that violating the boundary in the 0th dimension throws an error, and all other dimensions get truncated.
If you feed in a single string, it is interpreted as a packed boundary array if all of its characters are valid boundary specifiers (e.g. 'pet'), but as a single word-style specifier if they are not (e.g. 'forbid').
Where the source PDL is empty, all non-barfing boundary conditions are changed to truncation, since there is no data to reflect, extend, or mirror.
OUTPUT
The output broadcasts over both $index and $source. Because implicit broadcasting can happen in a couple of ways, a little thought is needed. The returned dimension list is stacked up like this:
(index broadcast dims), (index dims (size)), (source broadcast dims)The first few dims of the output correspond to the extra dims of $index (beyond the 0 dim). They allow you to pick out individual ranges from a large, broadcasted collection.
The middle few dims of the output correspond to the size dims specified in $size, and contain the range of values that is extracted at each location in $source. Every nonzero element of $size is copied to the dimension list here, so that if you feed in (for example) $size = [2,0,1] you get an index dim list of (2,1).
The last few dims of the output correspond to extra dims of $source beyond the number of dims indexed by $index. These dims act like ordinary broadcast dims, because adding more dims to $source just tacks extra dims on the end of the output. Each source broadcast dim ranges over the entire corresponding dim of $source.
Dataflow: Dataflow is bidirectional.
Examples: Here are basic examples of range operation, showing how to get ranges out of a small matrix. The first few examples show extraction and selection of individual chunks. The last example shows how to mark loci in the original matrix (using dataflow).
pdl> $src = 10*xvals(10,5)+yvals(10,5)
pdl> print $src->range([2,3]) # Cut out a single element
23
pdl> print $src->range([2,3],1) # Cut out a single 1x1 block
[
[23]
]
pdl> print $src->range([2,3], [2,1]) # Cut a 2x1 chunk
[
[23 33]
]
pdl> print $src->range([[2,3]],[2,1]) # Trivial list of 1 chunk
[
[
[23]
[33]
]
]
pdl> print $src->range([[2,3],[0,1]], [2,1]) # two 2x1 chunks
[
[
[23 1]
[33 11]
]
]
pdl> # A 2x2 collection of 2x1 chunks
pdl> print $src->range([[[1,1],[2,2]],[[2,3],[0,1]]],[2,1])
[
[
[
[11 22]
[23 1]
]
[
[21 32]
[33 11]
]
]
]
pdl> $src = xvals(5,3)*10+yvals(5,3)
pdl> print $src->range(3,1) # Broadcast over y dimension in $src
[
[30]
[31]
[32]
]
pdl> $src = zeroes(5,4);
pdl> $src->range(pdl([2,3],[0,1]),pdl(2,1)) .= xvals(2,2,1) + 1
pdl> print $src
[
[0 0 0 0 0]
[2 2 0 0 0]
[0 0 0 0 0]
[0 0 1 1 0]
]CAVEAT: It's quite possible to select multiple ranges that intersect. In that case, modifying the ranges doesn't have a guaranteed result in the original PDL -- the result is an arbitrary choice among the valid values. For some things that's OK; but for others it's not. In particular, this doesn't work:
pdl> $photon_list = PDL::RandVar->new->sample(500)->reshape(2,250)*10
pdl> $histogram = zeroes(10,10)
pdl> $histogram->range($photon_list,1)++; #not what you wantedThe reason is that if two photons land in the same bin, then that bin doesn't get incremented twice. (That may get fixed in a later version...)
PERMISSIVE RANGING: If $index has too many dimensions compared to $source, then $source is treated as though it had dummy dimensions of size 1, up to the required number of dimensions. These virtual dummy dimensions have the usual boundary conditions applied to them.
If the 0 dimension of $index is ludicrously large (if its size is more than 5 greater than the number of dims in the source PDL) then range will insist that you specify a size in every dimension, to make sure that you know what you're doing. That catches a common error with range usage: confusing the initial dim (which is usually small) with another index dim (perhaps of size 1000).
If the index variable is Empty, then range() always returns the Empty PDL. If the index variable is not Empty, indexing it always yields a boundary violation. All non-barfing conditions are treated as truncation, since there are no actual data to return.
EFFICIENCY: Because range isn't an affine transformation (it involves lookup into a list of N-D indices), it is somewhat memory-inefficient for long lists of ranges, and keeping dataflow open is much slower than for affine transformations (which don't have to copy data around).
Doing operations on small subfields of a large range is inefficient because the engine must flow the entire range back into the original PDL with every atomic perl operation, even if you only touch a single element. One way to speed up such code is to sever your range, so that PDL doesn't have to copy the data with each operation, then copy the elements explicitly at the end of your loop. Here's an example that labels each region in a range sequentially, using many small operations rather than a single xvals assignment:
### How to make a collection of small ops run fast with range...
$x = $data->range($index, $sizes, $bound)->sever;
$aa = $data->range($index, $sizes, $bound);
$x($_ - 1) .= $_ for 1..$x->nelem; # Lots of little ops
$aa .= $x;range is a perl front-end to a PP function, rangeb. Calling rangeb is marginally faster but requires that you include all arguments.
DEVEL NOTES
* index broadcast dimensions are effectively clumped internally. This makes it easier to loop over the index array but a little more brain-bending to tease out the algorithm.
Engine for "range"
Same calling convention as "range", but you must supply all parameters. rangeb is marginally faster as it makes a direct PP call, avoiding the perl argument-parsing step. EOD HandleBad => 1, TwoWay => 1, Lvalue => 1, P2Child => 1,
# rdim: dimensionality of each range (0 dim of index PDL) # nitems: total number of list elements (product of itdims) # itdim: number of index broadcast dimensions # ntsize: number of nonzero size dimensions # bsize: Number of independently specified boundary conditions # nsizes: Number of independently specified range dim sizes # sizes: array of range sizes, indexed (0..rdim-1). A zero element means # that the dimension is omitted from the child dim list. # itdims: Size of each index broadcast dimension, indexed (0..itdim-1) # corners: parent coordinates of each corner, running fastest over coord index. # (indexed 0 .. (nitems-1)*(rdim)+rdim-1) # boundary: Array containing all the boundary condition specs
Comp => 'PDL_Indx rdim;
PDL_Indx nitems;
PDL_Indx itdim;
PDL_Indx ntsize;
PDL_Indx bsize;
PDL_Indx nsizes;
PDL_Indx sizes[$COMP(rdim)];
PDL_Indx itdims[$COMP(itdim)];
PDL_Indx corners[$COMP(rdim) * $COMP(nitems)];
char boundary[$COMP(rdim)];
char size_pdl_destroy;
char ind_pdl_destroy;
',
MakeComp => pp_line_numbers(__LINE__, <<'EOD-MakeComp'),
pdl *size_pdl = NULL;
PDL_RETERROR(PDL_err, PDL->make_physdims(ind_pdl));
$COMP(size_pdl_destroy) = $COMP(ind_pdl_destroy) = 0;/* Generalized empties are ok, but not in the special 0 dim (the index vector) */ if(ind_pdl->dims[0] == 0) { $CROAK("can't handle Empty indices -- call range instead"); }
/*** * Ensure that the index is a PDL_Indx. If there's no loss of information, * just upgrade it -- otherwise, make a temporary copy. */ switch(ind_pdl->datatype) { default: /* Most types: */ ind_pdl = PDL->hard_copy(ind_pdl); /* copy and fall through */ if (!ind_pdl) $CROAK("Error in hard_copy"); $COMP(ind_pdl_destroy) = 1; case PDL_SB: case PDL_B: case PDL_S: case PDL_US: case PDL_L: case PDL_UL: case PDL_LL: case PDL_ULL: PDL_RETERROR(PDL_err, PDL->converttype(ind_pdl,PDL_IND)); /* convert in place. */ break; case PDL_IND: PDL_RETERROR(PDL_err, PDL->make_physical(ind_pdl)); break; }
/*** * Figure sizes of the COMP arrays and allocate them. */ { PDL_Indx i,nitems;
$COMP(rdim) = ind_pdl->ndims ? ind_pdl->dims[0] : 1;
for(i=nitems=1; i < ind_pdl->ndims; i++) /* Accumulate item list size */
nitems *= ind_pdl->dims[i];
$COMP(nitems) = nitems;
$COMP(itdim) = ind_pdl->ndims ? ind_pdl->ndims - 1 : 0;
$DOCOMPALLOC();
}/*** * Fill in the boundary condition array */ { char *bstr; STRLEN blen; bstr = SvPV(boundary_sv,blen);
if(blen == 0) {
/* If no boundary is specified then every dim gets forbidden */
PDL_Indx i;
for (i=0;i<$COMP(rdim);i++)
$COMP(boundary)[i] = 0;
} else {
PDL_Indx i;
for(i=0;i<$COMP(rdim);i++) {
switch(bstr[PDLMIN(i, blen-1)]) {
case '0': case 'f': case 'F': /* forbid */
$COMP(boundary)[i] = 0;
break;
case '1': case 't': case 'T': /* truncate */
$COMP(boundary)[i] = 1;
break;
case '2': case 'e': case 'E': /* extend */
$COMP(boundary)[i] = 2;
break;
case '3': case 'p': case 'P': /* periodic */
$COMP(boundary)[i] = 3;
break;
case '4': case 'm': case 'M': /* mirror */
$COMP(boundary)[i] = 4;
break;
default:
{
/* No need to check if i < blen -- this will barf out the
* first time it gets hit.
*/
$CROAK("Unknown boundary condition '%c' in range",bstr[i]);
}
break;
} // end of switch
}
}
}
/***
* Store the sizes of the index-broadcast dims
*/
{
PDL_Indx i;
PDL_Indx nd = ind_pdl->ndims - 1;
for(i=0; i < nd ; i++)
$COMP(itdims)[i] = ind_pdl->dims[i+1];
}/*** * Check and condition the size ndarray, and store sizes of the ranges */ { PDL_Indx i,ntsize;
if( (size_sv == NULL) || !SvOK(size_sv) ) {
// NO size was passed in, not normally executed even if you passed in no size to range(),
// as range() generates a size array...
for(i=0;i<$COMP(rdim);i++)
$COMP(sizes)[i] = 0;
} else {
/* Normal case with sizes present in a PDL */
if(!(size_pdl = PDL->SvPDLV(size_sv))) /* assignment */
$CROAK("Unable to convert size to a PDL");
if(size_pdl->nvals == 0) {
// no values in the size_pdl - Empty or Null. Just copy 0s to all the range dims
for(i=0;i<$COMP(rdim);i++)
$COMP(sizes)[i] = 0;
} else {
// Convert size PDL to PDL_IND to support indices
switch(size_pdl->datatype) {
default: /* Most types: */
size_pdl = PDL->hard_copy(size_pdl); /* copy and fall through */
if (!size_pdl) $CROAK("Error in hard_copy");
$COMP(size_pdl_destroy) = 1;
case PDL_SB: case PDL_B: case PDL_S: case PDL_US: case PDL_L: case PDL_UL: case PDL_LL: case PDL_ULL:
PDL_RETERROR(PDL_err, PDL->converttype(size_pdl,PDL_IND)); /* convert in place. */
break;
case PDL_IND:
PDL_RETERROR(PDL_err, PDL->make_physical(size_pdl));
break;
}
$COMP(nsizes) = size_pdl->nvals; /* Store for later permissiveness check */
/* Copy the sizes, or die if they're the wrong shape */
if(size_pdl->nvals == 1) {
for(i=0;i<$COMP(rdim);i++) {
$COMP(sizes)[i] = *((PDL_Indx *)(size_pdl->data));
}
/* Check for nonnegativity of sizes. The rdim>0 mask ensures that */
/* we don't barf on the Empty PDL (as an index). */
if( $COMP(rdim) > 0 && $COMP(sizes)[0] < 0 ) {
$CROAK("Negative range size is not allowed\n");
}
}
else if( size_pdl->nvals <= $COMP(rdim) && size_pdl->ndims == 1) {
for(i=0;i<$COMP(rdim);i++) {
$COMP(sizes)[i] = ( (i < size_pdl->nvals) ?
((PDL_Indx *)(size_pdl->data))[i] :
0
);
if($COMP(sizes)[i] < 0)
$CROAK("Negative range sizes are not allowed\n");
}
}
else {
$CROAK("Size must match index's 0th dim\n");
}
} /* end of nonempty size-ndarray code */
} /* end of defined-size-ndarray code */
/* Insert the number of nontrivial sizes (these get output dimensions) */
for(i=ntsize=0;i<$COMP(rdim);i++)
if($COMP(sizes)[i])
ntsize++;
$COMP(ntsize) = ntsize;
}/*** * Stash coordinates of the corners */
PDL_Indx j,k,ioff, iter[$COMP(itdim)]; /* initialize iterator to loop over index broadcasts */ PDL_Indx *cptr = iter; for(k=0;k<$COMP(itdim);k++) *(cptr++) = 0; cptr = $COMP(corners); do { /* accumulate offset into the index from the iterator */ for(k=ioff=0;k<$COMP(itdim);k++) ioff += iter[k] * ind_pdl->dimincs[k+1]; /* Loop over the 0th dim of index, copying coords. */ /* This is the natural place to check for permissive ranging; too */ /* bad we don't have access to the parent ndarray here... */ for(j=0;j<$COMP(rdim);j++) *(cptr++) = ((PDL_Indx *)(ind_pdl->data))[ioff + ind_pdl->dimincs[0] * j]; /* Increment the iterator -- the test increments, the body carries. */ for(k=0; k<$COMP(itdim) && (++(iter[k]))>=($COMP(itdims)[k]) ;k++) iter[k] = 0; } while(k<$COMP(itdim)); if ($COMP(ind_pdl_destroy)) PDL->destroy(ind_pdl); /* finished with our copy */ if ($COMP(size_pdl_destroy)) PDL->destroy(size_pdl); /* finished with our copy */ EOD-MakeComp
RedoDims => pp_line_numbers(__LINE__, <<'EOD-RedoDims'),
PDL_Indx stdim = $PDL(PARENT)->ndims - $COMP(rdim);
PDL_Indx dim,inc;
PDL_Indx i,rdvalid;
// Speed bump for ludicrous cases
if( $COMP(rdim) > $PDL(PARENT)->ndims+5 && $COMP(nsizes) != $COMP(rdim)) {
$CROAK(
"Ludicrous number of extra dims in range index; leaving child null.\n"
" (%"IND_FLAG" implicit dims is > 5; index has %"IND_FLAG" dims; source has %"IND_FLAG" dim%s.)\n"
" This often means that your index PDL is incorrect.\n"
" To avoid this message, allocate dummy dims in\n"
" the source or use %"IND_FLAG" dims in range's size field.\n",
$COMP(rdim)-$PDL(PARENT)->ndims,$COMP(rdim),$PDL(PARENT)->ndims,
$PDL(PARENT)->ndims>1?"s":"",$COMP(rdim)
);
}
if(stdim < 0)
stdim = 0;
/* Set dimensionality of child */
$SETNDIMS($COMP(itdim)+$COMP(ntsize)+stdim);
inc = 1;
/* Copy size dimensions to child, crunching as we go. */
dim = $COMP(itdim);
for(i=rdvalid=0;i<$COMP(rdim);i++) {
if($COMP(sizes)[i]) {
rdvalid++;
$PDL(CHILD)->dimincs[dim] = inc;
inc *= ($PDL(CHILD)->dims[dim++] = $COMP(sizes)[i]); /* assignment */
}
}
/* Copy index broadcast dimensions to child */
for(dim=0; dim<$COMP(itdim); dim++) {
$PDL(CHILD)->dimincs[dim] = inc;
inc *= ($PDL(CHILD)->dims[dim] = $COMP(itdims)[dim]); /* assignment */
}
/* Copy source broadcast dimensions to child */
dim = $COMP(itdim) + rdvalid;
for(i=0;i<stdim;i++) {
$PDL(CHILD)->dimincs[dim] = inc;
inc *= ($PDL(CHILD)->dims[dim++] = $PDL(PARENT)->dims[i+$COMP(rdim)]); /* assignment */
}
/* Cover bizarre case where the source PDL is empty - in that case, change */
/* all non-barfing boundary conditions to truncation, since we have no data */
/* to reflect, extend, or mirror. */
if($PDL(PARENT)->dims[0] == 0) {
for(dim=0; dim<$COMP(rdim); dim++) {
if($COMP(boundary)[dim])
$COMP(boundary)[dim] = 1; // force truncation
}
}$PDL(CHILD)->datatype = $PDL(PARENT)->datatype; $SETDIMS(); EOD-RedoDims
EquivCPOffsCode => pp_line_numbers(__LINE__, <<'EOD-EquivCPOffsCode'),
PDL_Indx *ip; /* vector iterator */
PDL_Indx *sp; /* size vector including stdims */
PDL_Indx *coords; /* current coordinates */PDL_Indx k; /* index */ PDL_Indx item; /* index broadcast iterator */ PDL_Indx pdim = $PDL(PARENT)->ndims; PDL_Indx rdim = $COMP(rdim); PDL_Indx prdim = PDLMIN(rdim,pdim); PDL_Indx iter2[pdim * 2 + rdim]; PDL_Indx *sizes = iter2 + pdim; coords = sizes + pdim;
/* Figure out size vector */ for(ip = $COMP(sizes), sp = sizes, k=0; k<rdim; k++) *(sp++) = *(ip++); for(; k < pdim; k++) *(sp++) = $PDL(PARENT)->dims[k];
/* Loop over all the ranges in the index list */ for(item=0; item<$COMP(nitems); item++) {
/* initialize in-range iterator to loop within each range */
for(ip = iter2, k=0; k<pdim; k++)
*(ip++) = 0;
do {
PDL_Indx poff = 0;
PDL_Indx coff;
PDL_Indx k2;
char trunc = 0; /* Flag used to skip truncation case */
/* Collect and boundary-check the current N-D coords */
for(k=0; k < prdim; k++){
PDL_Indx ck = iter2[k] + $COMP(corners)[ item * rdim + k ] ;
/* normal case */
if(ck < 0 || ck >= $PDL(PARENT)->dims[k]) {
switch($COMP(boundary)[k]) {
case 0: /* no boundary breakage allowed */
$CROAK("index out-of-bounds in range (index vector #%ld)",item);
break;
case 1: /* truncation */
trunc = 1;
break;
case 2: /* extension -- crop */
ck = (ck >= $PDL(PARENT)->dims[k]) ? $PDL(PARENT)->dims[k]-1 : 0;
break;
case 3: /* periodic -- mod it */
ck %= $PDL(PARENT)->dims[k];
if(ck < 0) /* Fix mod breakage in C */
ck += $PDL(PARENT)->dims[k];
break;
case 4: /* mirror -- reflect off the edges */
ck += $PDL(PARENT)->dims[k];
ck %= ($PDL(PARENT)->dims[k] * 2);
if(ck < 0) /* Fix mod breakage in C */
ck += $PDL(PARENT)->dims[k]*2;
ck -= $PDL(PARENT)->dims[k];
if(ck < 0) {
ck *= -1;
ck -= 1;
}
break;
default:
$CROAK("Unknown boundary condition in range -- bug alert!");
break;
}
}
coords[k] = ck;
}
/* Check extra dimensions -- pick up where k left off... */
for( ; k < rdim ; k++) {
/* Check for indexing off the end of the dimension list */
PDL_Indx ck = iter2[k] + $COMP(corners)[ item * rdim + k ] ;
switch($COMP(boundary)[k]) {
case 0: /* No boundary breakage allowed -- nonzero corners cause barfage */
if(ck != 0)
$CROAK("Too many dims in range index (and you've forbidden boundary violations)");
break;
case 1: /* truncation - just truncate if the corner is nonzero */
trunc |= (ck != 0);
break;
case 2: /* extension -- ignore the corner (same as 3) */
case 3: /* periodic -- ignore the corner */
case 4: /* mirror -- ignore the corner */
ck = 0;
break;
default:
$CROAK("Unknown boundary condition in range -- bug alert!");
break;
}
}
/* Find offsets into the child and parent arrays, from the N-D coords */
/* Note we only loop over real source dims (prdim) to accumulate -- */
/* because the offset is trivial and/or we're truncating for virtual */
/* dims caused by permissive ranging. */
coff = $PDL(CHILD)->dimincs[0] * item;
for(k2 = $COMP(itdim), poff = k = 0;
k < prdim;
k++) {
poff += coords[k]*$PDL(PARENT)->dimincs[k];
if($COMP(sizes)[k])
coff += iter2[k] * $PDL(CHILD)->dimincs[k2++];
}
/* Loop the copy over all the source broadcast dims (above rdim). */
do {
PDL_Indx poff1 = poff;
PDL_Indx coff1 = coff;
/* Accumulate the offset due to source broadcasting */
for(k2 = $COMP(itdim) + $COMP(ntsize), k = rdim;
k < pdim;
k++) {
poff1 += iter2[k] * $PDL(PARENT)->dimincs[k];
coff1 += iter2[k] * $PDL(CHILD)->dimincs[k2++];
}
/* Finally -- make the copy
* EQUIVCPTRUNC works like EQUIVCPOFFS but with checking for
* out-of-bounds conditions.
*/
$EQUIVCPTRUNC(coff1,poff1,trunc);
/* Increment the source broadcast iterator */
for( k=$COMP(rdim);
k < pdim && (++(iter2[k]) >= $PDL(PARENT)->dims[k]);
k++)
iter2[k] = 0;
} while(k < pdim); /* end of source-broadcast iteration */
/* Increment the in-range iterator */
for(k = 0;
k < $COMP(rdim) && (++(iter2[k]) >= $COMP(sizes)[k]);
k++)
iter2[k] = 0;
} while(k < $COMP(rdim)); /* end of main iteration */
} /* end of item do loop */
EOD-EquivCPOffsCode
);pp_def('rld', GenericTypes => [ppdefs_all], Pars=>'indx a(n); b(n); [o]c(m);', OtherPars=>'IV sumover_max => m', PMCode =>pp_line_numbers(__LINE__, <<'EOD'), sub PDL::rld { my ($x,$y) = @_; my ($c,$sm) = @_ == 3 ? ($_[2], $_[2]->dim(0)) : (PDL->null, $x->sumover->max->sclr); PDL::_rld_int($x,$y,$c,$sm); $c; } EOD Code => pp_line_numbers(__LINE__, <<'EOF'), PDL_Indx j=0; loop (n) %{ PDL_Indx jlim = j + $a(); $GENERIC(b) bv = $b(); loop (m=j:jlim) %{ $c() = bv; %} j = jlim; %} EOF Doc => <<'EOD' =for ref
Run-length decode a vector
Given a vector $x of the numbers of instances of values $y, run-length decode to $c. EOD );
pp_def('rle', GenericTypes => [ppdefs_all], Pars=>'c(n); indx [o]a(m=CALC($SIZE(n))); [o]b(m);', PMCode =>pp_line_numbers(__LINE__, <<'EOC'), sub PDL::rle { my $c = shift; my ($x,$y) = @_==2 ? @_ : (null,null); PDL::_rle_int($c,$x,$y); my $max_ind = ($c->ndims<2) ? ($x!=0)->sumover-1 : ($x!=0)->clump(1..$x->ndims-1)->sumover->max->sclr-1; return ($x->slice("0:$max_ind"),$y->slice("0:$max_ind")); } EOC Code =>pp_line_numbers(__LINE__, <<'EOF'), PDL_Indx j=0; $GENERIC(c) cv, clv = $c(n=>0); $b(m=>0) = clv; $a(m=>0) = 0; loop (n) %{ cv = $c(); if (cv == clv) { $a(m=>j)++; } else { j++; $b(m=>j) = clv = cv; $a(m=>j) = 1; } %} loop (m=j+1) %{ $a() = 0; $b() = 0; %} EOF Doc => <<'EOD' =for ref
Run-length encode a vector
Given vector $c, generate a vector $x with the number of each element, and a vector $y of the unique values. New in PDL 2.017, only the elements up to the first instance of 0 in $x are returned, which makes the common use case of a 1-dimensional $c simpler. For broadcast operation, $x and $y will be large enough to hold the largest row of $y, and only the elements up to the first instance of 0 in each row of $x should be considered.
$c = floor(4*random(10));
rle($c,$x=null,$y=null);
#or
($x,$y) = rle($c);
#for $c of shape [10, 4]:
$c = floor(4*random(10,4));
($x,$y) = rle($c);
#to see the results of each row one at a time:
foreach (0..$c->dim(1)-1){
my ($as,$bs) = ($x(:,($_)),$y(:,($_)));
my ($ta,$tb) = where($as,$bs,$as!=0); #only the non-zero elements of $x
print $c(:,($_)) . " rle==> " , ($ta,$tb) , "\trld==> " . rld($ta,$tb) . "\n";
}
# the inverse of (chance of all 6 3d6 rolls being >= each possible sum)
($nrolls, $ndice, $dmax) = (6, 3, 6);
($x, $x1) = (allaxisvals(($dmax) x $ndice)+1)->sumover->flat->qsort->rle;
$y = $x->cumusumover;
$yprob1x = $y->slice('-1:0')->double / $y->slice('(-1)');
$z = cat($x1, 1 / $yprob1x**$nrolls)->transpose;Run-length encode a set of vectors.
Higher-order rle(), for use with qsortvec().
Given set of vectors $c, generate a vector $a with the number of occurrences of each element (where an "element" is a vector of length $M occurring in $c), and a set of vectors $b containing the unique values. As for rle(), only the elements up to the first instance of 0 in $a should be considered.
Can be used together with clump() to run-length encode "values" of arbitrary dimensions. Can be used together with rotate(), cat(), append(), and qsortvec() to count N-grams over a 1d PDL.
See also: "rle", "qsortvec" in PDL::Ufunc, "uniqvec" in PDL::Primitive Contributed by Bryan Jurish <moocow@cpan.org>.
EOD );
pp_def('rldvec', Pars => 'indx a(uniqvals); b(M,uniqvals); [o]c(M,decodedvals)', OtherPars=>'IV sumover_max => decodedvals', PMCode =>pp_line_numbers(__LINE__, <<'EOC'), sub PDL::rldvec { my ($a,$b,$c) = @_; ($c,my $sm) = defined($c) ? ($c,$c->dim(1)) : (PDL->null,$a->sumover->max->sclr); PDL::_rldvec_int($a,$b,$c,$sm); return $c; } EOC Code =>pp_line_numbers(__LINE__, <<'EOC'), PDL_Indx cn=0; loop (uniqvals) %{ PDL_Indx cnlim = cn + $a(); loop (decodedvals=cn:cnlim,M) %{ $c() = $b(); %} cn = cnlim; %} EOC Doc =><<'EOD' =for ref
Run-length decode a set of vectors, akin to a higher-order rld().
Given a vector $a() of the number of occurrences of each row, and a set $b() of row-vectors each of length $M, run-length decode to $c().
Can be used together with clump() to run-length decode "values" of arbitrary dimensions.
See also: "rld". Contributed by Bryan Jurish <moocow@cpan.org>. EOD );
pp_def('rleseq', Pars => "c(N); indx [o]a(N); [o]b(N)", Code =>pp_line_numbers(__LINE__, <<'EOC'), PDL_Indx j=0; $GENERIC(c) coff; coff = $c(N=>0); $b(N=>0) = coff; $a(N=>0) = 0; loop (N) %{ if ($c() == coff+$a(N=>j)) { $a(N=>j)++; } else { j++; $b(N=>j) = coff = $c(); $a(N=>j) = 1; } %} loop (N=j+1) %{ $a() = 0; $b() = 0; %} EOC Doc =><<'EOD', =for ref
Run-length encode a vector of subsequences.
Given a vector of $c() of concatenated variable-length, variable-offset subsequences, generate a vector $a containing the length of each subsequence and a vector $b containing the subsequence offsets. As for rle(), only the elements up to the first instance of 0 in $a should be considered.
See also "rle". Contributed by Bryan Jurish <moocow@cpan.org>. EOD );
pp_def('rldseq', Pars => 'indx a(N); b(N); [o]c(M)', OtherPars=>'IV sumover_max => M', PMCode =>pp_line_numbers(__LINE__, <<'EOC'), sub PDL::rldseq { my ($a,$b,$c) = @_; ($c,my $sm) = defined($c) ? ($c,$c->dim(1)) : (PDL->null,$a->sumover->max->sclr); PDL::_rldseq_int($a,$b,$c,$sm); return $c; } EOC Code =>pp_line_numbers(__LINE__, <<'EOC'), PDL_Indx mi=0; loop (N) %{ PDL_Indx milim = mi + $a(), li = 0; loop (M=mi:milim) %{ $c() = $b() + li++; %} mi = milim; %} EOC Doc =><<'EOD' =for ref
Run-length decode a subsequence vector.
Given a vector $a() of sequence lengths and a vector $b() of corresponding offsets, decode concatenation of subsequences to $c(), as for:
$c = null;
$c = $c->append($b($_)+sequence($a->type,$a($_))) foreach (0..($N-1));See also: "rld". Contributed by Bryan Jurish <moocow@cpan.org>. EOD );
pp_add_exported('','rleND rldND'); pp_addpm(<<'EOF'); =head2 rleND
Signature: (data(@vdims,N); int [o]counts(N); [o]elts(@vdims,N))Run-length encode a set of (sorted) n-dimensional values.
Generalization of rle() and vv_rlevec(): given set of values $data, generate a vector $counts with the number of occurrences of each element (where an "element" is a matrix of dimensions @vdims occurring as a sequential run over the final dimension in $data), and a set of vectors $elts containing the elements which begin a run. Really just a wrapper for clump() and rlevec().
See also: "rle", "rlevec". Contributed by Bryan Jurish <moocow@cpan.org>.
rldND
Signature: (int counts(N); elts(@vdims,N); [o]data(@vdims,N);)Run-length decode a set of (sorted) n-dimensional values.
Generalization of rld() and rldvec(): given a vector $counts() of the number of occurrences of each @vdims-dimensioned element, and a set $elts() of @vdims-dimensioned elements, run-length decode to $data().
Really just a wrapper for clump() and rldvec().
See also: "rld", "rldvec". Contributed by Bryan Jurish <moocow@cpan.org>.
exchange two dimensions
Negative dimension indices count from the end.
The command
$y = $x->xchg(2,3);creates $y to be like $x except that the dimensions 2 and 3 are exchanged with each other i.e.
$y->at(5,3,2,8) == $x->at(5,3,8,2)reorder
Re-orders the dimensions of a PDL based on the supplied list.
Similar to the "xchg" method, this method re-orders the dimensions of a PDL. While the "xchg" method swaps the position of two dimensions, the reorder method can change the positions of many dimensions at once.
# Completely reverse the dimension order of a 6-Dim array.
$reOrderedPDL = $pdl->reorder(5,4,3,2,1,0);The argument to reorder is an array representing where the current dimensions should go in the new array. In the above usage, the argument to reorder (5,4,3,2,1,0) indicates that the old dimensions ($pdl's dims) should be re-arranged to make the new pdl ($reOrderPDL) according to the following:
Old Position New Position
------------ ------------
5 0
4 1
3 2
2 3
1 4
0 5You do not need to specify all dimensions, only a complete set starting at position 0. (Extra dimensions are left where they are). This means, for example, that you can reorder() the X and Y dimensions of an image, and not care whether it is an RGB image with a third dimension running across color plane.
Example:
pdl> $x = sequence(5,3,2); # Create a 3-d Array
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]
]
]
pdl> p $x->reorder(2,1,0); # Reverse the order of the 3-D PDL
[
[
[ 0 15]
[ 5 20]
[10 25]
]
[
[ 1 16]
[ 6 21]
[11 26]
]
[
[ 2 17]
[ 7 22]
[12 27]
]
[
[ 3 18]
[ 8 23]
[13 28]
]
[
[ 4 19]
[ 9 24]
[14 29]
]
]The above is a simple example that could be duplicated by calling $x->xchg(0,2), but it demonstrates the basic functionality of reorder.
As this is an index function, any modifications to the result PDL will change the parent.
move a dimension to another position
The command
$y = $x->mv(4,1);creates $y to be like $x except that the dimension 4 is moved to the place 1, so:
$y->at(1,2,3,4,5,6) == $x->at(1,5,2,3,4,6);The other dimensions are moved accordingly. Negative dimension indices count from the end.
using
Returns list of columns requested
line $pdl->using(1,2);Plot, as a line, column 1 of $pdl vs. column 2
pdl> $pdl = rcols("file");
pdl> line $pdl->using(1,2);meshgrid
Returns list of given 1-D vectors, but each expanded to match dims using "dummy" in PDL::Core.
meshgrid($vec1, $vec2, $vec3);print map $_->info, meshgrid(xvals(3), xvals(4), xvals(2));
# PDL: Double D [3,4,2] PDL: Double D [3,4,2] PDL: Double D [3,4,2]Returns an ndarray of lags to parent.
Usage:
I.e. if $x contains
[0,1,2,3,4,5,6,7]then
$y = $x->lags(0,2,2);is a (5,2) matrix
[2,3,4,5,6,7]
[0,1,2,3,4,5]This order of returned indices is kept because the function is called "lags" i.e. the nth lag is n steps behind the original.
$step and $nlags must be positive. $nthdim can be negative and will then be counted from the last dim backwards in the usual way (-1 = last dim).
Splits a dimension in the parent ndarray (opposite of clump). As of 2.076, throws exception if non-divisible nsp given, and can give negative nthdim which then counts backwards.
After
$y = $x->splitdim(2,3);the expression
$y->at(6,4,m,n,3,6) == $x->at(6,4,m+3*n)is always true (m has to be less than 3).
internal
Put some dimensions to a broadcastid.
$y = $x->broadcastI(0,1,5); # broadcast over dims 1,5 in id 1
EOD
P2Child => 1,
TwoWay => 1,
Lvalue => 1,
AffinePriv => 1,
OtherPars => "PDL_Indx id; PDL_Indx whichdims[]",
Comp => 'PDL_Indx nrealwhichdims',
MakeComp => pp_line_numbers(__LINE__, <<'EOF'),
PDL_Indx i,j;
$COMP(nrealwhichdims) = 0;
for(i=0; i<$COMP(whichdims_count); i++) {
for(j=i+1; j<$COMP(whichdims_count); j++)
if($COMP(whichdims)[i] == $COMP(whichdims)[j] &&
$COMP(whichdims)[i] != -1) {
$CROAK("duplicate arg %"IND_FLAG" %"IND_FLAG" %"IND_FLAG"",
i,j,$COMP(whichdims)[i]);
}
if($COMP(whichdims)[i] != -1) {
$COMP(nrealwhichdims) ++;
}
}
EOF
RedoDims => pp_line_numbers(__LINE__, <<'EOF'),
PDL_Indx nthc,i,j,flag;
$SETNDIMS($PDL(PARENT)->ndims);
$DOPRIVALLOC();
$PRIV(offs) = 0;
nthc=0;
for(i=0; i<$PDL(PARENT)->ndims; i++) {
flag=0;
if($PDL(PARENT)->nbroadcastids > $COMP(id) && $COMP(id) >= 0 &&
i == $PDL(PARENT)->broadcastids[$COMP(id])) {
nthc += $COMP(whichdims_count);
}
for(j=0; j<$COMP(whichdims_count); j++) {
if($COMP(whichdims)[j] == i) {flag=1; break;}
}
if(flag) {
continue;
}
$PDL(CHILD)->dims[nthc] = $PDL(PARENT)->dims[i];
$PRIV(incs)[nthc] = $PDL(PARENT)->dimincs[i];
nthc++;
}
for(i=0; i<$COMP(whichdims_count); i++) {
PDL_Indx cdim,pdim;
cdim = i +
($PDL(PARENT)->nbroadcastids > $COMP(id) && $COMP(id) >= 0?
$PDL(PARENT)->broadcastids[$COMP(id)] : $PDL(PARENT)->ndims)
- $COMP(nrealwhichdims);
pdim = $COMP(whichdims)[i];
if(pdim == -1) {
$PDL(CHILD)->dims[cdim] = 1;
$PRIV(incs)[cdim] = 0;
} else {
$PDL(CHILD)->dims[cdim] = $PDL(PARENT)->dims[pdim];
$PRIV(incs)[cdim] = $PDL(PARENT)->dimincs[pdim];
}
}
$SETDIMS();
PDL_RETERROR(PDL_err, PDL->reallocbroadcastids($PDL(CHILD),
PDLMAX($COMP(id)+1, $PDL(PARENT)->nbroadcastids)));
for(i=0; i<$PDL(CHILD)->nbroadcastids-1; i++) {
$PDL(CHILD)->broadcastids[i] =
($PDL(PARENT)->nbroadcastids > i ?
$PDL(PARENT)->broadcastids[i] : $PDL(PARENT)->ndims) +
(i <= $COMP(id) ? - $COMP(nrealwhichdims) :
$COMP(whichdims_count) - $COMP(nrealwhichdims));
}
$PDL(CHILD)->broadcastids[$PDL(CHILD)->nbroadcastids-1] = $PDL(CHILD)->ndims;
EOF
);pp_def('unbroadcast', Doc => <<'EOD', =for ref
All broadcasted dimensions are made real again.
See [TBD Doc] for details and examples.
dice
Dice rows/columns/planes out of a PDL using indexes for each dimension.
This function can be used to extract irregular subsets along many dimension of a PDL, e.g. only certain rows in an image, or planes in a cube. This can of course be done with the usual dimension tricks but this saves having to figure it out each time!
This method is similar in functionality to the "slice" method, but "slice" requires that contiguous ranges or ranges with constant offset be extracted. ( i.e. "slice" requires ranges of the form 1,2,3,4,5 or 2,4,6,8,10). Because of this restriction, "slice" is more memory efficient and slightly faster than dice
$slice = $data->dice([0,2,6],[2,1,6]); # Dicing a 2-D arrayThe arguments to dice are arrays (or 1D PDLs) for each dimension in the PDL. These arrays are used as indexes to which rows/columns/cubes,etc to dice-out (or extract) from the $data PDL.
Use X to select all indices along a given dimension (compare also mslice). As usual (in slicing methods) trailing dimensions can be omitted implying X'es for those.
pdl> $x = sequence(10,4)
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]
]
pdl> p $x->dice([1,2],[0,3]) # Select columns 1,2 and rows 0,3
[
[ 1 2]
[31 32]
]
pdl> p $x->dice(X,[0,3])
[
[ 0 1 2 3 4 5 6 7 8 9]
[30 31 32 33 34 35 36 37 38 39]
]
pdl> p $x->dice([0,2,5])
[
[ 0 2 5]
[10 12 15]
[20 22 25]
[30 32 35]
]As this is an index function, any modifications to the slice will change the parent (use the .= operator).
dice_axis
Dice rows/columns/planes from a single PDL axis (dimension) using index along a specified axis
This function can be used to extract irregular subsets along any dimension, e.g. only certain rows in an image, or planes in a cube. This can of course be done with the usual dimension tricks but this saves having to figure it out each time!
$slice = $data->dice_axis($axis,$index);pdl> $x = sequence(10,4)
pdl> $idx = pdl(1,2)
pdl> p $x->dice_axis(0,$idx) # Select columns
[
[ 1 2]
[11 12]
[21 22]
[31 32]
]
pdl> $t = $x->dice_axis(1,$idx) # Select rows
pdl> $t.=0
pdl> p $x
[
[ 0 1 2 3 4 5 6 7 8 9]
[ 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0]
[30 31 32 33 34 35 36 37 38 39]
]The trick to using this is that the index selects elements along the dimensions specified, so if you have a 2D image axis=0 will select certain X values - i.e. extract columns
As this is an index function, any modifications to the slice will change the parent.
$slice = $data->slice([2,3],'x',[2,2,0],"-1:1:-1", "*3");Extract rectangular slices of an ndarray, from a string specifier, an array ref specifier, or a combination.
slice is the main method for extracting regions of PDLs and manipulating their dimensionality. You can call it directly or via the NiceSlice source prefilter that extends Perl syntax to include array slicing.
slice can extract regions along each dimension of a source PDL, subsample or reverse those regions, dice each dimension by selecting a list of locations along it, or basic PDL indexing routine. The selected subfield remains connected to the original PDL via dataflow. In most cases this neither allocates more memory nor slows down subsequent operations on either of the two connected PDLs.
You pass in a list of arguments. Each term in the list controls the disposition of one axis of the source PDL and/or returned PDL. Each term can be a string-format cut specifier, a list ref that gives the same information without recourse to string manipulation, or a PDL with up to 1 dimension giving indices along that axis that should be selected.
If you want to pass in a single string specifier for the entire operation, you can pass in a comma-delimited list as the first argument. slice detects this condition and splits the string into a regular argument list. This calling style is fully backwards compatible with slice calls from before PDL 2.006.
STRING SYNTAX
If a particular argument to slice is a string, it is parsed as a selection, an affine slice, or a dummy dimension depending on the form. Leading or trailing whitespace in any part of each specifier is ignored (though it is not ignored within numbers).
'',:, orX-- keep-
The empty string,
:, orXcause the entire corresponding dimension to be kept unchanged. <n>-- selection-
A single number alone causes a single index to be selected from the corresponding dimension. The dimension is kept (and reduced to size 1) in the output.
(<n>)-- selection and collapse-
A single number in parenthesis causes a single index to be selected from the corresponding dimension. The dimension is discarded (completely eliminated) in the output.
<n>:<m>-- select an inclusive range-
Two numbers separated by a colon selects a range of values from the corresponding axis, e.g.
3:4selects elements 3 and 4 along the corresponding axis, and reduces that axis to size 2 in the output. Both numbers are regularized so that you can address the last element of the axis with an index of-1. If, after regularization, the two numbers are the same, then exactly one element gets selected (just like the<n>case). If, after regulariation, the second number is lower than the first, then the resulting slice counts down rather than up -- e.g.-1:0will return the entire axis, in reversed order. <n>:<m>:<s>-- select a range with explicit step-
If you include a third parameter, it is the stride of the extracted range. For example,
0:-1:2will sample every other element across the complete dimension. Specifying a stride of 1 prevents autoreversal -- so to ensure that your slice is *always* forward you can specify, e.g.,2:$n:1. In that case, an "impossible" slice gets an Empty PDL (with 0 elements along the corresponding dimension), so you can generate an Empty PDL with a slice of the form2:1:1. *<n>-- insert a dummy dimension-
Dummy dimensions aren't present in the original source and are "mocked up" to match dimensional slots, by repeating the data in the original PDL some number of times. An asterisk followed by a number produces a dummy dimension in the output, for example
*2will generate a dimension of size 2 at the corresponding location in the output dim list. Omitting the number (and using just an asterisk) inserts a dummy dimension of size 1.
ARRAY REF SYNTAX
If you feed in an ARRAY ref as a slice term, then it can have 0-3 elements. The first element is the start of the slice along the corresponding dim; the second is the end; and the third is the stepsize. Different combinations of inputs give the same flexibility as the string syntax.
[]- keep dim intact-
An empty ARRAY ref keeps the entire corresponding dim
[ 'X' ]- keep dim intact[ '*',$n ]- generate a dummy dim of size $n-
If $n is missing, you get a dummy dim of size 1.
[ $dex, 0, 0 ]- collapse and discard dim-
$dexmust be a single value. It is used to index the source, and the corresponding dimension is discarded. [ $start, $end ]- collect inclusive slice-
In the simple two-number case, you get a slice that runs up or down (as appropriate) to connect $start and $end.
[ $start, $end, $inc ]- collect inclusive slice-
The three-number case works exactly like the three-number string case above.
PDL args for dicing
If you pass in a 0- or 1-D PDL as a slicing argument, the corresponding dimension is "diced" -- you get one position along the corresponding dim, per element of the indexing PDL, e.g. $x->slice( pdl(3,4,9)) gives you elements 3, 4, and 9 along the 0 dim of $x.
Because dicing is not an affine transformation, it is slower than direct slicing even though the syntax is convenient.
$x->slice('1:3'); # return the second to fourth elements of $x
$x->slice('3:1'); # reverse the above
$x->slice('-2:1'); # return last-but-one to second elements of $x
$x->slice([1,3]); # Same as above three calls, but using array ref syntax
$x->slice([3,1]);
$x->slice([-2,1]);
EOD-slice
PMCode => pp_line_numbers(__LINE__, <<'EOD-slice'),
sub PDL::slice :lvalue {
my ($source, @others) = @_;
for my $i(0..$#others) {
my $idx = $others[$i];
if (ref $idx eq 'ARRAY') {
my @arr = map UNIVERSAL::isa($_, 'PDL') ? $_->flat->at(0) : $_, @{$others[$i]};
$others[$i] = \@arr;
next;
}
next if !( blessed($idx) && $idx->isa('PDL') );
# Deal with dicing. This is lame and slow compared to the
# faster slicing, but works okay. We loop over each argument,
# and if it's a PDL we dispatch it in the most straightforward
# way. Single-element and zero-element PDLs are trivial and get
# converted into slices for faster handling later.
barf("slice: dicing parameters must be at most 1D (arg $i)\n")
if $idx->ndims > 1;
my $nlm = $idx->nelem;
if($nlm > 1) {
#### More than one element - we have to dice (darn it).
$source = $source->mv($i,0)->index1d($idx)->mv(0,$i);
$others[$i] = '';
}
elsif($nlm) {
#### One element - convert to a regular slice.
$others[$i] = $idx->flat->at(0);
}
else {
#### Zero elements -- force an extended empty.
$others[$i] = "1:0:1";
}
}
PDL::_slice_int($source,my $o=$source->initialize,\@others);
$o;
}
EOD-slice
P2Child => 1,
OtherPars => 'pdl_slice_args *arglist;',
#
# Comp stash definitions:
# nargs - number of args in original call
# odim[] - maps argno to output dim (or -1 for squished dims)
# idim[] - maps argno to input dim (or -1 for dummy dims)
# odim_top - one more than the highest odim encountered
# idim_top - one more than the highest idim encountered
# start[] - maps argno to start index of slice range (inclusive)
# inc[] - maps argno to increment of slice range
# end[] - maps argno to end index of slice range (inclusive)
#
Comp => 'PDL_Indx nargs;
PDL_Indx odim[$COMP(nargs)];
PDL_Indx idim[$COMP(nargs)];
PDL_Indx idim_top;
PDL_Indx odim_top;
PDL_Indx start[$COMP(nargs)];
PDL_Indx inc[$COMP(nargs)];
PDL_Indx end[$COMP(nargs)];
',
AffinePriv => 1,
TwoWay => 1,
MakeComp => pp_line_numbers(__LINE__, <<'SLICE-MC'),
PDL_Indx nargs = 0;
pdl_slice_args *argsptr = arglist;
while (argsptr) nargs++, argsptr = argsptr->next;
$COMP(nargs) = nargs;
$DOCOMPALLOC();
PDL_Indx i, idim, odim;
argsptr = arglist;
for(odim=idim=i=0; i<nargs; i++) {
/* Copy parsed values into the limits */
$COMP(start)[i] = argsptr->start;
$COMP(end)[i] = argsptr->end;
$COMP(inc)[i] = argsptr->inc;
/* Deal with dimensions */
$COMP(odim)[i] = argsptr->squish ? -1 : odim++;
$COMP(idim)[i] = argsptr->dummy ? -1 : idim++;
argsptr = argsptr->next;
} /* end of arg-parsing loop */
$COMP(idim_top) = idim;
$COMP(odim_top) = odim;
SLICE-MC
RedoDims => pp_line_numbers(__LINE__, <<'EOF'),
PDL_Indx i;
PDL_Indx o_ndims_extra = PDLMAX(0, $PDL(PARENT)->ndims - $COMP(idim_top));
/* slurped dims from the arg parsing, plus any extra broadcast dims */
$SETNDIMS( $COMP(odim_top) + o_ndims_extra );
$DOPRIVALLOC();
$PRIV(offs) = 0; /* Offset vector to start of slice */
for(i=0; i<$COMP(nargs); i++) {
/** Belt-and-suspenders **/
if (($COMP(idim)[i] < 0) && ($COMP(odim)[i] < 0))
$CROAK("Hmmm, both dummy and squished -- this can never happen. I quit.");
/** First handle dummy dims since there's no input from the parent **/
if ($COMP(idim)[i] < 0) {
/* dummy dim - offset or diminc. */
$PDL(CHILD)->dims[ $COMP(odim)[i] ] = $COMP(end)[i] - $COMP(start)[i] + 1;
$PRIV(incs)[ $COMP(odim)[i] ] = 0;
} else {
/** This is not a dummy dim -- deal with a regular slice along it. **/
/** Get parent dim size for this idim, and/or allow permissive slicing **/
PDL_Indx pdsize = $COMP(idim)[i] < $PDL(PARENT)->ndims ?
$PDL(PARENT)->dims[$COMP(idim)[i]] : 1;
PDL_Indx start = $COMP(start)[i];
PDL_Indx end = $COMP(end)[i];
if (
/** Trap special case: full slices of an empty dim are empty **/
(pdsize==0 && start==0 && end==-1 && ($COMP(inc)[i] == 0 || $COMP(inc)[i] == 1))
||
/* the values given when PDL::slice gets empty ndarray for index */
(start==1 && end==0 && $COMP(inc)[i] == 1)
) {
$PDL(CHILD)->dims[$COMP(odim)[i]] = 0;
$PRIV(incs)[$COMP(odim)[i]] = 0;
} else {
/** Regularize and bounds-check the start location **/
if (start < 0)
start += pdsize;
if (start < 0 || start >= pdsize) {
if (i >= $PDL(PARENT)->ndims) {
$CROAK("slice has too many dims (indexes dim %"IND_FLAG"; highest is %"IND_FLAG")",i,$PDL(PARENT)->ndims-1);
} else {
$CROAK("slice starts out of bounds in pos %"IND_FLAG" (start is %"IND_FLAG"; end is %"IND_FLAG"; inc is %"IND_FLAG"; source dim %"IND_FLAG" runs 0 to %"IND_FLAG")",i,start,end,$COMP(inc)[i],$COMP(idim)[i],pdsize-1);
}
}
if ($COMP(odim)[i] < 0) {
/* squished dim - just update the offset. */
/* start is always defined and regularized if we are here here, since */
/* both idim[i] and odim[i] can't be <0 */
$PRIV(offs) += start * $PDL(PARENT)->dimincs[ $COMP(idim)[i] ];
} else /* normal operation */ {
/** Regularize and bounds-check the end location **/
if (end<0) end += pdsize;
if (end < 0 || end >= pdsize)
$CROAK("slice ends out of bounds in pos %"IND_FLAG" (end is %"IND_FLAG"; source dim %"IND_FLAG" runs 0 to %"IND_FLAG")",i,end,$COMP(idim)[i],pdsize-1);
/* regularize inc */
PDL_Indx inc = $COMP(inc)[i];
if (!inc)
inc = (start <= end) ? 1 : -1;
$PDL(CHILD)->dims[ $COMP(odim)[i] ] = PDLMAX(0, (end - start + inc) / inc);
$PRIV(incs)[ $COMP(odim)[i] ] = $PDL(PARENT)->dimincs[ $COMP(idim)[i] ] * inc;
$PRIV(offs) += start * $PDL(PARENT)->dimincs[ $COMP(idim)[i] ];
} /* end of normal slice case */
} /* end of trapped strange slice case */
} /* end of non-dummy slice case */
} /* end of nargs loop */
/* Now fill in broadcast dimensions as needed. idim and odim persist from the parsing loop */
/* up above. */
for(i=0; i<o_ndims_extra; i++) {
$PDL(CHILD)->dims[ $COMP(odim_top) + i ] = $PDL(PARENT)->dims[ $COMP(idim_top) + i ];
$PRIV(incs)[ $COMP(odim_top) + i ] = $PDL(PARENT)->dimincs[ $COMP(idim_top) + i ];
}
$SETDIMS();
EOF
);pp_def('diagonal', P2Child => 1, TwoWay => 1, AffinePriv => 1, OtherPars => 'PDL_Indx whichdims[]', MakeComp => pp_line_numbers(__LINE__-1, ' if ($COMP(whichdims_count) < 1) $CROAK("must have at least 1 dimension"); qsort($COMP(whichdims), $COMP(whichdims_count), sizeof(PDL_Indx), cmp_pdll); '), CHeader => <<'END', static int cmp_pdll(const void *a_,const void *b_) { PDL_Indx *a = (PDL_Indx *)a_; PDL_Indx *b=(PDL_Indx *)b_; if(*a>*b) return 1; else if(*a==*b) return 0; else return -1; } END RedoDims => pp_line_numbers(__LINE__-1, ' PDL_Indx nthp,nthc,nthd, cd = $COMP(whichdims)[0]; $SETNDIMS($PDL(PARENT)->ndims-$COMP(whichdims_count)+1); $DOPRIVALLOC(); $PRIV(offs) = 0; if ($COMP(whichdims)[$COMP(whichdims_count)-1] >= $PDL(PARENT)->ndims || $COMP(whichdims)[0] < 0) $CROAK("dim out of range"); nthd=0; nthc=0; for(nthp=0; nthp<$PDL(PARENT)->ndims; nthp++) if (nthd < $COMP(whichdims_count) && nthp == $COMP(whichdims)[nthd]) { if (!nthd) { $PDL(CHILD)->dims[cd] = $PDL(PARENT)->dims[cd]; nthc++; $PRIV(incs)[cd] = 0; } if (nthd && $COMP(whichdims)[nthd] == $COMP(whichdims)[nthd-1]) $CROAK("dims must be unique"); nthd++; /* advance pointer into whichdims */ if($PDL(CHILD)->dims[cd] != $PDL(PARENT)->dims[nthp]) { $CROAK("Different dims %"IND_FLAG" and %"IND_FLAG"", $PDL(CHILD)->dims[cd], $PDL(PARENT)->dims[nthp]); } $PRIV(incs)[cd] += $PDL(PARENT)->dimincs[nthp]; } else { $PRIV(incs)[nthc] = $PDL(PARENT)->dimincs[nthp]; $PDL(CHILD)->dims[nthc] = $PDL(PARENT)->dims[nthp]; nthc++; } $SETDIMS(); '), PMCode =>pp_line_numbers(__LINE__, <<'EOD'), sub PDL::diagonal :lvalue { shift->_diagonal_int(my $o=PDL->null, \@_); $o } EOD Doc => <<'EOD', =for ref
Returns the multidimensional diagonal over the specified dimensions.
The diagonal is placed at the first (by number) dimension that is diagonalized. The other diagonalized dimensions are removed. So if $x has dimensions (5,3,5,4,6,5) then after
$d = $x->diagonal(dim1, dim2,...)$y = $x->diagonal(0,2,5);the ndarray $y has dimensions (5,3,4,6) and $y->at(2,1,0,1) refers to $x->at(2,1,2,0,1,2).
NOTE: diagonal doesn't handle broadcastids correctly. XXX FIX
pdl> $x = zeroes(3,3,3);
pdl> ($y = $x->diagonal(0,1))++;
pdl> p $x
[
[
[1 0 0]
[0 1 0]
[0 0 1]
]
[
[1 0 0]
[0 1 0]
[0 0 1]
]
[
[1 0 0]
[0 1 0]
[0 0 1]
]
]BUGS
For the moment, you can't slice one of the zero-length dims of an empty ndarray. It is not clear how to implement this in a way that makes sense.
Many types of index errors are reported far from the indexing operation that caused them. This is caused by the underlying architecture: slice() sets up a mapping between variables, but that mapping isn't tested for correctness until it is used (potentially much later).
AUTHOR
Copyright (C) 1997 Tuomas J. Lukka. Contributions by Craig DeForest, deforest@boulder.swri.edu. Documentation contributions by David Mertens. All rights reserved. There is no warranty. You are allowed to redistribute this software / documentation under certain conditions. For details, see the file COPYING in the PDL distribution. If this file is separated from the PDL distribution, the copyright notice should be included in the file.