"Fossies" - the Fresh Open Source Software Archive

Member "PDL-2.080/Basic/Slices/slices.pd" (19 May 2022, 83885 Bytes) of package /linux/misc/PDL-2.080.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. See also the latest Fossies "Diffs" side-by-side code changes report for "slices.pd": 2.079_vs_2.080.

    1 use PDL::Types qw(ppdefs_all);
    2 use strict;
    3 use warnings;
    4 
    5 pp_addpm({At => 'Top'},<< 'EOD');
    6 
    7 =head1 NAME
    8 
    9 PDL::Slices -- Indexing, slicing, and dicing
   10 
   11 =head1 SYNOPSIS
   12 
   13   use PDL;
   14   $x = ones(3,3);
   15   $y = $x->slice('-1:0,(1)');
   16   $c = $x->dummy(2);
   17 
   18 
   19 =head1 DESCRIPTION
   20 
   21 This package provides many of the powerful PerlDL core index
   22 manipulation routines.  These routines mostly allow two-way data flow,
   23 so you can modify your data in the most convenient representation.
   24 For example, you can make a 1000x1000 unit matrix with
   25 
   26  $x = zeroes(1000,1000);
   27  $x->diagonal(0,1) ++;
   28 
   29 which is quite efficient. See L<PDL::Indexing> and L<PDL::Tips> for
   30 more examples.
   31 
   32 Slicing is so central to the PDL language that a special compile-time
   33 syntax has been introduced to handle it compactly; see L<PDL::NiceSlice>
   34 for details.
   35 
   36 PDL indexing and slicing functions usually include two-way data flow,
   37 so that you can separate the actions of reshaping your data structures
   38 and modifying the data themselves.  Two special methods, L</copy> and
   39 L</sever>, help you control the data flow connection between related
   40 variables.
   41 
   42  $y = $x->slice("1:3"); # Slice maintains a link between $x and $y.
   43  $y += 5;               # $x is changed!
   44 
   45 If you want to force a physical copy and no data flow, you can copy or
   46 sever the slice expression:
   47 
   48  $y = $x->slice("1:3")->copy;
   49  $y += 5;               # $x is not changed.
   50 
   51  $y = $x->slice("1:3")->sever;
   52  $y += 5;               # $x is not changed.
   53 
   54 The difference between C<sever> and C<copy> is that sever acts on (and
   55 returns) its argument, while copy produces a disconnected copy.  If you
   56 say
   57 
   58  $y = $x->slice("1:3");
   59  $c = $y->sever;
   60 
   61 then the variables C<$y> and C<$c> point to the same object but with
   62 C<-E<gt>copy> they would not.
   63 
   64 =cut
   65 
   66 use strict;
   67 use warnings;
   68 use PDL::Core ':Internal';
   69 use Scalar::Util 'blessed';
   70 
   71 EOD
   72 
   73 # $::PP_VERBOSE=1;
   74 
   75 pp_addhdr(<<'EOH');
   76 
   77 #ifdef _MSC_VER
   78 #define strtoll _strtoi64
   79 #endif
   80 
   81 EOH
   82 
   83 my $doc = <<'EOD';
   84 =for ref
   85 
   86 C<index>, C<index1d>, and C<index2d> provide rudimentary index indirection.
   87 
   88 =for example
   89 
   90  $c = index($source,$ind);
   91  $c = index1d($source,$ind);
   92  $c = index2d($source2,$ind1,$ind2);
   93 
   94 use the C<$ind> variables as indices to look up values in C<$source>.
   95 The three routines broadcast slightly differently.
   96 
   97 =over 3
   98 
   99 =item *
  100 
  101 C<index> uses direct broadcasting for 1-D indexing across the 0 dim
  102 of C<$source>.  It can broadcast over source broadcast dims or index broadcast
  103 dims, but not (easily) both: If C<$source> has more than 1
  104 dimension and C<$ind> has more than 0 dimensions, they must agree in
  105 a broadcasting sense.
  106 
  107 =item * 
  108 
  109 C<index1d> uses a single active dim in C<$ind> to produce a list of
  110 indexed values in the 0 dim of the output - it is useful for
  111 collapsing C<$source> by indexing with a single row of values along
  112 C<$source>'s 0 dimension.  The output has the same number of dims as
  113 C<$source>.  The 0 dim of the output has size 1 if C<$ind> is a
  114 scalar, and the same size as the 0 dim of C<$ind> if it is not. If
  115 C<$ind> and C<$source> both have more than 1 dim, then all dims higher
  116 than 0 must agree in a broadcasting sense.
  117 
  118 =item * 
  119 
  120 C<index2d> works like C<index> but uses separate ndarrays for X and Y
  121 coordinates.  For more general N-dimensional indexing, see the
  122 L<PDL::NiceSlice> syntax or L<PDL::Slices> (in particular C<slice>,
  123 C<indexND>, and C<range>).
  124 
  125 =back 
  126 
  127 These functions are two-way, i.e. after
  128 
  129  $c = $x->index(pdl[0,5,8]);
  130  $c .= pdl [0,2,4];
  131 
  132 the changes in C<$c> will flow back to C<$x>.
  133 
  134 C<index> provids simple broadcasting:  multiple-dimensioned arrays are treated
  135 as collections of 1-D arrays, so that
  136 
  137  $x = xvals(10,10)+10*yvals(10,10);
  138  $y = $x->index(3);
  139  $c = $x->index(9-xvals(10));
  140 
  141 puts a single column from C<$x> into C<$y>, and puts a single element
  142 from each column of C<$x> into C<$c>.  If you want to extract multiple
  143 columns from an array in one operation, see L</dice> or
  144 L</indexND>.
  145 
  146 =cut
  147 
  148 EOD
  149 
  150 my $index_init =
  151        'register PDL_Indx this_ind = $ind();
  152         if( PDL_IF_BAD($ISBADVAR(this_ind,ind) ||,) this_ind<0 || this_ind>=$SIZE(n) ) {
  153            $CROAK("invalid index %"IND_FLAG" (valid range 0..%"IND_FLAG")",
  154                 this_ind,$SIZE(n)-1);
  155         }';
  156 pp_def(
  157        'index',
  158        GenericTypes => [ppdefs_all],
  159        HandleBad => 1,
  160        DefaultFlow => 1,
  161        TwoWay => 1,
  162        Pars => 'a(n); indx ind(); [oca] c();',
  163        Code =>
  164        $index_init . ' $c() = $a(n => this_ind);',
  165        BackCode =>
  166        $index_init . ' $a(n => this_ind) = $c();',
  167        Doc => $doc,
  168        BadDoc =>
  169        'index barfs if any of the index values are bad.',
  170        );
  171 
  172 pp_def(
  173        'index1d',
  174        GenericTypes => [ppdefs_all],
  175        HandleBad => 1,
  176        DefaultFlow => 1,
  177        TwoWay => 1,
  178        Pars => 'a(n); indx ind(m); [oca] c(m);',
  179        Code => q{
  180          PDL_Indx i;
  181          for(i=0;i<$SIZE(m);i++) {
  182                  PDL_Indx this_ind = $ind(m=>i);
  183                  PDL_IF_BAD(if( $ISBADVAR(this_ind, ind) ) {
  184                      $SETBAD(c(m=>i));
  185                  } else,) {
  186                    if( this_ind<0 || this_ind >= $SIZE(n) ) {
  187                      $CROAK("invalid index %"IND_FLAG" at pos %"IND_FLAG" (valid range 0..%"IND_FLAG")",
  188                       this_ind, i, $SIZE(n)-1);
  189                    }
  190                    $c(m=>i) = $a(n=>this_ind);
  191                  }
  192          }
  193         },
  194        BackCode => q{
  195          PDL_Indx i;
  196          for(i=0;i<$SIZE(m);i++) {
  197                  PDL_Indx this_ind = $ind(m=>i);
  198                  PDL_IF_BAD(if( $ISBADVAR(this_ind, ind) ) {
  199                    /* do nothing */
  200                  } else,) {
  201                    if( this_ind<0 || this_ind >= $SIZE(n) ) {
  202                      $CROAK("invalid index %"IND_FLAG" at pos %"IND_FLAG" (valid range 0..%"IND_FLAG")",
  203                       this_ind, i, $SIZE(n)-1);
  204                    }
  205                    $a(n=>this_ind) = $c(m=>i);
  206                  }
  207          }
  208         },
  209        Doc => $doc,
  210        BadDoc =>
  211        'index1d propagates BAD index elements to the output variable.'
  212        );
  213 
  214 my $index2d_init =
  215        'register PDL_Indx this_ind_a = $inda(),this_ind_b = $indb();
  216         if( PDL_IF_BAD($ISBADVAR(this_ind_a,inda) ||,) this_ind_a<0 || this_ind_a>=$SIZE(na) )
  217            $CROAK("invalid x-index %"IND_FLAG" (valid range 0..%"IND_FLAG")",
  218                 this_ind_a,$SIZE(na)-1);
  219         if( PDL_IF_BAD($ISBADVAR(this_ind_b,indb) ||,) this_ind_b<0 || this_ind_b>=$SIZE(nb) )
  220            $CROAK("invalid y-index %"IND_FLAG" (valid range 0..%"IND_FLAG")",
  221                 this_ind_b,$SIZE(nb)-1);
  222         ';
  223 pp_def(
  224        'index2d',
  225        GenericTypes => [ppdefs_all],
  226        HandleBad => 1,
  227        DefaultFlow => 1,
  228        TwoWay => 1,
  229        Pars => 'a(na,nb); indx inda(); indx indb(); [oca] c();',
  230        Code =>
  231        $index2d_init . ' $c() = $a(na => this_ind_a, nb => this_ind_b);',
  232        BackCode =>
  233        $index2d_init . ' $a(na => this_ind_a, nb => this_ind_b) = $c();',
  234        Doc => $doc,
  235        BadDoc =>
  236        'index2d barfs if either of the index values are bad.',
  237 );
  238 
  239 # indexND: CED 2-Aug-2002
  240 pp_add_exported('','indexND indexNDb');
  241 pp_addpm(<<'EOD-indexND');
  242 =head2 indexNDb
  243 
  244 =for ref
  245 
  246   Backwards-compatibility alias for indexND
  247 
  248 =head2 indexND
  249 
  250 =for ref
  251 
  252   Find selected elements in an N-D ndarray, with optional boundary handling
  253 
  254 =for example
  255 
  256   $out = $source->indexND( $index, [$method] )
  257 
  258   $source = 10*xvals(10,10) + yvals(10,10);
  259   $index  = pdl([[2,3],[4,5]],[[6,7],[8,9]]);
  260   print $source->indexND( $index );
  261 
  262   [
  263    [23 45]
  264    [67 89]
  265   ]
  266 
  267 IndexND collapses C<$index> by lookup into C<$source>.  The
  268 0th dimension of C<$index> is treated as coordinates in C<$source>, and
  269 the return value has the same dimensions as the rest of C<$index>.
  270 The returned elements are looked up from C<$source>.  Dataflow
  271 works -- propagated assignment flows back into C<$source>.
  272 
  273 IndexND and IndexNDb were originally separate routines but they are both
  274 now implemented as a call to L</range>, and have identical syntax to
  275 one another.
  276 
  277 SEE ALSO:
  278 
  279 L<PDL::Primitive/whichND> returns N-D indices into a multidimensional
  280 PDL, suitable for feeding to this.
  281 
  282 =cut
  283 
  284 sub PDL::indexND {
  285         my($source,$index, $boundary) = @_;
  286         return PDL::range($source,$index,undef,$boundary);
  287 }
  288 
  289 *PDL::indexNDb = \&PDL::indexND;
  290 
  291 EOD-indexND
  292 
  293 pp_addpm(<<'EOD-range');
  294 
  295 sub PDL::range {
  296   my($source,$ind,$sz,$bound) = @_;
  297 
  298 # Convert to indx type up front (also handled in rangeb if necessary)
  299   my $index = (ref $ind && UNIVERSAL::isa($ind,'PDL') && $ind->type eq 'indx') ? $ind : indx($ind);
  300   my $size = defined($sz) ? PDL->pdl($sz) : undef;
  301 
  302 
  303   # Handle empty PDL case: return a properly constructed Empty.
  304   if($index->isempty) {
  305       my @sdims= $source->dims;
  306       splice(@sdims, 0, $index->dim(0) + ($index->dim(0)==0)); # added term is to treat Empty[0] like a single empty coordinate
  307       unshift(@sdims, $size->list) if(defined($size));
  308       return PDL->new_from_specification(0 x ($index->ndims-1), @sdims);
  309   }
  310 
  311 
  312   $index = $index->dummy(0,1) unless $index->ndims;
  313 
  314 
  315   # Pack boundary string if necessary
  316   if(defined $bound) {
  317     if(ref $bound eq 'ARRAY') {
  318       my ($s,$el);
  319       foreach $el(@$bound) {
  320         barf "Illegal boundary value '$el' in range"
  321           unless( $el =~ m/^([0123fFtTeEpPmM])/ );
  322         $s .= $1;
  323       }
  324       $bound = $s;
  325     }
  326     elsif($bound !~ m/^[0123ftepx]+$/  && $bound =~ m/^([0123ftepx])/i ) {
  327       $bound = $1;
  328     }
  329   }
  330 
  331   no warnings; # shut up about passing undef into rangeb
  332   $source->rangeb($index,$size,$bound);
  333 }
  334 EOD-range
  335 
  336 pp_def(
  337         'rangeb',
  338         OtherPars => 'pdl *ind_pdl; SV *size; SV *boundary_sv',
  339         Doc => <<'EOD',
  340 =for ref
  341 
  342 Engine for L</range>
  343 
  344 =for example
  345 
  346 Same calling convention as L</range>, but you must supply all
  347 parameters.  C<rangeb> is marginally faster as it makes a direct PP call,
  348 avoiding the perl argument-parsing step.
  349 
  350 =head2 range
  351 
  352 =for ref
  353 
  354 Extract selected chunks from a source ndarray, with boundary conditions
  355 
  356 =for example
  357 
  358         $out = $source->range($index,[$size,[$boundary]])
  359 
  360 Returns elements or rectangular slices of the original ndarray, indexed by
  361 the C<$index> ndarray.  C<$source> is an N-dimensional ndarray, and C<$index> is
  362 an ndarray whose first dimension has size up to N.  Each row of C<$index> is
  363 treated as coordinates of a single value or chunk from C<$source>, specifying
  364 the location(s) to extract.
  365 
  366 If you specify a single index location, then range is essentially an expensive
  367 slice, with controllable boundary conditions.
  368 
  369 B<INPUTS>
  370 
  371 C<$index> and C<$size> can be ndarrays or array refs such as you would
  372 feed to L<zeroes|PDL::Core/zeroes> and its ilk.  If C<$index>'s 0th dimension
  373 has size higher than the number of dimensions in C<$source>, then
  374 C<$source> is treated as though it had trivial dummy dimensions of
  375 size 1, up to the required size to be indexed by C<$index> -- so if
  376 your source array is 1-D and your index array is a list of 3-vectors,
  377 you get two dummy dimensions of size 1 on the end of your source array.
  378 
  379 You can extract single elements or N-D rectangular ranges from C<$source>,
  380 by setting C<$size>.  If C<$size> is undef or zero, then you get a single
  381 sample for each row of C<$index>.  This behavior is similar to
  382 L</indexNDb>, which is in fact implemented as a call to L</range>.
  383 
  384 If C<$size> is positive then you get a range of values from C<$source> at
  385 each location, and the output has extra dimensions allocated for them.
  386 C<$size> can be a scalar, in which case it applies to all dimensions, or an
  387 N-vector, in which case each element is applied independently to the
  388 corresponding dimension in C<$source>.  See below for details.
  389 
  390 C<$boundary> is a number, string, or list ref indicating the type of
  391 boundary conditions to use when ranges reach the edge of C<$source>.  If you
  392 specify no boundary conditions the default is to forbid boundary violations
  393 on all axes.  If you specify exactly one boundary condition, it applies to
  394 all axes.  If you specify more (as elements of a list ref, or as a packed
  395 string, see below), then they apply to dimensions in the order in which they
  396 appear, and the last one applies to all subsequent dimensions.  (This is
  397 less difficult than it sounds; see the examples below).
  398 
  399 =over 3
  400 
  401 =item 0 (synonyms: 'f','forbid') B<(default)>
  402 
  403 Ranges are not allowed to cross the boundary of the original PDL.  Disallowed
  404 ranges throw an error.  The errors are thrown at evaluation time, not
  405 at the time of the range call (this is the same behavior as L</slice>).
  406 
  407 =item 1 (synonyms: 't','truncate')
  408 
  409 Values outside the original ndarray get BAD if you've got bad value
  410 support compiled into your PDL and set the badflag for the source PDL;
  411 or 0 if you haven't (you must set the badflag if you want BADs for out
  412 of bound values, otherwise you get 0).  Reverse dataflow works OK for
  413 the portion of the child that is in-bounds.  The out-of-bounds part of
  414 the child is reset to (BAD|0) during each dataflow operation, but
  415 execution continues.
  416 
  417 =item 2 (synonyms: 'e','x','extend')
  418 
  419 Values that would be outside the original ndarray point instead to the
  420 nearest allowed value within the ndarray.  See the CAVEAT below on
  421 mappings that are not single valued.
  422 
  423 =item 3 (synonyms: 'p','periodic')
  424 
  425 Periodic boundary conditions apply: the numbers in $index are applied,
  426 strict-modulo the corresponding dimensions of $source.  This is equivalent to
  427 duplicating the $source ndarray throughout N-D space.  See the CAVEAT below
  428 about mappings that are not single valued.
  429 
  430 =item 4 (synonyms: 'm','mirror')
  431 
  432 Mirror-reflection periodic boundary conditions apply.  See the CAVEAT
  433 below about mappings that are not single valued.
  434 
  435 =back
  436 
  437 The boundary condition identifiers all begin with unique characters, so
  438 you can feed in multiple boundary conditions as either a list ref or a
  439 packed string.  (The packed string is marginally faster to run).  For
  440 example, the four expressions [0,1], ['forbid','truncate'], ['f','t'],
  441 and 'ft' all specify that violating the boundary in the 0th dimension
  442 throws an error, and all other dimensions get truncated.
  443 
  444 If you feed in a single string, it is interpreted as a packed boundary
  445 array if all of its characters are valid boundary specifiers (e.g. 'pet'),
  446 but as a single word-style specifier if they are not (e.g. 'forbid').
  447 
  448 B<OUTPUT>
  449 
  450 The output broadcasts over both C<$index> and C<$source>.  Because implicit
  451 broadcasting can happen in a couple of ways, a little thought is needed.  The
  452 returned dimension list is stacked up like this:
  453 
  454    (index broadcast dims), (index dims (size)), (source broadcast dims)
  455 
  456 The first few dims of the output correspond to the extra dims of
  457 C<$index> (beyond the 0 dim). They allow you to pick out individual
  458 ranges from a large, broadcasted collection.
  459 
  460 The middle few dims of the output correspond to the size dims
  461 specified in C<$size>, and contain the range of values that is extracted
  462 at each location in C<$source>.  Every nonzero element of C<$size> is copied to
  463 the dimension list here, so that if you feed in (for example) C<$size
  464 = [2,0,1]> you get an index dim list of C<(2,1)>.
  465 
  466 The last few dims of the output correspond to extra dims of C<$source> beyond
  467 the number of dims indexed by C<$index>.  These dims act like ordinary
  468 broadcast dims, because adding more dims to C<$source> just tacks extra dims
  469 on the end of the output.  Each source broadcast dim ranges over the entire
  470 corresponding dim of C<$source>.
  471 
  472 B<Dataflow>: Dataflow is bidirectional.
  473 
  474 B<Examples>:
  475 Here are basic examples of C<range> operation, showing how to get
  476 ranges out of a small matrix.  The first few examples show extraction
  477 and selection of individual chunks.  The last example shows
  478 how to mark loci in the original matrix (using dataflow).
  479 
  480  pdl> $src = 10*xvals(10,5)+yvals(10,5)
  481  pdl> print $src->range([2,3])    # Cut out a single element
  482  23
  483  pdl> print $src->range([2,3],1)  # Cut out a single 1x1 block
  484  [
  485   [23]
  486  ]
  487  pdl> print $src->range([2,3], [2,1]) # Cut a 2x1 chunk
  488  [
  489   [23 33]
  490  ]
  491  pdl> print $src->range([[2,3]],[2,1]) # Trivial list of 1 chunk
  492  [
  493   [
  494    [23]
  495    [33]
  496   ]
  497  ]
  498  pdl> print $src->range([[2,3],[0,1]], [2,1])   # two 2x1 chunks
  499  [
  500   [
  501    [23  1]
  502    [33 11]
  503   ]
  504  ]
  505  pdl> # A 2x2 collection of 2x1 chunks
  506  pdl> print $src->range([[[1,1],[2,2]],[[2,3],[0,1]]],[2,1])
  507  [
  508   [
  509    [
  510     [11 22]
  511     [23  1]
  512    ]
  513    [
  514     [21 32]
  515     [33 11]
  516    ]
  517   ]
  518  ]
  519  pdl> $src = xvals(5,3)*10+yvals(5,3)
  520  pdl> print $src->range(3,1)  # Broadcast over y dimension in $src
  521  [
  522   [30]
  523   [31]
  524   [32]
  525  ]
  526 
  527  pdl> $src = zeroes(5,4);
  528  pdl> $src->range(pdl([2,3],[0,1]),pdl(2,1)) .= xvals(2,2,1) + 1
  529  pdl> print $src
  530  [
  531   [0 0 0 0 0]
  532   [2 2 0 0 0]
  533   [0 0 0 0 0]
  534   [0 0 1 1 0]
  535  ]
  536 
  537 B<CAVEAT>: It's quite possible to select multiple ranges that
  538 intersect.  In that case, modifying the ranges doesn't have a
  539 guaranteed result in the original PDL -- the result is an arbitrary
  540 choice among the valid values.  For some things that's OK; but for
  541 others it's not. In particular, this doesn't work:
  542 
  543     pdl> $photon_list = new PDL::RandVar->sample(500)->reshape(2,250)*10
  544     pdl> $histogram = zeroes(10,10)
  545     pdl> $histogram->range($photon_list,1)++;  #not what you wanted
  546 
  547 The reason is that if two photons land in the same bin, then that bin
  548 doesn't get incremented twice.  (That may get fixed in a later version...)
  549 
  550 B<PERMISSIVE RANGING>: If C<$index> has too many dimensions compared
  551 to C<$source>, then $source is treated as though it had dummy
  552 dimensions of size 1, up to the required number of dimensions.  These
  553 virtual dummy dimensions have the usual boundary conditions applied to
  554 them.
  555 
  556 If the 0 dimension of C<$index> is ludicrously large (if its size is
  557 more than 5 greater than the number of dims in the source PDL) then
  558 range will insist that you specify a size in every dimension, to make
  559 sure that you know what you're doing.  That catches a common error with
  560 range usage: confusing the initial dim (which is usually small) with another
  561 index dim (perhaps of size 1000).
  562 
  563 If the index variable is Empty, then range() always returns the Empty PDL.
  564 If the index variable is not Empty, indexing it always yields a boundary
  565 violation.  All non-barfing conditions are treated as truncation, since
  566 there are no actual data to return.
  567 
  568 B<EFFICIENCY>: Because C<range> isn't an affine transformation (it
  569 involves lookup into a list of N-D indices), it is somewhat
  570 memory-inefficient for long lists of ranges, and keeping dataflow open
  571 is much slower than for affine transformations (which don't have to copy
  572 data around).
  573 
  574 Doing operations on small subfields of a large range is inefficient
  575 because the engine must flow the entire range back into the original
  576 PDL with every atomic perl operation, even if you only touch a single element.
  577 One way to speed up such code is to sever your range, so that PDL
  578 doesn't have to copy the data with each operation, then copy the
  579 elements explicitly at the end of your loop.  Here's an example that
  580 labels each region in a range sequentially, using many small
  581 operations rather than a single xvals assignment:
  582 
  583   ### How to make a collection of small ops run fast with range...
  584   $x =  $data->range($index, $sizes, $bound)->sever;
  585   $aa = $data->range($index, $sizes, $bound);
  586   $x($_ - 1) .= $_ for 1..$x->nelem;    # Lots of little ops
  587   $aa .= $x;
  588 
  589 C<range> is a perl front-end to a PP function, C<rangeb>.  Calling
  590 C<rangeb> is marginally faster but requires that you include all arguments.
  591 
  592 DEVEL NOTES
  593 
  594 * index broadcast dimensions are effectively clumped internally.  This
  595 makes it easier to loop over the index array but a little more brain-bending
  596 to tease out the algorithm.
  597 
  598 =cut
  599 EOD
  600         HandleBad => 1,
  601         TwoWay => 1,
  602         P2Child => 1,
  603 
  604 # rdim: dimensionality of each range (0 dim of index PDL)
  605 #
  606 # ntsize: number of nonzero size dimensions
  607 # sizes:  array of range sizes, indexed (0..rdim-1).  A zero element means
  608 #         that the dimension is omitted from the child dim list.
  609 # corners: parent coordinates of each corner, running fastest over coord index.
  610 #       (indexed 0 .. (nitems-1)*(rdim)+rdim-1)
  611 # nitems: total number of list elements   (product of itdims)
  612 # itdim:  number of index broadcast dimensions
  613 # itdims: Size of each index broadcast dimension, indexed (0..itdim-1)
  614 #
  615 # bsize: Number of independently specified boundary conditions
  616 # nsizes: Number of independently specified range dim sizes
  617 # boundary: Array containing all the boundary condition specs
  618 # indord: Order/size of the indexing dim (0th dim of $index)
  619 
  620         Comp => 'PDL_Indx rdim;
  621                  PDL_Indx nitems;
  622                  PDL_Indx itdim;
  623                  PDL_Indx ntsize;
  624                  PDL_Indx bsize;
  625                  PDL_Indx nsizes;
  626                  PDL_Indx sizes[$COMP(rdim)];
  627                  PDL_Indx itdims[$COMP(itdim)];
  628                  PDL_Indx corners[$COMP(rdim) * $COMP(nitems)];
  629                  char boundary[$COMP(rdim)];
  630                  char size_pdl_destroy;
  631                  char ind_pdl_destroy;
  632                  ',
  633         MakeComp => <<'EOD-MakeComp',
  634 pdl *size_pdl;
  635 PDL_RETERROR(PDL_err, PDL->make_physdims(ind_pdl));
  636 $COMP(size_pdl_destroy) = $COMP(ind_pdl_destroy) = 0;
  637 
  638 /* Generalized empties are ok, but not in the special 0 dim (the index vector) */
  639 if(ind_pdl->dims[0] == 0)
  640     { $CROAK("can't handle Empty indices -- call range instead"); }
  641 
  642 /***
  643  * Ensure that the index is a PDL_Indx.  If there's no loss of information,
  644  * just upgrade it -- otherwise, make a temporary copy.
  645  */
  646 switch(ind_pdl->datatype) {
  647  default:                              /* Most types: */
  648    ind_pdl = PDL->hard_copy(ind_pdl);  /*   copy and fall through */
  649    if (!ind_pdl) $CROAK("Error in hard_copy");
  650    $COMP(ind_pdl_destroy) = 1;
  651  case PDL_SB: case PDL_B: case PDL_S: case PDL_US: case PDL_L: case PDL_UL: case PDL_LL: case PDL_ULL:
  652    PDL_RETERROR(PDL_err, PDL->converttype(ind_pdl,PDL_IND)); /* convert in place. */
  653    break;
  654  case PDL_IND:
  655    PDL_RETERROR(PDL_err, PDL->make_physical(ind_pdl));
  656    break;
  657 }
  658 
  659 /***
  660  * Figure sizes of the COMP arrrays and allocate them.
  661  */
  662 {
  663   PDL_Indx i,nitems;
  664 
  665   $COMP(rdim) = ind_pdl->ndims ? ind_pdl->dims[0] : 1;
  666   for(i=nitems=1; i < ind_pdl->ndims; i++)  /* Accumulate item list size */
  667     nitems *= ind_pdl->dims[i];
  668   $COMP(nitems) = nitems;
  669   $COMP(itdim) = ind_pdl->ndims ? ind_pdl->ndims - 1 : 0;
  670   $DOCOMPALLOC();
  671 }
  672 
  673 /***
  674  * Fill in the boundary condition array
  675  */
  676 {
  677   char *bstr;
  678   STRLEN blen;
  679   bstr = SvPV(boundary_sv,blen);
  680 
  681   if(blen == 0) {
  682     /* If no boundary is specified then every dim gets forbidden */
  683     PDL_Indx i;
  684     for (i=0;i<$COMP(rdim);i++)
  685       $COMP(boundary[i]) = 0;
  686   } else {
  687     PDL_Indx i;
  688     for(i=0;i<$COMP(rdim);i++) {
  689       switch(bstr[i < blen ? i : blen-1 ]) {
  690       case '0': case 'f': case 'F':               /* forbid */
  691         $COMP(boundary[i]) = 0;
  692         break;
  693       case '1': case 't': case 'T':               /* truncate */
  694         $COMP(boundary[i]) = 1;
  695         break;
  696       case '2': case 'e': case 'E':               /* extend */
  697         $COMP(boundary[i]) = 2;
  698         break;
  699       case '3': case 'p': case 'P':               /* periodic */
  700         $COMP(boundary[i]) = 3;
  701         break;
  702       case '4': case 'm': case 'M':               /* mirror */
  703         $COMP(boundary[i]) = 4;
  704         break;
  705       default:
  706         {
  707           /* No need to check if i < blen -- this will barf out the
  708            * first time it gets hit.
  709            */
  710           $CROAK("Unknown boundary condition '%c' in range",bstr[i]);
  711         }
  712         break;
  713       } // end of switch
  714     }
  715   }
  716 }
  717 /***
  718  * Store the sizes of the index-broadcast dims
  719  */
  720 {
  721   PDL_Indx i;
  722   PDL_Indx nd = ind_pdl->ndims - 1;
  723   for(i=0; i < nd ; i++)
  724     $COMP(itdims[i]) = ind_pdl->dims[i+1];
  725 }
  726 
  727 /***
  728  * Check and condition the size ndarray, and store sizes of the ranges
  729  */
  730 {
  731   PDL_Indx i,ntsize;
  732 
  733   if( (size == NULL) || (size == &PL_sv_undef) ) {
  734     // NO size was passed in, not normally executed even if you passed in no size to range(),
  735     // as range() generates a size array...
  736     for(i=0;i<$COMP(rdim);i++)
  737           $COMP(sizes[i]) = 0;
  738 
  739   } else {
  740     /* Normal case with sizes present in a PDL */
  741 
  742     if(!(size_pdl = PDL->SvPDLV(size))) /* assignment */
  743       $CROAK("Unable to convert size to a PDL");
  744 
  745     if(size_pdl->nvals == 0) {
  746       // no values in the size_pdl - Empty or Null.  Just copy 0s to all the range dims
  747       for(i=0;i<$COMP(rdim);i++)
  748         $COMP(sizes[i]) = 0;
  749 
  750     } else {
  751 
  752       // Convert size PDL to PDL_IND to support indices
  753       switch(size_pdl->datatype) {
  754       default:                              /* Most types: */
  755         size_pdl = PDL->hard_copy(size_pdl);  /*   copy and fall through */
  756         if (!size_pdl) $CROAK("Error in hard_copy");
  757         $COMP(size_pdl_destroy) = 1;
  758       case PDL_SB: case PDL_B: case PDL_S: case PDL_US: case PDL_L: case PDL_UL: case PDL_LL: case PDL_ULL:
  759         PDL_RETERROR(PDL_err, PDL->converttype(size_pdl,PDL_IND)); /* convert in place. */
  760         break;
  761       case PDL_IND:
  762         PDL_RETERROR(PDL_err, PDL->make_physical(size_pdl));
  763         break;
  764       }
  765 
  766       $COMP(nsizes) = size_pdl->nvals; /* Store for later permissiveness check */
  767 
  768       /* Copy the sizes, or die if they're the wrong shape */
  769       if(size_pdl->nvals == 1) {
  770         for(i=0;i<$COMP(rdim);i++) {
  771           $COMP(sizes[i]) = *((PDL_Indx *)(size_pdl->data));
  772         }
  773 
  774         /* Check for nonnegativity of sizes.  The rdim>0 mask ensures that */
  775         /* we don't barf on the Empty PDL (as an index). */
  776         if( $COMP(rdim) > 0 && $COMP(sizes[0]) < 0 ) {
  777           $CROAK("Negative range size is not allowed\n");
  778         }
  779       }
  780       else if( size_pdl->nvals <= $COMP(rdim) && size_pdl->ndims == 1) {
  781         for(i=0;i<$COMP(rdim);i++) {
  782           $COMP(sizes[i]) = (   (i < size_pdl->nvals) ?
  783                                 ((PDL_Indx *)(size_pdl->data))[i] :
  784                                 0
  785                             );
  786           if($COMP(sizes[i]) < 0)
  787                 $CROAK("Negative range sizes are not allowed\n");
  788         }
  789       }
  790       else {
  791         $CROAK("Size must match index's 0th dim\n");
  792       }
  793 
  794     } /* end of nonempty size-ndarray code */
  795   } /* end of defined-size-ndarray code */
  796 
  797   /* Insert the number of nontrivial sizes (these get output dimensions) */
  798   for(i=ntsize=0;i<$COMP(rdim);i++)
  799     if($COMP(sizes[i]))
  800       ntsize++;
  801   $COMP(ntsize) = ntsize;
  802 }
  803 
  804 /***
  805  * Stash coordinates of the corners
  806  */
  807 
  808 PDL_Indx j,k,ioff, iter[$COMP(itdim)];
  809 /* initialize iterator to loop over index broadcasts */
  810 PDL_Indx *cptr = iter;
  811 for(k=0;k<$COMP(itdim);k++)
  812   *(cptr++) = 0;
  813 cptr = $COMP(corners);
  814 do {
  815   /* accumulate offset into the index from the iterator */
  816   for(k=ioff=0;k<$COMP(itdim);k++)
  817     ioff += iter[k] * ind_pdl->dimincs[k+1];
  818   /* Loop over the 0th dim of index, copying coords. */
  819   /* This is the natural place to check for permissive ranging; too */
  820   /* bad we don't have access to the parent ndarray here... */
  821   for(j=0;j<$COMP(rdim);j++)
  822       *(cptr++) = ((PDL_Indx *)(ind_pdl->data))[ioff + ind_pdl->dimincs[0] * j];
  823   /* Increment the iterator -- the test increments, the body carries. */
  824   for(k=0; k<$COMP(itdim) && (++(iter[k]))>=($COMP(itdims)[k]) ;k++)
  825     iter[k] = 0;
  826 } while(k<$COMP(itdim));
  827 if ($COMP(ind_pdl_destroy))
  828   PDL->destroy(ind_pdl); /* finished with our copy */
  829 if ($COMP(size_pdl_destroy))
  830   PDL->destroy(size_pdl); /* finished with our copy */
  831 EOD-MakeComp
  832 
  833         RedoDims => <<'EOD-RedoDims' ,
  834 PDL_Indx stdim = $PDL(PARENT)->ndims - $COMP(rdim);
  835 PDL_Indx dim,inc;
  836 PDL_Indx i,rdvalid;
  837 
  838   // Speed bump for ludicrous cases
  839   if( $COMP(rdim) > $PDL(PARENT)->ndims+5 && $COMP(nsizes) != $COMP(rdim)) {
  840     $CROAK(
  841       "Ludicrous number of extra dims in range index; leaving child null.\n"
  842       "  (%"IND_FLAG" implicit dims is > 5; index has %"IND_FLAG" dims; source has %"IND_FLAG" dim%s.)\n"
  843       "  This often means that your index PDL is incorrect.\n"
  844       "  To avoid this message, allocate dummy dims in\n"
  845       "  the source or use %"IND_FLAG" dims in range's size field.\n",
  846       $COMP(rdim)-$PDL(PARENT)->ndims,$COMP(rdim),$PDL(PARENT)->ndims,
  847       $PDL(PARENT)->ndims>1?"s":"",$COMP(rdim)
  848     );
  849   }
  850 
  851   if(stdim < 0)
  852     stdim = 0;
  853 
  854   /* Set dimensionality of child */
  855   $PDL(CHILD)->ndims = $COMP(itdim) + $COMP(ntsize) + stdim;
  856   $SETNDIMS($COMP(itdim)+$COMP(ntsize)+stdim);
  857 
  858   inc = 1;
  859   /* Copy size dimensions to child, crunching as we go. */
  860   dim = $COMP(itdim);
  861   for(i=rdvalid=0;i<$COMP(rdim);i++) {
  862     if($COMP(sizes[i])) {
  863       rdvalid++;
  864       $PDL(CHILD)->dimincs[dim] = inc;
  865       inc *= ($PDL(CHILD)->dims[dim++] = $COMP(sizes[i])); /* assignment */
  866     }
  867   }
  868 
  869   /* Copy index broadcast dimensions to child */
  870   for(dim=0; dim<$COMP(itdim); dim++) {
  871     $PDL(CHILD)->dimincs[dim] = inc;
  872     inc *= ($PDL(CHILD)->dims[dim] = $COMP(itdims[dim])); /* assignment */
  873   }
  874 
  875   /* Copy source broadcast dimensions to child */
  876   dim = $COMP(itdim) + rdvalid;
  877   for(i=0;i<stdim;i++) {
  878     $PDL(CHILD)->dimincs[dim] = inc;
  879     inc *= ($PDL(CHILD)->dims[dim++] = $PDL(PARENT)->dims[i+$COMP(rdim)]); /* assignment */
  880   }
  881 
  882   /* Cover bizarre case where the source PDL is empty - in that case, change */
  883   /* all non-barfing boundary conditions to truncation, since we have no data */
  884   /* to reflect, extend, or mirror. */
  885   if($PDL(PARENT)->dims[0] == 0) {
  886     for(dim=0; dim<$COMP(rdim); dim++) {
  887       if($COMP(boundary[dim]))
  888         $COMP(boundary[dim]) = 1; // force truncation
  889     }
  890   }
  891 
  892 $PDL(CHILD)->datatype = $PDL(PARENT)->datatype;
  893 $SETDIMS();
  894 EOD-RedoDims
  895 
  896         EquivCPOffsCode => <<'EOD-EquivCPOffsCode',
  897 PDL_Indx *ip;  /* vector iterator */
  898 PDL_Indx *sp; /* size vector including stdims */
  899 PDL_Indx *coords;     /* current coordinates */
  900 
  901 PDL_Indx k;           /* index */
  902 PDL_Indx item;        /* index broadcast iterator */
  903 PDL_Indx pdim = $PDL(PARENT)->ndims;
  904 PDL_Indx rdim = $COMP(rdim);
  905 PDL_Indx prdim = PDLMIN(rdim,pdim);
  906 PDL_Indx iter2[pdim * 2 + rdim];
  907 PDL_Indx *sizes  = iter2 + pdim;
  908 coords = sizes + pdim;
  909 
  910 /* Figure out size vector */
  911 for(ip = $COMP(sizes), sp = sizes, k=0; k<rdim; k++)
  912    *(sp++) = *(ip++);
  913 for(; k < pdim; k++)
  914    *(sp++) = $PDL(PARENT)->dims[k];
  915 
  916 
  917 /* Loop over all the ranges in the index list */
  918 for(item=0; item<$COMP(nitems); item++) {
  919 
  920   /* initialize in-range iterator to loop within each range */
  921   for(ip = iter2, k=0; k<pdim; k++)
  922     *(ip++) = 0;
  923 
  924   do {
  925     PDL_Indx poff = 0;
  926     PDL_Indx coff;
  927     PDL_Indx k2;
  928     char trunc = 0;       /* Flag used to skip truncation case */
  929 
  930     /* Collect and boundary-check the current N-D coords */
  931     for(k=0; k < prdim; k++){
  932 
  933       PDL_Indx ck = iter2[k] + $COMP(corners[ item * rdim + k  ]) ;
  934 
  935       /* normal case */
  936         if(ck < 0 || ck >= $PDL(PARENT)->dims[k]) {
  937           switch($COMP(boundary[k])) {
  938           case 0: /* no boundary breakage allowed */
  939             $CROAK("index out-of-bounds in range (index vector #%ld)",item);
  940             break;
  941           case 1: /* truncation */
  942             trunc = 1;
  943             break;
  944           case 2: /* extension -- crop */
  945             ck = (ck >= $PDL(PARENT)->dims[k]) ? $PDL(PARENT)->dims[k]-1 : 0;
  946             break;
  947           case 3: /* periodic -- mod it */
  948             ck %= $PDL(PARENT)->dims[k];
  949             if(ck < 0)   /* Fix mod breakage in C */
  950               ck += $PDL(PARENT)->dims[k];
  951             break;
  952           case 4: /* mirror -- reflect off the edges */
  953             ck += $PDL(PARENT)->dims[k];
  954             ck %= ($PDL(PARENT)->dims[k] * 2);
  955             if(ck < 0) /* Fix mod breakage in C */
  956               ck += $PDL(PARENT)->dims[k]*2;
  957             ck -= $PDL(PARENT)->dims[k];
  958             if(ck < 0) {
  959               ck *= -1;
  960               ck -= 1;
  961             }
  962             break;
  963           default:
  964             $CROAK("Unknown boundary condition in range -- bug alert!");
  965             break;
  966           }
  967         }
  968 
  969       coords[k] = ck;
  970 
  971     }
  972 
  973     /* Check extra dimensions -- pick up where k left off... */
  974     for( ; k < rdim ; k++) {
  975       /* Check for indexing off the end of the dimension list */
  976 
  977       PDL_Indx ck = iter2[k] + $COMP(corners[ item * rdim + k  ]) ;
  978 
  979       switch($COMP(boundary[k])) {
  980           case 0: /* No boundary breakage allowed -- nonzero corners cause barfage */
  981             if(ck != 0)
  982                $CROAK("Too many dims in range index (and you've forbidden boundary violations)");
  983             break;
  984           case 1: /* truncation - just truncate if the corner is nonzero */
  985             trunc |= (ck != 0);
  986             break;
  987           case 2: /* extension -- ignore the corner (same as 3) */
  988           case 3: /* periodic  -- ignore the corner */
  989           case 4: /* mirror -- ignore the corner */
  990             ck = 0;
  991             break;
  992           default:
  993             $CROAK("Unknown boundary condition in range -- bug alert!");
  994             break;
  995         }
  996     }
  997 
  998     /* Find offsets into the child and parent arrays, from the N-D coords */
  999     /* Note we only loop over real source dims (prdim) to accumulate -- */
 1000     /* because the offset is trivial and/or we're truncating for virtual */
 1001     /* dims caused by permissive ranging. */
 1002     coff = $PDL(CHILD)->dimincs[0] * item;
 1003     for(k2 = $COMP(itdim), poff = k = 0;
 1004         k < prdim;
 1005         k++) {
 1006       poff += coords[k]*$PDL(PARENT)->dimincs[k];
 1007       if($COMP(sizes[k]))
 1008         coff += iter2[k] * $PDL(CHILD)->dimincs[k2++];
 1009     }
 1010 
 1011     /* Loop the copy over all the source broadcast dims (above rdim). */
 1012     do {
 1013       PDL_Indx poff1 = poff;
 1014       PDL_Indx coff1 = coff;
 1015 
 1016       /* Accumulate the offset due to source broadcasting */
 1017       for(k2 = $COMP(itdim) + $COMP(ntsize), k = rdim;
 1018           k < pdim;
 1019           k++) {
 1020         poff1 += iter2[k] * $PDL(PARENT)->dimincs[k];
 1021         coff1 += iter2[k] * $PDL(CHILD)->dimincs[k2++];
 1022       }
 1023 
 1024       /* Finally -- make the copy
 1025        * EQUIVCPTRUNC works like EQUIVCPOFFS but with checking for
 1026        * out-of-bounds conditions.
 1027        */
 1028       $EQUIVCPTRUNC(coff1,poff1,trunc);
 1029 
 1030       /* Increment the source broadcast iterator */
 1031       for( k=$COMP(rdim);
 1032            k < pdim && (++(iter2[k]) >= $PDL(PARENT)->dims[k]);
 1033            k++)
 1034         iter2[k] = 0;
 1035     } while(k < pdim); /* end of source-broadcast iteration */
 1036 
 1037     /* Increment the in-range iterator */
 1038     for(k = 0;
 1039         k < $COMP(rdim) && (++(iter2[k]) >= $COMP(sizes[k]));
 1040         k++)
 1041       iter2[k] = 0;
 1042   } while(k < $COMP(rdim)); /* end of main iteration */
 1043 } /* end of item do loop */
 1044 EOD-EquivCPOffsCode
 1045 );
 1046 
 1047 pp_def(
 1048         'rld',
 1049         GenericTypes => [ppdefs_all],
 1050         Pars=>'indx a(n); b(n); [o]c(m);',
 1051         PMCode =><<'EOD',
 1052 sub PDL::rld {
 1053   my ($x,$y) = @_;
 1054   my ($c);
 1055   if ($#_ == 2) {
 1056     $c = $_[2];
 1057   } else {
 1058 # XXX Need to improve emulation of broadcasting in auto-generating c
 1059     my ($size) = $x->sumover->max->sclr;
 1060     my (@dims) = $x->dims;
 1061     shift @dims;
 1062     $c = $y->zeroes($size,@dims);
 1063   }
 1064   &PDL::_rld_int($x,$y,$c);
 1065   $c;
 1066 }
 1067 EOD
 1068         Code=>'
 1069           PDL_Indx i,j=0,an;
 1070           $GENERIC(b) bv;
 1071           loop (n) %{
 1072             an = $a();
 1073             bv = $b();
 1074             for (i=0;i<an;i++) {
 1075               $c(m=>j) = bv;
 1076               j++;
 1077             }
 1078           %}',
 1079         Doc => <<'EOD'
 1080 =for ref
 1081 
 1082 Run-length decode a vector
 1083 
 1084 Given a vector C<$x> of the numbers of instances of values C<$y>, run-length
 1085 decode to C<$c>.
 1086 
 1087 =for example
 1088 
 1089  rld($x,$y,$c=null);
 1090 
 1091 =cut
 1092 
 1093 EOD
 1094 );
 1095 
 1096 pp_def(
 1097         'rle',
 1098         GenericTypes => [ppdefs_all],
 1099         Pars=>'c(n); indx [o]a(m); [o]b(m);',
 1100         RedoDimsCode=>'$SIZE(m)=$SIZE(n);',
 1101         PMCode=><<'EOC',
 1102 sub PDL::rle {
 1103   my $c = shift;
 1104   my ($x,$y) = @_==2 ? @_ : (null,null);
 1105   PDL::_rle_int($c,$x,$y);
 1106   my $max_ind = ($c->ndims<2) ? ($x!=0)->sumover-1 :
 1107                                 ($x!=0)->clump(1..$x->ndims-1)->sumover->max->sclr-1;
 1108   return ($x->slice("0:$max_ind"),$y->slice("0:$max_ind"));
 1109 }
 1110 EOC
 1111         Code=>'
 1112           PDL_Indx j=0;
 1113           $GENERIC(c) cv, clv;
 1114           clv = $c(n=>0);
 1115           $b(m=>0) = clv;
 1116           $a(m=>0) = 0;
 1117           loop (n) %{
 1118             cv = $c();
 1119             if (cv == clv) {
 1120               $a(m=>j)++;
 1121             } else {
 1122               j++;
 1123               $b(m=>j) = clv = cv;
 1124               $a(m=>j) = 1;
 1125             }
 1126           %}
 1127           for (j++;j<$SIZE(m);j++) {
 1128             $a(m=>j) = 0;
 1129             $b(m=>j) = 0;
 1130           }
 1131         ',
 1132         Doc => <<'EOD'
 1133 =for ref
 1134 
 1135 Run-length encode a vector
 1136 
 1137 Given vector C<$c>, generate a vector C<$x> with the number of each
 1138 element, and a vector C<$y> of the unique values.  New in PDL 2.017,
 1139 only the elements up to the first instance of C<0> in C<$x> are
 1140 returned, which makes the common use case of a 1-dimensional C<$c> simpler.
 1141 For broadcast operation, C<$x> and C<$y> will be large enough
 1142 to hold the largest row of C<$y>, and only the elements up to the
 1143 first instance of C<0> in each row of C<$x> should be considered.
 1144 
 1145 =for example
 1146 
 1147  $c = floor(4*random(10));
 1148  rle($c,$x=null,$y=null);
 1149  #or
 1150  ($x,$y) = rle($c);
 1151 
 1152  #for $c of shape [10, 4]:
 1153  $c = floor(4*random(10,4));
 1154  ($x,$y) = rle($c);
 1155 
 1156  #to see the results of each row one at a time:
 1157  foreach (0..$c->dim(1)-1){
 1158   my ($as,$bs) = ($x(:,($_)),$y(:,($_)));
 1159   my ($ta,$tb) = where($as,$bs,$as!=0); #only the non-zero elements of $x
 1160   print $c(:,($_)) . " rle==> " , ($ta,$tb) , "\trld==> " . rld($ta,$tb) . "\n";
 1161  }
 1162 
 1163  # the inverse of (chance of all 6 3d6 rolls being >= each possible sum)
 1164  ($nrolls, $ndice, $dmax) = (6, 3, 6);
 1165  ($x, $x1) = (allaxisvals(($dmax) x $ndice)+1)->sumover->flat->qsort->rle;
 1166  $y = $x->cumusumover;
 1167  $yprob1x = $y->slice('-1:0')->double / $y->slice('(-1)');
 1168  $z = cat($x1, 1 / $yprob1x**$nrolls)->transpose;
 1169 
 1170 =cut
 1171 
 1172 EOD
 1173 );
 1174 
 1175 pp_def('rlevec',
 1176        Pars => "c(M,N); indx [o]a(N); [o]b(M,N)",
 1177        Code =><<'EOC',
 1178   PDL_Indx cn,bn=0, sn=$SIZE(N), matches;
 1179   loop (M) %{ $b(N=>0)=$c(N=>0); %}
 1180   $a(N=>0) = 1;
 1181   for (cn=1; cn<sn; cn++) {
 1182      matches=1;
 1183      loop (M) %{
 1184        if ($c(N=>cn) != $b(N=>bn)) {
 1185          matches=0;
 1186          break;
 1187        }
 1188      %}
 1189      if (matches) {
 1190        $a(N=>bn)++;
 1191      } else {
 1192        bn++;
 1193        loop (M) %{ $b(N=>bn) = $c(N=>cn); %}
 1194        $a(N=>bn) = 1;
 1195      }
 1196    }
 1197    for (bn++; bn<sn; bn++) {
 1198      $a(N=>bn) = 0;
 1199      loop (M) %{ $b(N=>bn) = 0; %}
 1200    }
 1201 EOC
 1202        Doc =><<'EOD',
 1203 =for ref
 1204 
 1205 Run-length encode a set of vectors.
 1206 
 1207 Higher-order rle(), for use with qsortvec().
 1208 
 1209 Given set of vectors $c, generate a vector $a with the number of occurrences of each element
 1210 (where an "element" is a vector of length $M ocurring in $c),
 1211 and a set of vectors $b containing the unique values.
 1212 As for rle(), only the elements up to the first instance of 0 in $a should be considered.
 1213 
 1214 Can be used together with clump() to run-length encode "values" of arbitrary dimensions.
 1215 Can be used together with rotate(), cat(), append(), and qsortvec() to count N-grams
 1216 over a 1d PDL.
 1217 
 1218 See also: L</rle>, L<PDL::Ufunc/qsortvec>, L<PDL::Primitive/uniqvec>
 1219 Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
 1220 
 1221 EOD
 1222 );
 1223 
 1224 pp_def('rldvec',
 1225        Pars => 'indx a(N); b(M,N); [o]c(M,N)',
 1226        PMCode=><<'EOC',
 1227 sub PDL::rldvec {
 1228   my ($a,$b,$c) = @_;
 1229   if (!defined($c)) {
 1230 # XXX Need to improve emulation of threading in auto-generating c
 1231     my ($rowlen) = $b->dim(0);
 1232     my ($size) = $a->sumover->max;
 1233     my (@dims) = $a->dims;
 1234     shift(@dims);
 1235     $c = $b->zeroes($b->type,$rowlen,$size,@dims);
 1236   }
 1237   &PDL::_rldvec_int($a,$b,$c);
 1238   return $c;
 1239 }
 1240 EOC
 1241        Code =><<'EOC',
 1242   PDL_Indx i,nrows,bn,cn=0, sn=$SIZE(N);
 1243   for (bn=0; bn<sn; bn++) {
 1244     nrows = $a(N=>bn);
 1245     for (i=0; i<nrows; i++) {
 1246       loop (M) %{ $c(N=>cn) = $b(N=>bn); %}
 1247       cn++;
 1248     }
 1249    }
 1250 EOC
 1251        Doc =><<'EOD'
 1252 =for ref
 1253 
 1254 Run-length decode a set of vectors, akin to a higher-order rld().
 1255 
 1256 Given a vector $a() of the number of occurrences of each row, and a set $c()
 1257 of row-vectors each of length $M, run-length decode to $c().
 1258 
 1259 Can be used together with clump() to run-length decode "values" of arbitrary dimensions.
 1260 
 1261 See also: L</rld>.
 1262 Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
 1263 EOD
 1264 );
 1265 
 1266 pp_def('rleseq',
 1267        Pars => "c(N); indx [o]a(N); [o]b(N)",
 1268        Code=><<'EOC',
 1269   PDL_Indx j=0, sizeN=$SIZE(N);
 1270   $GENERIC(c) coff;
 1271   coff     = $c(N=>0);
 1272   $b(N=>0) = coff;
 1273   $a(N=>0) = 0;
 1274   loop (N) %{
 1275     if ($c() == coff+$a(N=>j)) {
 1276       $a(N=>j)++;
 1277     } else {
 1278       j++;
 1279       $b(N=>j) = coff = $c();
 1280       $a(N=>j) = 1;
 1281     }
 1282   %}
 1283   for (j++; j<sizeN; j++) {
 1284     $a(N=>j) = 0;
 1285     $b(N=>j) = 0;
 1286   }
 1287 EOC
 1288        Doc =><<'EOD',
 1289 =for ref
 1290 
 1291 Run-length encode a vector of subsequences.
 1292 
 1293 Given a vector of $c() of concatenated variable-length, variable-offset subsequences,
 1294 generate a vector $a containing the length of each subsequence
 1295 and a vector $b containing the subsequence offsets.
 1296 As for rle(), only the elements up to the first instance of 0 in $a should be considered.
 1297 
 1298 See also L</rle>.
 1299 Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
 1300 EOD
 1301 );
 1302 
 1303 pp_def('rldseq',
 1304        Pars => 'indx a(N); b(N); [o]c(M)',
 1305        PMCode=><<'EOC',
 1306 sub PDL::rldseq {
 1307   my ($a,$b,$c) = @_;
 1308   if (!defined($c)) {
 1309     my $size   = $a->sumover->max;
 1310     my (@dims) = $a->dims;
 1311     shift(@dims);
 1312     $c = $b->zeroes($b->type,$size,@dims);
 1313   }
 1314   &PDL::_rldseq_int($a,$b,$c);
 1315   return $c;
 1316 }
 1317 EOC
 1318        Code =><<'EOC',
 1319   size_t mi=0;
 1320   loop (N) %{
 1321     size_t     len = $a(), li;
 1322     for (li=0; li < len; ++li, ++mi) {
 1323       $c(M=>mi) = $b() + li;
 1324     }
 1325   %}
 1326 EOC
 1327        Doc =><<'EOD'
 1328 =for ref
 1329 
 1330 Run-length decode a subsequence vector.
 1331 
 1332 Given a vector $a() of sequence lengths
 1333 and a vector $b() of corresponding offsets,
 1334 decode concatenation of subsequences to $c(),
 1335 as for:
 1336 
 1337  $c = null;
 1338  $c = $c->append($b($_)+sequence($a->type,$a($_))) foreach (0..($N-1));
 1339 
 1340 See also: L</rld>.
 1341 Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
 1342 EOD
 1343 );
 1344 
 1345 pp_add_exported('','rleND rldND');
 1346 pp_addpm(<<'EOF');
 1347 =head2 rleND
 1348 
 1349 =for sig
 1350 
 1351   Signature: (data(@vdims,N); int [o]counts(N); [o]elts(@vdims,N))
 1352 
 1353 =for ref
 1354 
 1355 Run-length encode a set of (sorted) n-dimensional values.
 1356 
 1357 Generalization of rle() and vv_rlevec():
 1358 given set of values $data, generate a vector $counts with the number of occurrences of each element
 1359 (where an "element" is a matrix of dimensions @vdims ocurring as a sequential run over the
 1360 final dimension in $data), and a set of vectors $elts containing the elements which begin a run.
 1361 Really just a wrapper for clump() and rlevec().
 1362 
 1363 See also: L</rle>, L</rlevec>.
 1364 Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
 1365 
 1366 =cut
 1367 
 1368 *PDL::rleND = \&rleND;
 1369 sub rleND {
 1370   my $data   = shift;
 1371   my @vdimsN = $data->dims;
 1372 
 1373   ##-- construct output pdls
 1374   my $counts = $#_ >= 0 ? $_[0] : zeroes(long, $vdimsN[$#vdimsN]);
 1375   my $elts   = $#_ >= 1 ? $_[1] : zeroes($data->type, @vdimsN);
 1376 
 1377   ##-- guts: call rlevec()
 1378   rlevec($data->clump($#vdimsN), $counts, $elts->clump($#vdimsN));
 1379 
 1380   return ($counts,$elts);
 1381 }
 1382 
 1383 =head2 rldND
 1384 
 1385 =for sig
 1386 
 1387   Signature: (int counts(N); elts(@vdims,N); [o]data(@vdims,N);)
 1388 
 1389 =for ref
 1390 
 1391 Run-length decode a set of (sorted) n-dimensional values.
 1392 
 1393 Generalization of rld() and rldvec():
 1394 given a vector $counts() of the number of occurrences of each @vdims-dimensioned element,
 1395 and a set $elts() of @vdims-dimensioned elements, run-length decode to $data().
 1396 
 1397 Really just a wrapper for clump() and rldvec().
 1398 
 1399 See also: L</rld>, L</rldvec>.
 1400 Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
 1401 
 1402 =cut
 1403 
 1404 *PDL::rldND = \&rldND;
 1405 sub rldND {
 1406   my ($counts,$elts) = (shift,shift);
 1407   my @vdimsN        = $elts->dims;
 1408 
 1409   ##-- construct output pdl
 1410   my ($data);
 1411   if ($#_ >= 0) { $data = $_[0]; }
 1412   else {
 1413     my $size      = $counts->sumover->max; ##-- get maximum size for Nth-dimension for small encodings
 1414     my @countdims = $counts->dims;
 1415     shift(@countdims);
 1416     $data         = zeroes($elts->type, @vdimsN, @countdims);
 1417   }
 1418 
 1419   ##-- guts: call rldvec()
 1420   rldvec($counts, $elts->clump($#vdimsN), $data->clump($#vdimsN));
 1421 
 1422   return $data;
 1423 }
 1424 EOF
 1425 
 1426 # the perl wrapper clump is now defined in Core.pm
 1427 # this is just the low level interface
 1428 pp_def(
 1429         '_clump_int',
 1430         OtherPars => 'PDL_Indx n',
 1431         P2Child => 1,
 1432         RedoDims => 'PDL_Indx i, d1;
 1433 
 1434              /* truncate overly long clumps to just clump existing dimensions */
 1435          if($COMP(n) > $PDL(PARENT)->ndims)
 1436                         $COMP(n) = $PDL(PARENT)->ndims;
 1437 
 1438          if($COMP(n) < -1)
 1439                         $COMP(n) = $PDL(PARENT)->ndims + $COMP(n) + 1;
 1440 
 1441                  PDL_Indx nrem = ($COMP(n) == -1 ? $PDL(PARENT)->broadcastids[0] : $COMP(n));
 1442                  $SETNDIMS($PDL(PARENT)->ndims - nrem + 1);
 1443                  d1=1;
 1444                  for(i=0; i<nrem; i++) {
 1445                         d1 *= $PDL(PARENT)->dims[i];
 1446                  }
 1447                  $PDL(CHILD)->dims[0] = d1;
 1448                  for(; i<$PDL(PARENT)->ndims; i++) {
 1449                         $PDL(CHILD)->dims[i-nrem+1] = $PDL(PARENT)->dims[i];
 1450                  }
 1451                  $SETDIMS();
 1452                  $SETDELTABROADCASTIDS(1-nrem);
 1453                  ',
 1454         EquivCPOffsCode => '
 1455                 PDL_Indx i;
 1456                 for(i=0; i<$PDL(CHILD)->nvals; i++) {
 1457                         $EQUIVCPOFFS(i,i);
 1458                 }
 1459                 ',
 1460         TwoWay => 1,
 1461         Doc => 'internal',
 1462 );
 1463 
 1464 pp_def(
 1465         'xchg',
 1466         OtherPars => 'PDL_Indx n1; PDL_Indx n2;',
 1467         TwoWay => 1,
 1468         P2Child => 1,
 1469         AffinePriv => 1,
 1470         EquivDimCheck => 'if ($COMP(n1) <0)
 1471                                 $COMP(n1) += $PDL(PARENT)->broadcastids[0];
 1472                           if ($COMP(n2) <0)
 1473                                 $COMP(n2) += $PDL(PARENT)->broadcastids[0];
 1474           if (PDLMIN($COMP(n1),$COMP(n2)) <0 ||
 1475               PDLMAX($COMP(n1),$COMP(n2)) >= $PDL(PARENT)->broadcastids[0])
 1476                 $CROAK("One of dims %"IND_FLAG", %"IND_FLAG" out of range: should be 0<=dim<%"IND_FLAG"",
 1477                     $COMP(n1),$COMP(n2),$PDL(PARENT)->broadcastids[0]);',
 1478         EquivPDimExpr => '(($CDIM == $COMP(n1)) ? $COMP(n2) : ($CDIM == $COMP(n2)) ? $COMP(n1) : $CDIM)',
 1479         Doc => <<'EOD',
 1480 =for ref
 1481 
 1482 exchange two dimensions
 1483 
 1484 Negative dimension indices count from the end.
 1485 
 1486 The command
 1487 
 1488 =for example
 1489 
 1490  $y = $x->xchg(2,3);
 1491 
 1492 creates C<$y> to be like C<$x> except that the dimensions 2 and 3
 1493 are exchanged with each other i.e.
 1494 
 1495  $y->at(5,3,2,8) == $x->at(5,3,8,2)
 1496 
 1497 =cut
 1498 
 1499 EOD
 1500 );
 1501 
 1502 pp_addpm(<< 'EOD');
 1503 
 1504 =head2 reorder
 1505 
 1506 =for ref
 1507 
 1508 Re-orders the dimensions of a PDL based on the supplied list.
 1509 
 1510 Similar to the L</xchg> method, this method re-orders the dimensions
 1511 of a PDL. While the L</xchg> method swaps the position of two dimensions,
 1512 the reorder method can change the positions of many dimensions at
 1513 once.
 1514 
 1515 =for usage
 1516 
 1517  # Completely reverse the dimension order of a 6-Dim array.
 1518  $reOrderedPDL = $pdl->reorder(5,4,3,2,1,0);
 1519 
 1520 The argument to reorder is an array representing where the current dimensions
 1521 should go in the new array. In the above usage, the argument to reorder
 1522 C<(5,4,3,2,1,0)>
 1523 indicates that the old dimensions (C<$pdl>'s dims) should be re-arranged to make the
 1524 new pdl (C<$reOrderPDL>) according to the following:
 1525 
 1526    Old Position   New Position
 1527    ------------   ------------
 1528    5              0
 1529    4              1
 1530    3              2
 1531    2              3
 1532    1              4
 1533    0              5
 1534 
 1535 You do not need to specify all dimensions, only a complete set
 1536 starting at position 0.  (Extra dimensions are left where they are).
 1537 This means, for example, that you can reorder() the X and Y dimensions of
 1538 an image, and not care whether it is an RGB image with a third dimension running
 1539 across color plane.
 1540 
 1541 =for example
 1542 
 1543 Example:
 1544 
 1545  pdl> $x = sequence(5,3,2);       # Create a 3-d Array
 1546  pdl> p $x
 1547  [
 1548   [
 1549    [ 0  1  2  3  4]
 1550    [ 5  6  7  8  9]
 1551    [10 11 12 13 14]
 1552   ]
 1553   [
 1554    [15 16 17 18 19]
 1555    [20 21 22 23 24]
 1556    [25 26 27 28 29]
 1557   ]
 1558  ]
 1559  pdl> p $x->reorder(2,1,0); # Reverse the order of the 3-D PDL
 1560  [
 1561   [
 1562    [ 0 15]
 1563    [ 5 20]
 1564    [10 25]
 1565   ]
 1566   [
 1567    [ 1 16]
 1568    [ 6 21]
 1569    [11 26]
 1570   ]
 1571   [
 1572    [ 2 17]
 1573    [ 7 22]
 1574    [12 27]
 1575   ]
 1576   [
 1577    [ 3 18]
 1578    [ 8 23]
 1579    [13 28]
 1580   ]
 1581   [
 1582    [ 4 19]
 1583    [ 9 24]
 1584    [14 29]
 1585   ]
 1586  ]
 1587 
 1588 The above is a simple example that could be duplicated by calling
 1589 C<$x-E<gt>xchg(0,2)>, but it demonstrates the basic functionality of reorder.
 1590 
 1591 As this is an index function, any modifications to the
 1592 result PDL will change the parent.
 1593 
 1594 =cut
 1595 
 1596 sub PDL::reorder {
 1597         my ($pdl,@newDimOrder) = @_;
 1598 
 1599         my $arrayMax = $#newDimOrder;
 1600 
 1601         #Error Checking:
 1602         if( $pdl->getndims < scalar(@newDimOrder) ){
 1603                 my $errString = "PDL::reorder: Number of elements (".scalar(@newDimOrder).") in newDimOrder array exceeds\n";
 1604                 $errString .= "the number of dims in the supplied PDL (".$pdl->getndims.")";
 1605                 barf($errString);
 1606         }
 1607 
 1608         # Check to make sure all the dims are within bounds
 1609         for my $i(0..$#newDimOrder) {
 1610           my $dim = $newDimOrder[$i];
 1611           if($dim < 0 || $dim > $#newDimOrder) {
 1612               my $errString = "PDL::reorder: Dim index $newDimOrder[$i] out of range in position $i\n(range is 0-$#newDimOrder)";
 1613               barf($errString);
 1614           }
 1615         }
 1616 
 1617         # Checking that they are all present and also not duplicated is done by broadcast() [I think]
 1618 
 1619         # a quicker way to do the reorder
 1620         return $pdl->broadcast(@newDimOrder)->unbroadcast(0);
 1621 }
 1622 
 1623 EOD
 1624 
 1625 pp_addhdr(<<'EOF');
 1626 #define EQUIVDIM(dima,dimb,cdim,inc) \
 1627   ((cdim < PDLMIN(dima,dimb) || cdim > PDLMAX(dima,dimb)) ? \
 1628     cdim : ((cdim == dimb) ? dima : cdim + inc))
 1629 EOF
 1630 pp_def(
 1631         'mv',
 1632         OtherPars => 'PDL_Indx n1; PDL_Indx n2;',
 1633         TwoWay => 1,
 1634         P2Child => 1,
 1635         AffinePriv => 1,
 1636         EquivDimCheck => 'if ($COMP(n1) <0)
 1637                                 $COMP(n1) += $PDL(PARENT)->broadcastids[0];
 1638                           if ($COMP(n2) <0)
 1639                                 $COMP(n2) += $PDL(PARENT)->broadcastids[0];
 1640           if (PDLMIN($COMP(n1),$COMP(n2)) <0 ||
 1641               PDLMAX($COMP(n1),$COMP(n2)) >= $PDL(PARENT)->broadcastids[0])
 1642                 $CROAK("One of dims %"IND_FLAG", %"IND_FLAG" out of range: should be 0<=dim<%"IND_FLAG"",
 1643                     $COMP(n1),$COMP(n2),$PDL(PARENT)->broadcastids[0]);',
 1644         EquivPDimExpr => '(
 1645           $COMP(n1) == $COMP(n2) ? $CDIM :
 1646           $COMP(n1) < $COMP(n2) ? EQUIVDIM($COMP(n1),$COMP(n2),$CDIM,1) :
 1647           EQUIVDIM($COMP(n1),$COMP(n2),$CDIM,-1)
 1648         )',
 1649         Doc => << 'EOD',
 1650 =for ref
 1651 
 1652 move a dimension to another position
 1653 
 1654 The command
 1655 
 1656 =for example
 1657 
 1658  $y = $x->mv(4,1);
 1659 
 1660 creates C<$y> to be like C<$x> except that the dimension 4 is moved to the
 1661 place 1, so:
 1662 
 1663  $y->at(1,2,3,4,5,6) == $x->at(1,5,2,3,4,6);
 1664 
 1665 The other dimensions are moved accordingly.
 1666 Negative dimension indices count from the end.
 1667 =cut
 1668 EOD
 1669 );
 1670 
 1671 pp_addhdr << 'EOH';
 1672 #define sign(x) ( (x) < 0 ? -1 : 1)
 1673 EOH
 1674 
 1675 pp_addpm(<<'EOD'
 1676 
 1677 =head2 using
 1678 
 1679 =for ref
 1680 
 1681 Returns array of column numbers requested
 1682 
 1683 =for usage
 1684 
 1685  line $pdl->using(1,2);
 1686 
 1687 Plot, as a line, column 1 of C<$pdl> vs. column 2
 1688 
 1689 =for example
 1690 
 1691  pdl> $pdl = rcols("file");
 1692  pdl> line $pdl->using(1,2);
 1693 
 1694 =cut
 1695 
 1696 *using = \&PDL::using;
 1697 sub PDL::using {
 1698   my ($x,@ind)=@_;
 1699   @ind = list $ind[0] if (blessed($ind[0]) && $ind[0]->isa('PDL'));
 1700   foreach (@ind) {
 1701     $_ = $x->slice("($_)");
 1702   }
 1703   @ind;
 1704 }
 1705 
 1706 EOD
 1707 );
 1708 
 1709 pp_add_exported('', 'using');
 1710 
 1711 pp_addhdr(<<END
 1712 static int cmp_pdll(const void *a_,const void *b_) {
 1713         PDL_Indx *a = (PDL_Indx *)a_; PDL_Indx *b=(PDL_Indx *)b_;
 1714         if(*a>*b) return 1;
 1715         else if(*a==*b) return 0;
 1716         else return -1;
 1717 }
 1718 END
 1719 );
 1720 
 1721 pp_def(
 1722         'diagonal',
 1723         P2Child => 1,
 1724         TwoWay => 1,
 1725         AffinePriv => 1,
 1726         OtherPars => 'PDL_Indx whichdims[]',
 1727         MakeComp => pp_line_numbers(__LINE__-1, '
 1728                 if ($COMP(whichdims_count) < 1)
 1729                   $CROAK("must have at least 1 dimension");
 1730                 qsort($COMP(whichdims), $COMP(whichdims_count), sizeof(PDL_Indx),
 1731                         cmp_pdll);
 1732         '),
 1733         RedoDims => pp_line_numbers(__LINE__-1, '
 1734                 PDL_Indx nthp,nthc,nthd, cd = $COMP(whichdims[0]);
 1735                 $SETNDIMS($PDL(PARENT)->ndims-$COMP(whichdims_count)+1);
 1736                 $DOPRIVALLOC();
 1737                 $PRIV(offs) = 0;
 1738                 if ($COMP(whichdims)[$COMP(whichdims_count)-1] >= $PDL(PARENT)->ndims ||
 1739                         $COMP(whichdims)[0] < 0)
 1740                         $CROAK("dim out of range");
 1741                 nthd=0; nthc=0;
 1742                 for(nthp=0; nthp<$PDL(PARENT)->ndims; nthp++)
 1743                         if (nthd < $COMP(whichdims_count) &&
 1744                             nthp == $COMP(whichdims)[nthd]) {
 1745                                 if (!nthd) {
 1746                                         $PDL(CHILD)->dims[cd] = $PDL(PARENT)->dims[cd];
 1747                                         nthc++;
 1748                                         $PRIV(incs)[cd] = 0;
 1749                                 }
 1750                                 if (nthd && $COMP(whichdims)[nthd] ==
 1751                                     $COMP(whichdims)[nthd-1])
 1752                                        $CROAK("dims must be unique");
 1753                                 nthd++; /* advance pointer into whichdims */
 1754                                 if($PDL(CHILD)->dims[cd] !=
 1755                                     $PDL(PARENT)->dims[nthp]) {
 1756                                         $CROAK("Different dims %"IND_FLAG" and %"IND_FLAG"",
 1757                                                 $PDL(CHILD)->dims[cd],
 1758                                                 $PDL(PARENT)->dims[nthp]);
 1759                                 }
 1760                                 $PRIV(incs)[cd] += $PDL(PARENT)->dimincs[nthp];
 1761                         } else {
 1762                                 $PRIV(incs)[nthc] = $PDL(PARENT)->dimincs[nthp];
 1763                                 $PDL(CHILD)->dims[nthc] = $PDL(PARENT)->dims[nthp];
 1764                                 nthc++;
 1765                         }
 1766                 $SETDIMS();
 1767         '),
 1768         PMCode => << 'EOD',
 1769 sub PDL::diagonal { shift->_diagonal_int(my $o=PDL->null, \@_); $o }
 1770 EOD
 1771         Doc => << 'EOD',
 1772 =for ref
 1773 
 1774 Returns the multidimensional diagonal over the specified dimensions.
 1775 
 1776 The diagonal is placed at the first (by number) dimension that is
 1777 diagonalized.
 1778 The other diagonalized dimensions are removed. So if C<$x> has dimensions
 1779 C<(5,3,5,4,6,5)> then after
 1780 
 1781 =for usage
 1782 
 1783  $d = $x->diagonal(dim1, dim2,...)
 1784 
 1785 =for example
 1786 
 1787  $y = $x->diagonal(0,2,5);
 1788 
 1789 the ndarray C<$y> has dimensions C<(5,3,4,6)> and
 1790 C<$y-E<gt>at(2,1,0,1)> refers
 1791 to C<$x-E<gt>at(2,1,2,0,1,2)>.
 1792 
 1793 NOTE: diagonal doesn't handle broadcastids correctly. XXX FIX
 1794 
 1795  pdl> $x = zeroes(3,3,3);
 1796  pdl> ($y = $x->diagonal(0,1))++;
 1797  pdl> p $x
 1798  [
 1799   [
 1800    [1 0 0]
 1801    [0 1 0]
 1802    [0 0 1]
 1803   ]
 1804   [
 1805    [1 0 0]
 1806    [0 1 0]
 1807    [0 0 1]
 1808   ]
 1809   [
 1810    [1 0 0]
 1811    [0 1 0]
 1812    [0 0 1]
 1813   ]
 1814  ]
 1815 =cut
 1816 EOD
 1817 );
 1818 
 1819 pp_def(
 1820         'lags',
 1821         Doc => <<'EOD',
 1822 =for ref
 1823 
 1824 Returns an ndarray of lags to parent.
 1825 
 1826 Usage:
 1827 
 1828 =for usage
 1829 
 1830   $lags = $x->lags($nthdim,$step,$nlags);
 1831 
 1832 I.e. if C<$x> contains
 1833 
 1834  [0,1,2,3,4,5,6,7]
 1835 
 1836 then
 1837 
 1838 =for example
 1839 
 1840  $y = $x->lags(0,2,2);
 1841 
 1842 is a (5,2) matrix
 1843 
 1844  [2,3,4,5,6,7]
 1845  [0,1,2,3,4,5]
 1846 
 1847 This order of returned indices is kept because the function is
 1848 called "lags" i.e. the nth lag is n steps behind the original.
 1849 
 1850 C<$step> and C<$nlags> must be positive. C<$nthdim> can be
 1851 negative and will then be counted from the last dim backwards
 1852 in the usual way (-1 = last dim).
 1853 =cut
 1854 EOD
 1855         P2Child => 1,
 1856         TwoWay => 1,
 1857         AffinePriv => 1,
 1858         OtherPars => join('', map "PDL_Indx $_;", qw(nthdim step n)),
 1859         RedoDims => '
 1860                 PDL_Indx i;
 1861                 if ($COMP(nthdim) < 0)  /* the usual conventions */
 1862                    $COMP(nthdim) += $PDL(PARENT)->ndims;
 1863                 if ($COMP(nthdim) < 0 || $COMP(nthdim) >= $PDL(PARENT)->ndims)
 1864                    $CROAK("dim out of range");
 1865                 if ($COMP(n) < 1)
 1866                    $CROAK("number of lags must be positive");
 1867                 if ($COMP(step) < 1)
 1868                    $CROAK("step must be positive");
 1869                 $PRIV(offs) = 0;
 1870                 $SETNDIMS($PDL(PARENT)->ndims+1);
 1871                 $DOPRIVALLOC();
 1872                 for(i=0; i<$COMP(nthdim); i++) {
 1873                         $PDL(CHILD)->dims[i] = $PDL(PARENT)->dims[i];
 1874                         $PRIV(incs)[i] = $PDL(PARENT)->dimincs[i];
 1875                 }
 1876                 $PDL(CHILD)->dims[i] = $PDL(PARENT)->dims[i] - $COMP(step) * ($COMP(n)-1);
 1877                 if ($PDL(CHILD)->dims[i] < 1)
 1878                   $CROAK("product of step size and number of lags too large");
 1879                 $PDL(CHILD)->dims[i+1] = $COMP(n);
 1880                 $PRIV(incs)[i] = ($PDL(PARENT)->dimincs[i]);
 1881                 $PRIV(incs)[i+1] = - $PDL(PARENT)->dimincs[i] * $COMP(step);
 1882                 $PRIV(offs) += ($PDL(CHILD)->dims[i+1] - 1) * (-$PRIV(incs)[i+1]);
 1883                 i++;
 1884                 for(; i<$PDL(PARENT)->ndims; i++) {
 1885                         $PDL(CHILD)->dims[i+1] = $PDL(PARENT)->dims[i];
 1886                         $PRIV(incs)[i+1] = $PDL(PARENT)->dimincs[i];
 1887                 }
 1888                 $SETDIMS();
 1889         '
 1890 );
 1891 
 1892 pp_def(
 1893         'splitdim',
 1894         Doc => <<'EOD',
 1895 =for ref
 1896 
 1897 Splits a dimension in the parent ndarray (opposite of L<clump|PDL::Core/clump>).
 1898 As of 2.076, throws exception if non-divisible C<nsp> given, and can
 1899 give negative C<nthdim> which then counts backwards.
 1900 
 1901 =for example
 1902 
 1903 After
 1904 
 1905  $y = $x->splitdim(2,3);
 1906 
 1907 the expression
 1908 
 1909  $y->at(6,4,m,n,3,6) == $x->at(6,4,m+3*n)
 1910 
 1911 is always true (C<m> has to be less than 3).
 1912 =cut
 1913 EOD
 1914         P2Child => 1,
 1915         TwoWay => 1,
 1916         OtherPars => join('', map "PDL_Indx $_;", qw(nthdim nsp)),
 1917         AffinePriv => 1,
 1918         RedoDims => '
 1919                 PDL_Indx i = $COMP(nthdim);
 1920                 PDL_Indx nsp = $COMP(nsp);
 1921                 if(nsp == 0) {$CROAK("Cannot split to 0\n");}
 1922                 if (i < 0)
 1923                         i = $COMP(nthdim) += $PDL(PARENT)->ndims;
 1924                 if (i < 0 || i >= $PDL(PARENT)->ndims)
 1925                         $CROAK("nthdim %"IND_FLAG" after adjusting for negative must not be negative or greater or equal to number of dims %"IND_FLAG"\n",
 1926                                 i, $PDL(PARENT)->ndims);
 1927                 if (nsp > $PDL(PARENT)->dims[i])
 1928                         $CROAK("nsp %"IND_FLAG" cannot be greater than dim %"IND_FLAG"\n",
 1929                                 nsp, $PDL(PARENT)->dims[i]);
 1930                 if (($PDL(PARENT)->dims[i] % nsp) != 0)
 1931                         $CROAK("nsp %"IND_FLAG" non-divisible into dim %"IND_FLAG"\n",
 1932                                 nsp, $PDL(PARENT)->dims[i]);
 1933                 $PRIV(offs) = 0;
 1934                 $SETNDIMS($PDL(PARENT)->ndims+1);
 1935                 $DOPRIVALLOC();
 1936                 for(i=0; i<$COMP(nthdim); i++) {
 1937                         $PDL(CHILD)->dims[i] = $PDL(PARENT)->dims[i];
 1938                         $PRIV(incs)[i] = $PDL(PARENT)->dimincs[i];
 1939                 }
 1940                 $PDL(CHILD)->dims[i] = $COMP(nsp);
 1941                 $PDL(CHILD)->dims[i+1] = $PDL(PARENT)->dims[i] / $COMP(nsp);
 1942                 $PRIV(incs)[i] = $PDL(PARENT)->dimincs[i];
 1943                 $PRIV(incs)[i+1] = $PDL(PARENT)->dimincs[i] * $COMP(nsp);
 1944                 i++;
 1945                 for(; i<$PDL(PARENT)->ndims; i++) {
 1946                         $PDL(CHILD)->dims[i+1] = $PDL(PARENT)->dims[i];
 1947                         $PRIV(incs)[i+1] = $PDL(PARENT)->dimincs[i];
 1948                 }
 1949                 $SETDIMS();
 1950         ',
 1951 );
 1952 
 1953 my $rotate_code = '
 1954         PDL_Indx i,j;
 1955         PDL_Indx n_size = $SIZE(n);
 1956         if (n_size == 0)
 1957           $CROAK("can not shift zero size ndarray (n_size is zero)");
 1958         j = ($shift()) % n_size;
 1959         if (j < 0)
 1960                 j += n_size;
 1961         for(i=0; i<n_size; i++,j++) {
 1962             if (j == n_size)
 1963                j = 0;
 1964 ';
 1965 pp_def('rotate',
 1966         GenericTypes => [ppdefs_all],
 1967         Doc => <<'EOD',
 1968 =for ref
 1969 
 1970 Shift vector elements along with wrap. Flows data back&forth.
 1971 =cut
 1972 EOD
 1973         Pars=>'x(n); indx shift(); [oca]y(n)',
 1974         DefaultFlow => 1,
 1975         TwoWay => 1,
 1976         Code=>$rotate_code.'
 1977             $y(n=>j) = $x(n=>i);
 1978         }',
 1979         BackCode=>$rotate_code.'
 1980             $x(n=>i) = $y(n=>j);
 1981         }
 1982         '
 1983 );
 1984 
 1985 # This is a bit tricky. Hope I haven't missed any cases.
 1986 pp_def(
 1987         'broadcastI',
 1988         Doc => <<'EOD',
 1989 =for ref
 1990 
 1991 internal
 1992 
 1993 Put some dimensions to a broadcastid.
 1994 
 1995 =for example
 1996 
 1997  $y = $x->broadcastI(0,1,5); # broadcast over dims 1,5 in id 1
 1998 
 1999 =cut
 2000 
 2001 EOD
 2002         P2Child => 1,
 2003         TwoWay => 1,
 2004         AffinePriv => 1,
 2005         OtherPars => "PDL_Indx id; PDL_Indx whichdims[]",
 2006         Comp => 'PDL_Indx nrealwhichdims',
 2007         MakeComp => '
 2008                 PDL_Indx i,j;
 2009                 $COMP(nrealwhichdims) = 0;
 2010                 for(i=0; i<$COMP(whichdims_count); i++) {
 2011                         for(j=i+1; j<$COMP(whichdims_count); j++)
 2012                                 if($COMP(whichdims[i]) == $COMP(whichdims[j]) &&
 2013                                    $COMP(whichdims[i]) != -1) {
 2014                                 $CROAK("duplicate arg %"IND_FLAG" %"IND_FLAG" %"IND_FLAG"",
 2015                                         i,j,$COMP(whichdims[i]));
 2016                         }
 2017                         if($COMP(whichdims)[i] != -1) {
 2018                                 $COMP(nrealwhichdims) ++;
 2019                         }
 2020                 }
 2021                 ',
 2022         RedoDims => '
 2023                 PDL_Indx nthc,i,j,flag;
 2024                 $SETNDIMS($PDL(PARENT)->ndims);
 2025                 $DOPRIVALLOC();
 2026                 $PRIV(offs) = 0;
 2027                 nthc=0;
 2028                 for(i=0; i<$PDL(PARENT)->ndims; i++) {
 2029                         flag=0;
 2030                         if($PDL(PARENT)->nbroadcastids > $COMP(id) && $COMP(id) >= 0 &&
 2031                            i == $PDL(PARENT)->broadcastids[$COMP(id])) {
 2032                            nthc += $COMP(whichdims_count);
 2033                         }
 2034                         for(j=0; j<$COMP(whichdims_count); j++) {
 2035                                 if($COMP(whichdims[j] == i)) {flag=1; break;}
 2036                         }
 2037                         if(flag) {
 2038                                 continue;
 2039                         }
 2040                         $PDL(CHILD)->dims[nthc] = $PDL(PARENT)->dims[i];
 2041                         $PRIV(incs[nthc]) = $PDL(PARENT)->dimincs[i];
 2042                         nthc++;
 2043                 }
 2044                 for(i=0; i<$COMP(whichdims_count); i++) {
 2045                         PDL_Indx cdim,pdim;
 2046                         cdim = i +
 2047                          ($PDL(PARENT)->nbroadcastids > $COMP(id) && $COMP(id) >= 0?
 2048                           $PDL(PARENT)->broadcastids[$COMP(id)] : $PDL(PARENT)->ndims)
 2049                           - $COMP(nrealwhichdims);
 2050                         pdim = $COMP(whichdims[i]);
 2051                         if(pdim == -1) {
 2052                                 $PDL(CHILD)->dims[cdim] = 1;
 2053                                 $PRIV(incs[cdim]) = 0;
 2054                         } else {
 2055                                 $PDL(CHILD)->dims[cdim] = $PDL(PARENT)->dims[pdim];
 2056                                 $PRIV(incs[cdim]) = $PDL(PARENT)->dimincs[pdim];
 2057                         }
 2058                 }
 2059                 $SETDIMS();
 2060                 PDL_RETERROR(PDL_err, PDL->reallocbroadcastids($PDL(CHILD),
 2061                                 PDLMAX($COMP(id)+1, $PDL(PARENT)->nbroadcastids)));
 2062                 for(i=0; i<$PDL(CHILD)->nbroadcastids-1; i++) {
 2063                         $PDL(CHILD)->broadcastids[i] =
 2064                          ($PDL(PARENT)->nbroadcastids > i ?
 2065                           $PDL(PARENT)->broadcastids[i] : $PDL(PARENT)->ndims) +
 2066                          (i <= $COMP(id) ? - $COMP(nrealwhichdims) :
 2067                           $COMP(whichdims_count) - $COMP(nrealwhichdims));
 2068                 }
 2069                 $PDL(CHILD)->broadcastids[$PDL(CHILD)->nbroadcastids-1] = $PDL(CHILD)->ndims;
 2070                 ',
 2071 );
 2072 
 2073 pp_def(
 2074         'unbroadcast',
 2075         Doc => <<'EOD',
 2076 =for ref
 2077 
 2078 All broadcasted dimensions are made real again.
 2079 
 2080 See [TBD Doc] for details and examples.
 2081 =cut
 2082 EOD
 2083         P2Child => 1,
 2084         TwoWay => 1,
 2085         AffinePriv => 1,
 2086         OtherPars => 'PDL_Indx atind;',
 2087         RedoDims => '
 2088                 PDL_Indx i;
 2089                 $SETNDIMS($PDL(PARENT)->ndims);
 2090                 $DOPRIVALLOC();
 2091                 $PRIV(offs) = 0;
 2092                 for(i=0; i<$PDL(PARENT)->ndims; i++) {
 2093                         PDL_Indx corc;
 2094                         if(i<$COMP(atind)) {
 2095                                 corc = i;
 2096                         } else if(i < $PDL(PARENT)->broadcastids[0]) {
 2097                                 corc = i + $PDL(PARENT)->ndims-$PDL(PARENT)->broadcastids[0];
 2098                         } else {
 2099                                 corc = i - $PDL(PARENT)->broadcastids[0] + $COMP(atind);
 2100                         }
 2101                         $PDL(CHILD)->dims[corc] = $PDL(PARENT)->dims[i];
 2102                         $PRIV(incs[corc]) = $PDL(PARENT)->dimincs[i];
 2103                 }
 2104                 $SETDIMS();
 2105         ',
 2106 );
 2107 
 2108 
 2109 pp_add_exported('', 'dice dice_axis');
 2110 pp_addpm(<<'EOD');
 2111 
 2112 =head2 dice
 2113 
 2114 =for ref
 2115 
 2116 Dice rows/columns/planes out of a PDL using indexes for
 2117 each dimension.
 2118 
 2119 This function can be used to extract irregular subsets
 2120 along many dimension of a PDL, e.g. only certain rows in an image,
 2121 or planes in a cube. This can of course be done with
 2122 the usual dimension tricks but this saves having to
 2123 figure it out each time!
 2124 
 2125 This method is similar in functionality to the L</slice>
 2126 method, but L</slice> requires that contiguous ranges or ranges
 2127 with constant offset be extracted. ( i.e. L</slice> requires
 2128 ranges of the form C<1,2,3,4,5> or C<2,4,6,8,10>). Because of this
 2129 restriction, L</slice> is more memory efficient and slightly faster
 2130 than dice
 2131 
 2132 =for usage
 2133 
 2134  $slice = $data->dice([0,2,6],[2,1,6]); # Dicing a 2-D array
 2135 
 2136 The arguments to dice are arrays (or 1D PDLs) for each dimension
 2137 in the PDL. These arrays are used as indexes to which rows/columns/cubes,etc
 2138 to dice-out (or extract) from the C<$data> PDL.
 2139 
 2140 Use C<X> to select all indices along a given dimension (compare also
 2141 L<mslice|PDL::Core/mslice>). As usual (in slicing methods) trailing
 2142 dimensions can be omitted implying C<X>'es for those.
 2143 
 2144 =for example
 2145 
 2146  pdl> $x = sequence(10,4)
 2147  pdl> p $x
 2148  [
 2149   [ 0  1  2  3  4  5  6  7  8  9]
 2150   [10 11 12 13 14 15 16 17 18 19]
 2151   [20 21 22 23 24 25 26 27 28 29]
 2152   [30 31 32 33 34 35 36 37 38 39]
 2153  ]
 2154  pdl> p $x->dice([1,2],[0,3]) # Select columns 1,2 and rows 0,3
 2155  [
 2156   [ 1  2]
 2157   [31 32]
 2158  ]
 2159  pdl> p $x->dice(X,[0,3])
 2160  [
 2161   [ 0  1  2  3  4  5  6  7  8  9]
 2162   [30 31 32 33 34 35 36 37 38 39]
 2163  ]
 2164  pdl> p $x->dice([0,2,5])
 2165  [
 2166   [ 0  2  5]
 2167   [10 12 15]
 2168   [20 22 25]
 2169   [30 32 35]
 2170  ]
 2171 
 2172 As this is an index function, any modifications to the
 2173 slice will change the parent (use the C<.=> operator).
 2174 
 2175 =cut
 2176 
 2177 sub PDL::dice {
 2178 
 2179         my $self = shift;
 2180         my @dim_indexes = @_;  # array of dimension indexes
 2181 
 2182         # Check that the number of dim indexes <=
 2183         #    number of dimensions in the PDL
 2184         my $no_indexes = scalar(@dim_indexes);
 2185         my $noDims = $self->getndims;
 2186         barf("PDL::dice: Number of index arrays ($no_indexes) not equal to the dimensions of the PDL ($noDims")
 2187                          if $no_indexes > $noDims;
 2188         my $index;
 2189         my $pdlIndex;
 2190         my $outputPDL=$self;
 2191         my $indexNo = 0;
 2192 
 2193         # Go thru each index array and dice the input PDL:
 2194         foreach $index(@dim_indexes){
 2195                 $outputPDL = $outputPDL->dice_axis($indexNo,$index)
 2196                         unless !ref $index && $index eq 'X';
 2197 
 2198                 $indexNo++;
 2199         }
 2200 
 2201         return $outputPDL;
 2202 }
 2203 *dice = \&PDL::dice;
 2204 
 2205 =head2 dice_axis
 2206 
 2207 =for ref
 2208 
 2209 Dice rows/columns/planes from a single PDL axis (dimension)
 2210 using index along a specified axis
 2211 
 2212 This function can be used to extract irregular subsets
 2213 along any dimension, e.g. only certain rows in an image,
 2214 or planes in a cube. This can of course be done with
 2215 the usual dimension tricks but this saves having to
 2216 figure it out each time!
 2217 
 2218 =for usage
 2219 
 2220  $slice = $data->dice_axis($axis,$index);
 2221 
 2222 =for example
 2223 
 2224  pdl> $x = sequence(10,4)
 2225  pdl> $idx = pdl(1,2)
 2226  pdl> p $x->dice_axis(0,$idx) # Select columns
 2227  [
 2228   [ 1  2]
 2229   [11 12]
 2230   [21 22]
 2231   [31 32]
 2232  ]
 2233  pdl> $t = $x->dice_axis(1,$idx) # Select rows
 2234  pdl> $t.=0
 2235  pdl> p $x
 2236  [
 2237   [ 0  1  2  3  4  5  6  7  8  9]
 2238   [ 0  0  0  0  0  0  0  0  0  0]
 2239   [ 0  0  0  0  0  0  0  0  0  0]
 2240   [30 31 32 33 34 35 36 37 38 39]
 2241  ]
 2242 
 2243 The trick to using this is that the index selects
 2244 elements along the dimensions specified, so if you
 2245 have a 2D image C<axis=0> will select certain C<X> values
 2246 - i.e. extract columns
 2247 
 2248 As this is an index function, any modifications to the
 2249 slice will change the parent.
 2250 
 2251 =cut
 2252 
 2253 sub PDL::dice_axis {
 2254   my($self,$axis,$idx) = @_;
 2255   my $ix = PDL->topdl($idx);
 2256   barf("dice_axis: index must be <=1D") if $ix->getndims > 1;
 2257   return $self->mv($axis,0)->index1d($ix)->mv(0,$axis);
 2258 }
 2259 *dice_axis = \&PDL::dice_axis;
 2260 EOD
 2261 
 2262 ##############################
 2263 # 'slice' is now implemented as a small Perl wrapper around
 2264 # a PP call.  This permits unification of the former slice,
 2265 # and dice into a single call.  At the moment, dicing
 2266 # is implemented a bit kludgily (it is detected in the Perl
 2267 # front-end), but it is serviceable.
 2268 #  --CED 12-Sep-2013
 2269 pp_def(
 2270         'slice',
 2271         Doc => <<'EOD-slice',
 2272 =for usage
 2273 
 2274   $slice = $data->slice([2,3],'x',[2,2,0],"-1:1:-1", "*3");
 2275 
 2276 =for ref
 2277 
 2278 Extract rectangular slices of an ndarray, from a string specifier,
 2279 an array ref specifier, or a combination.
 2280 
 2281 C<slice> is the main method for extracting regions of PDLs and
 2282 manipulating their dimensionality.  You can call it directly or
 2283 via the L<NiceSlice|PDL::NiceSlice> source prefilter that extends
 2284 Perl syntax to include array slicing.
 2285 
 2286 C<slice> can extract regions along each dimension of a source PDL,
 2287 subsample or reverse those regions, dice each dimension by selecting a
 2288 list of locations along it, or basic PDL indexing routine.  The
 2289 selected subfield remains connected to the original PDL via dataflow.
 2290 In most cases this neither allocates more memory nor slows down
 2291 subsequent operations on either of the two connected PDLs.
 2292 
 2293 You pass in a list of arguments.  Each term in the list controls
 2294 the disposition of one axis of the source PDL and/or returned PDL.
 2295 Each term can be a string-format cut specifier, a list ref that
 2296 gives the same information without recourse to string manipulation,
 2297 or a PDL with up to 1 dimension giving indices along that axis that
 2298 should be selected.
 2299 
 2300 If you want to pass in a single string specifier for the entire
 2301 operation, you can pass in a comma-delimited list as the first
 2302 argument.  C<slice> detects this condition and splits the string
 2303 into a regular argument list.  This calling style is fully
 2304 backwards compatible with C<slice> calls from before PDL 2.006.
 2305 
 2306 B<STRING SYNTAX>
 2307 
 2308 If a particular argument to C<slice> is a string, it is parsed as a
 2309 selection, an affine slice, or a dummy dimension depending on the
 2310 form.  Leading or trailing whitespace in any part of each specifier is
 2311 ignored (though it is not ignored within numbers).
 2312 
 2313 =over 3
 2314 
 2315 =item C<< '' >>, C<< : >>, or C<< X >> -- keep
 2316 
 2317 The empty string, C<:>, or C<X> cause the entire corresponding
 2318 dimension to be kept unchanged.
 2319 
 2320 
 2321 =item C<< <n> >> -- selection
 2322 
 2323 A single number alone causes a single index to be selected from the
 2324 corresponding dimension.  The dimension is kept (and reduced to size
 2325 1) in the output.
 2326 
 2327 =item C<< (<n>) >> -- selection and collapse
 2328 
 2329 A single number in parenthesis causes a single index to be selected
 2330 from the corresponding dimension.  The dimension is discarded
 2331 (completely eliminated) in the output.
 2332 
 2333 =item C<< <n>:<m> >> -- select an inclusive range
 2334 
 2335 Two numbers separated by a colon selects a range of values from the
 2336 corresponding axis, e.g. C<< 3:4 >> selects elements 3 and 4 along the
 2337 corresponding axis, and reduces that axis to size 2 in the output.
 2338 Both numbers are regularized so that you can address the last element
 2339 of the axis with an index of C< -1 >.  If, after regularization, the
 2340 two numbers are the same, then exactly one element gets selected (just
 2341 like the C<< <n> >> case).  If, after regulariation, the second number
 2342 is lower than the first, then the resulting slice counts down rather
 2343 than up -- e.g. C<-1:0> will return the entire axis, in reversed
 2344 order.
 2345 
 2346 =item C<< <n>:<m>:<s> >> -- select a range with explicit step
 2347 
 2348 If you include a third parameter, it is the stride of the extracted
 2349 range.  For example, C<< 0:-1:2 >> will sample every other element
 2350 across the complete dimension.  Specifying a stride of 1 prevents
 2351 autoreversal -- so to ensure that your slice is *always* forward
 2352 you can specify, e.g., C<< 2:$n:1 >>.  In that case, an "impossible"
 2353 slice gets an Empty PDL (with 0 elements along the corresponding
 2354 dimension), so you can generate an Empty PDL with a slice of the
 2355 form C<< 2:1:1 >>.
 2356 
 2357 =item C<< *<n> >> -- insert a dummy dimension
 2358 
 2359 Dummy dimensions aren't present in the original source and are
 2360 "mocked up" to match dimensional slots, by repeating the data
 2361 in the original PDL some number of times.  An asterisk followed
 2362 by a number produces a dummy dimension in the output, for
 2363 example C<< *2 >> will generate a dimension of size 2 at
 2364 the corresponding location in the output dim list.  Omitting
 2365 the number (and using just an asterisk) inserts a dummy dimension
 2366 of size 1.
 2367 
 2368 =back
 2369 
 2370 B<ARRAY REF SYNTAX>
 2371 
 2372 If you feed in an ARRAY ref as a slice term, then it can have
 2373 0-3 elements.  The first element is the start of the slice along
 2374 the corresponding dim; the second is the end; and the third is
 2375 the stepsize.  Different combinations of inputs give the same
 2376 flexibility as the string syntax.
 2377 
 2378 =over 3
 2379 
 2380 =item C<< [] >> - keep dim intact
 2381 
 2382 An empty ARRAY ref keeps the entire corresponding dim
 2383 
 2384 =item C<< [ 'X' ] >> - keep dim intact
 2385 
 2386 =item C<< [ '*',$n ] >> - generate a dummy dim of size $n
 2387 
 2388 If $n is missing, you get a dummy dim of size 1.
 2389 
 2390 =item C<< [ $dex, , 0 ] >> - collapse and discard dim
 2391 
 2392 C<$dex> must be a single value.  It is used to index
 2393 the source, and the corresponding dimension is discarded.
 2394 
 2395 =item C<< [ $start, $end ] >> - collect inclusive slice
 2396 
 2397 In the simple two-number case, you get a slice that runs
 2398 up or down (as appropriate) to connect $start and $end.
 2399 
 2400 =item C<< [ $start, $end, $inc ] >> - collect inclusive slice
 2401 
 2402 The three-number case works exactly like the three-number
 2403 string case above.
 2404 
 2405 =back
 2406 
 2407 B<PDL args for dicing>
 2408 
 2409 If you pass in a 0- or 1-D PDL as a slicing argument, the
 2410 corresponding dimension is "diced" -- you get one position
 2411 along the corresponding dim, per element of the indexing PDL,
 2412 e.g. C<< $x->slice( pdl(3,4,9)) >> gives you elements 3, 4, and
 2413 9 along the 0 dim of C<< $x >>.
 2414 
 2415 Because dicing is not an affine transformation, it is slower than
 2416 direct slicing even though the syntax is convenient.
 2417 
 2418 =for example
 2419 
 2420  $x->slice('1:3');  #  return the second to fourth elements of $x
 2421  $x->slice('3:1');  #  reverse the above
 2422  $x->slice('-2:1'); #  return last-but-one to second elements of $x
 2423 
 2424  $x->slice([1,3]);  # Same as above three calls, but using array ref syntax
 2425  $x->slice([3,1]);
 2426  $x->slice([-2,1]);
 2427 EOD-slice
 2428         PMCode => <<'EOD-slice',
 2429 sub PDL::slice {
 2430     my ($source, @others) = @_;
 2431     for my $i(0..$#others) {
 2432       my $idx = $others[$i];
 2433       if (ref $idx eq 'ARRAY') {
 2434         my @arr = map UNIVERSAL::isa($_, 'PDL') ? $_->flat->at(0) : $_, @{$others[$i]};
 2435         $others[$i] = \@arr;
 2436         next;
 2437       }
 2438       next if !( blessed($idx) && $idx->isa('PDL') );
 2439       # Deal with dicing.  This is lame and slow compared to the
 2440       # faster slicing, but works okay.  We loop over each argument,
 2441       # and if it's a PDL we dispatch it in the most straightforward
 2442       # way.  Single-element and zero-element PDLs are trivial and get
 2443       # converted into slices for faster handling later.
 2444       barf("slice: dicing parameters must be at most 1D (arg $i)\n")
 2445         if $idx->ndims > 1;
 2446       my $nlm = $idx->nelem;
 2447       if($nlm > 1) {
 2448          #### More than one element - we have to dice (darn it).
 2449          my $n = $source->getndims;
 2450          $source = $source->mv($i,0)->index1d($idx)->mv(0,$i);
 2451          $others[$i] = '';
 2452       }
 2453       elsif($nlm) {
 2454          #### One element - convert to a regular slice.
 2455          $others[$i] = $idx->flat->at(0);
 2456       }
 2457       else {
 2458          #### Zero elements -- force an extended empty.
 2459          $others[$i] = "1:0:1";
 2460       }
 2461     }
 2462     PDL::_slice_int($source,my $o=$source->initialize,\@others);
 2463     $o;
 2464 }
 2465 EOD-slice
 2466         P2Child => 1,
 2467         OtherPars => 'pdl_slice_args *arglist;',
 2468 #
 2469 # Comp stash definitions:
 2470 #  nargs - number of args in original call
 2471 #  odim[]   - maps argno to output dim (or -1 for squished dims)
 2472 #  idim[]   - maps argno to input dim  (or -1 for dummy dims)
 2473 #  odim_top - one more than the highest odim encountered
 2474 #  idim_top - one more than the highest idim encountered
 2475 #  start[]  - maps argno to start index of slice range (inclusive)
 2476 #  inc[]    - maps argno to increment of slice range
 2477 #  end[]    - maps argno to end index of slice range (inclusive)
 2478 #
 2479         Comp => 'PDL_Indx nargs;
 2480                  PDL_Indx odim[$COMP(nargs)];
 2481                  PDL_Indx idim[$COMP(nargs)];
 2482                  PDL_Indx idim_top;
 2483                  PDL_Indx odim_top;
 2484                  PDL_Indx start[$COMP(nargs)];
 2485                  PDL_Indx inc[$COMP(nargs)];
 2486                  PDL_Indx end[$COMP(nargs)];
 2487                  ',
 2488         AffinePriv => 1,
 2489         TwoWay => 1,
 2490         MakeComp => <<'SLICE-MC'
 2491                 PDL_Indx nargs = 0;
 2492                 pdl_slice_args *argsptr = arglist;
 2493                 while (argsptr) nargs++, argsptr = argsptr->next;
 2494                 $COMP(nargs) = nargs;
 2495                 $DOCOMPALLOC();
 2496                 PDL_Indx i, idim, odim, imax;
 2497                 argsptr = arglist;
 2498                 for(odim=idim=i=0; i<nargs; i++) {
 2499                   /* Copy parsed values into the limits */
 2500                   $COMP(start)[i] = argsptr->start;
 2501                   $COMP(end)[i]   = argsptr->end;
 2502                   $COMP(inc)[i]   = argsptr->inc;
 2503                   /* Deal with dimensions */
 2504                   $COMP(odim)[i]  = argsptr->squish ? -1 : odim++;
 2505                   $COMP(idim)[i]  = argsptr->dummy ? -1 : idim++;
 2506                   argsptr = argsptr->next;
 2507                 } /* end of arg-parsing loop */
 2508                 $COMP(idim_top) = idim;
 2509                 $COMP(odim_top) = odim;
 2510 SLICE-MC
 2511            ,
 2512            RedoDims => q{
 2513                   PDL_Indx i;
 2514           PDL_Indx PDIMS;
 2515                   PDL_Indx o_ndims_extra = PDLMAX(0, $PDL(PARENT)->ndims - $COMP(idim_top));
 2516 
 2517                   /* slurped dims from the arg parsing, plus any extra broadcast dims */
 2518                   $SETNDIMS( $COMP(odim_top) + o_ndims_extra );
 2519                   $DOPRIVALLOC();
 2520                   $PRIV(offs) = 0;  /* Offset vector to start of slice */
 2521 
 2522                   for(i=0; i<$COMP(nargs); i++) {
 2523                       /** Belt-and-suspenders **/
 2524                       if( ($COMP(idim[i]) < 0)  && ($COMP(odim[i]) < 0)  ) {
 2525                         PDL->changed($PDL(CHILD), PDL_PARENTDIMSCHANGED, 0);
 2526                         $CROAK("Hmmm, both dummy and squished -- this can never happen.  I quit.");
 2527                       }
 2528 
 2529                       /** First handle dummy dims since there's no input from the parent **/
 2530                       if( $COMP(idim[i]) < 0 ) {
 2531                          /* dummy dim - offset or diminc. */
 2532                          $PDL(CHILD)->dims[ $COMP(odim[i]) ] = $COMP(end[i]) - $COMP(start[i]) + 1;
 2533                          $PRIV( incs[ $COMP(odim[i]) ] ) = 0;
 2534                       } else {
 2535                         /** This is not a dummy dim -- deal with a regular slice along it.     **/
 2536                         /** Get parent dim size for this idim, and/or allow permissive slicing **/
 2537             PDL_Indx pdsize = $COMP(idim[i]) < $PDL(PARENT)->ndims ?
 2538                 $PDL(PARENT)->dims[$COMP(idim)[i]] : 1;
 2539             PDL_Indx start = $COMP(start)[i];
 2540             PDL_Indx end = $COMP(end)[i];
 2541             if(
 2542               /** Trap special case: full slices of an empty dim are empty **/
 2543                           (pdsize==0 && start==0 && end==-1 && $COMP(inc[i]) == 0)
 2544                           ||
 2545                           /* the values given when PDL::slice gets empty ndarray for index */
 2546                           (start==1 && end==0 && $COMP(inc[i]) == 1)
 2547                         ) {
 2548               $PDL(CHILD)->dims[$COMP(odim)[i]] = 0;
 2549               $PRIV(incs)[$COMP(odim)[i]] = 0;
 2550             } else {
 2551               /** Regularize and bounds-check the start location **/
 2552               if(start < 0)
 2553                 start += pdsize;
 2554               if( start < 0 || start >= pdsize ) {
 2555                 PDL->changed($PDL(CHILD), PDL_PARENTDIMSCHANGED, 0);
 2556                 if(i >= $PDL(PARENT)->ndims) {
 2557                   $CROAK("slice has too many dims (indexes dim %"IND_FLAG"; highest is %"IND_FLAG")",i,$PDL(PARENT)->ndims-1);
 2558                 } else {
 2559               $CROAK("slice starts out of bounds in pos %"IND_FLAG" (start is %"IND_FLAG"; source dim %"IND_FLAG" runs 0 to %"IND_FLAG")",i,start,$COMP(idim[i]),pdsize-1);
 2560                 }
 2561               }
 2562               if( $COMP(odim[i]) < 0) {
 2563                 /* squished dim - just update the offset. */
 2564                 /* start is always defined and regularized if we are here here, since */
 2565                 /* both idim[i] and odim[i] can't be <0 */
 2566                 $PRIV(offs) += start * $PDL(PARENT)->dimincs[ $COMP(idim[i]) ];
 2567               } else /* normal operation */ {
 2568                 /** Regularize and bounds-check the end location **/
 2569                 if(end<0) end += pdsize;
 2570                 if( end < 0 || end >= pdsize ) {
 2571                   PDL->changed($PDL(CHILD), PDL_PARENTDIMSCHANGED, 0);
 2572                   $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);
 2573                 }
 2574                 /* regularize inc */
 2575                 PDL_Indx inc = $COMP(inc)[i];
 2576                 if(!inc)
 2577                   inc = (start <= end) ? 1 : -1;
 2578                 $PDL(CHILD)->dims[ $COMP(odim)[i] ] = PDLMAX(0, (end - start + inc) / inc);
 2579                 $PRIV(  incs )[ $COMP(odim)[i] ] = $PDL(PARENT)->dimincs[ $COMP(idim)[i] ] * inc;
 2580                 $PRIV(offs) += start * $PDL(PARENT)->dimincs[ $COMP(idim)[i] ];
 2581               } /* end of normal slice case */
 2582             } /* end of trapped strange slice case */
 2583                       } /* end of non-dummy slice case */
 2584                   } /* end of nargs loop */
 2585 
 2586                   /* Now fill in broadcast dimensions as needed.  idim and odim persist from the parsing loop */
 2587                   /* up above. */
 2588                   for(i=0; i<o_ndims_extra; i++) {
 2589                     $PDL(CHILD)->dims[ $COMP(odim_top) + i ] = $PDL(PARENT)->dims[ $COMP(idim_top) + i ];
 2590                     $PRIV(incs)[ $COMP(odim_top) + i ] = $PDL(PARENT)->dimincs[ $COMP(idim_top) + i ];
 2591                   }
 2592 
 2593                 $SETDIMS();
 2594 
 2595         } # end of RedoDims for slice
 2596 );
 2597 
 2598 
 2599 pp_addpm({At => 'Bot'},<< 'EOD');
 2600 
 2601 =head1 BUGS
 2602 
 2603 For the moment, you can't slice one of the zero-length dims of an
 2604 empty ndarray.  It is not clear how to implement this in a way that makes
 2605 sense.
 2606 
 2607 Many types of index errors are reported far from the indexing
 2608 operation that caused them.  This is caused by the underlying architecture:
 2609 slice() sets up a mapping between variables, but that mapping isn't
 2610 tested for correctness until it is used (potentially much later).
 2611 
 2612 =head1 AUTHOR
 2613 
 2614 Copyright (C) 1997 Tuomas J. Lukka.  Contributions by
 2615 Craig DeForest, deforest@boulder.swri.edu.
 2616 Documentation contributions by David Mertens.
 2617 All rights reserved. There is no warranty. You are allowed
 2618 to redistribute this software / documentation under certain
 2619 conditions. For details, see the file COPYING in the PDL
 2620 distribution. If this file is separated from the PDL distribution,
 2621 the copyright notice should be included in the file.
 2622 
 2623 =cut
 2624 
 2625 EOD
 2626 
 2627 
 2628 pp_done();