"Fossies" - the Fresh Open Source Software Archive

Member "PDL-2.025/Basic/Slices/slices.pd" (19 Nov 2020, 120135 Bytes) of package /linux/misc/PDL-2.025.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.024_vs_2.025.

    1 pp_addpm({At => 'Top'},<< 'EOD');
    2 
    3 =head1 NAME
    4 
    5 PDL::Slices -- Indexing, slicing, and dicing
    6 
    7 =head1 SYNOPSIS
    8 
    9   use PDL;
   10   $x = ones(3,3);
   11   $y = $x->slice('-1:0,(1)');
   12   $c = $x->dummy(2);
   13 
   14 
   15 =head1 DESCRIPTION
   16 
   17 This package provides many of the powerful PerlDL core index
   18 manipulation routines.  These routines mostly allow two-way data flow,
   19 so you can modify your data in the most convenient representation.
   20 For example, you can make a 1000x1000 unit matrix with
   21 
   22  $x = zeroes(1000,1000);
   23  $x->diagonal(0,1) ++;
   24 
   25 which is quite efficient. See L<PDL::Indexing> and L<PDL::Tips> for
   26 more examples.
   27 
   28 Slicing is so central to the PDL language that a special compile-time
   29 syntax has been introduced to handle it compactly; see L<PDL::NiceSlice>
   30 for details.
   31 
   32 PDL indexing and slicing functions usually include two-way data flow,
   33 so that you can separate the actions of reshaping your data structures
   34 and modifying the data themselves.  Two special methods, L<copy|copy> and
   35 L<sever|sever>, help you control the data flow connection between related
   36 variables.
   37 
   38  $y = $x->slice("1:3"); # Slice maintains a link between $x and $y.
   39  $y += 5;               # $x is changed!
   40 
   41 If you want to force a physical copy and no data flow, you can copy or
   42 sever the slice expression:
   43 
   44  $y = $x->slice("1:3")->copy;
   45  $y += 5;               # $x is not changed.
   46 
   47  $y = $x->slice("1:3")->sever;
   48  $y += 5;               # $x is not changed.
   49 
   50 The difference between C<sever> and C<copy> is that sever acts on (and
   51 returns) its argument, while copy produces a disconnected copy.  If you
   52 say
   53 
   54  $y = $x->slice("1:3");
   55  $c = $y->sever;
   56 
   57 then the variables C<$y> and C<$c> point to the same object but with
   58 C<-E<gt>copy> they would not.
   59 
   60 =cut
   61 
   62 use PDL::Core ':Internal';
   63 use Scalar::Util 'blessed';
   64 
   65 EOD
   66 
   67 =head1
   68 FUNCTIONS
   69 =cut
   70 
   71 # $::PP_VERBOSE=1;
   72 
   73 pp_addhdr(<<'EOH');
   74 
   75 #ifdef _MSC_VER
   76 #if _MSC_VER < 1300
   77 #define strtoll strtol
   78 #else
   79 #define strtoll _strtoi64
   80 #endif
   81 #endif
   82 
   83 EOH
   84 
   85 pp_add_boot(
   86 "  PDL->readdata_affine = pdl_readdata_affineinternal;\n" .
   87 "     PDL->writebackdata_affine = pdl_writebackdata_affineinternal;\n"
   88 );
   89 
   90 ## Several routines use the 'Dims' and 'ParentInds'
   91 ## rules - these currently do nothing
   92 
   93 pp_def(
   94        'affineinternal',
   95        HandleBad => 1,
   96        AffinePriv => 1,
   97        DefaultFlow => 1,
   98        P2Child => 1,
   99        NoPdlThread => 1,
  100        ReadDataFuncName => "pdl_readdata_affineinternal",
  101        WriteBackDataFuncName => "pdl_writebackdata_affineinternal",
  102        MakeComp => '$CROAK("AFMC MUSTNT BE CALLED");',
  103        RedoDims => '$CROAK("AFRD MUSTNT BE CALLED");',
  104        EquivCPOffsCode => '
  105                 PDL_Indx i; PDL_Indx poffs=$PRIV(offs); int nd;
  106                 for(i=0; i<$CHILD_P(nvals); i++) {
  107                         $EQUIVCPOFFS(i,poffs);
  108                         for(nd=0; nd<$CHILD_P(ndims); nd++) {
  109                                 poffs += $PRIV(incs[nd]);
  110                                 if( (nd<$CHILD_P(ndims)-1 &&
  111                                      (i+1)%$CHILD_P(dimincs[nd+1])) ||
  112                                    nd == $CHILD_P(ndims)-1)
  113                                         break;
  114                                 poffs -= $PRIV(incs[nd]) *
  115                                         $CHILD_P(dims[nd]);
  116                         }
  117                 }',
  118        Doc => undef,    # 'internal',
  119 );
  120 
  121 =head2 s_identity
  122 =cut
  123 
  124 $doc = <<'DOC';
  125 =for ref
  126 
  127 Internal vaffine identity function.
  128 
  129 =cut
  130 
  131 DOC
  132 
  133 pp_def(
  134         's_identity',
  135         HandleBad => 1,
  136         P2Child => 1,
  137         NoPdlThread => 1,
  138         DefaultFlow => 1,
  139         OtherPars => '',
  140         Reversible => 1,
  141         Dims => '$COPYDIMS();',
  142         ParentInds => '$COPYINDS();',
  143         Identity => 1,
  144         Doc => $doc,
  145 );
  146 
  147 =head2 index, index1d, index2d
  148 
  149 =cut
  150 
  151 $doc = <<'EOD';
  152 =for ref
  153 
  154 C<index>, C<index1d>, and C<index2d> provide rudimentary index indirection.
  155 
  156 =for example
  157 
  158  $c = index($source,$ind);
  159  $c = index1d($source,$ind);
  160  $c = index2d($source2,$ind1,$ind2);
  161 
  162 use the C<$ind> variables as indices to look up values in C<$source>.
  163 The three routines thread slightly differently.
  164 
  165 =over 3
  166 
  167 =item * 
  168 
  169 C<index> uses direct threading for 1-D indexing across the 0 dim
  170 of C<$source>.  It can thread over source thread dims or index thread
  171 dims, but not (easily) both: If C<$source> has more than 1
  172 dimension and C<$ind> has more than 0 dimensions, they must agree in
  173 a threading sense.
  174 
  175 =item * 
  176 
  177 C<index1d> uses a single active dim in C<$ind> to produce a list of
  178 indexed values in the 0 dim of the output - it is useful for
  179 collapsing C<$source> by indexing with a single row of values along
  180 C<$source>'s 0 dimension.  The output has the same number of dims as
  181 C<$source>.  The 0 dim of the output has size 1 if C<$ind> is a
  182 scalar, and the same size as the 0 dim of C<$ind> if it is not. If
  183 C<$ind> and C<$source> both have more than 1 dim, then all dims higher
  184 than 0 must agree in a threading sense.
  185 
  186 =item * 
  187 
  188 C<index2d> works like C<index> but uses separate piddles for X and Y
  189 coordinates.  For more general N-dimensional indexing, see the
  190 L<PDL::NiceSlice|PDL::NiceSlice> syntax or L<PDL::Slices|PDL::Slices> (in particular C<slice>,
  191 C<indexND>, and C<range>).
  192 
  193 =back 
  194 
  195 These functions are two-way, i.e. after
  196 
  197  $c = $x->index(pdl[0,5,8]);
  198  $c .= pdl [0,2,4];
  199 
  200 the changes in C<$c> will flow back to C<$x>.
  201 
  202 C<index> provids simple threading:  multiple-dimensioned arrays are treated
  203 as collections of 1-D arrays, so that
  204 
  205  $x = xvals(10,10)+10*yvals(10,10);
  206  $y = $x->index(3);
  207  $c = $x->index(9-xvals(10));
  208 
  209 puts a single column from C<$x> into C<$y>, and puts a single element
  210 from each column of C<$x> into C<$c>.  If you want to extract multiple
  211 columns from an array in one operation, see L<dice|/dice> or
  212 L<indexND|/indexND>.
  213 
  214 =cut
  215 
  216 EOD
  217 
  218 my $index_init_good =
  219        'register PDL_Indx foo = $ind();
  220         if( foo<0 || foo>=$SIZE(n) ) {
  221            barf("PDL::index: invalid index %d (valid range 0..%d)",
  222                 foo,$SIZE(n)-1);
  223         }';
  224 my $index_init_bad =
  225        'register PDL_Indx foo = $ind();
  226         if( $ISBADVAR(foo,ind) || foo<0 || foo>=$SIZE(n) ) {
  227            barf("PDL::index: invalid index %d (valid range 0..%d)",
  228                 foo,$SIZE(n)-1);
  229         }';
  230 
  231 pp_def(
  232        'index',
  233        HandleBad => 1,
  234        DefaultFlow => 1,
  235        Reversible => 1,
  236        Pars => 'a(n); indx ind(); [oca] c();',
  237        Code =>
  238        $index_init_good . ' $c() = $a(n => foo);',
  239        BadCode =>
  240        $index_init_bad . ' $c() = $a(n => foo);',
  241        BackCode =>
  242        $index_init_good . ' $a(n => foo) = $c();',
  243        BadBackCode =>
  244        $index_init_bad . ' $a(n => foo) = $c();',
  245        Doc => $doc,
  246        BadDoc =>
  247        'index barfs if any of the index values are bad.',
  248        );
  249 
  250 pp_def(
  251        'index1d',
  252        HandleBad => 1,
  253        DefaultFlow => 1,
  254        Reversible => 1,
  255        Pars => 'a(n); indx ind(m); [oca] c(m);',
  256        Code =>
  257        q{
  258          PDL_Indx i;
  259          for(i=0;i<$SIZE(m);i++) {
  260                   PDL_Indx foo = $ind(m=>i);
  261                   if( foo<0 || foo >= $SIZE(n) ) {
  262                    barf("PDL::index1d: invalid index %d at pos %d (valid range 0..%d)",
  263                      foo, i, $SIZE(n)-1);
  264                   }
  265                   $c(m=>i) = $a(n=>foo);
  266          }
  267         },
  268         BadCode =>
  269         q{
  270          PDL_Indx i;
  271          for(i=0;i<$SIZE(m);i++) {
  272                  PDL_Indx foo = $ind(m=>i);
  273                  if( $ISBADVAR(foo, ind) ) {
  274                      $SETBAD(c(m=>i));
  275                  } else {
  276                    if( foo<0 || foo >= $SIZE(n) ) {
  277                      barf("PDL::index1d: invalid/bad index %d at pos %d (valid range 0..%d)",
  278                       foo, i, $SIZE(n)-1);
  279                    }
  280                    $c(m=>i) = $a(n=>foo);
  281                  }
  282          }
  283         },
  284        BackCode => q{
  285          PDL_Indx i;
  286          for(i=0;i<$SIZE(m);i++) {
  287                 PDL_Indx foo = $ind(m=>i);
  288                 if( foo<0 || foo >= $SIZE(n) ) {
  289                  barf("PDL::index1d: invalid index %d at pos %d (valid range 0..%d)",
  290                    foo, i, $SIZE(n)-1);
  291                 }
  292                 $a(n=>foo) = $c(m=>i);
  293          }
  294         },
  295        BadBackCode => q{
  296          PDL_Indx i;
  297          for(i=0;i<$SIZE(m);i++) {
  298                  PDL_Indx foo = $ind(m=>i);
  299                  if( $ISBADVAR(foo, ind) ) {
  300                    /* do nothing */
  301                  } else {
  302                    if( foo<0 || foo >= $SIZE(n) ) {
  303                      barf("PDL::index1d: invalid/bad index %d at pos %d (valid range 0..%d)",
  304                       foo, i, $SIZE(n)-1);
  305                    }
  306                    $a(n=>foo) = $c(m=>i);
  307                  }
  308          }
  309         },
  310        Doc => $doc,
  311        BadDoc =>
  312        'index1d propagates BAD index elements to the output variable.'
  313        );
  314 
  315 my $index2d_init_good =
  316        'register PDL_Indx fooa,foob;
  317         fooa = $inda();
  318         if( fooa<0 || fooa>=$SIZE(na) ) {
  319            barf("PDL::index: invalid x-index %d (valid range 0..%d)",
  320                 fooa,$SIZE(na)-1);
  321         }
  322         foob = $indb();
  323         if( foob<0 || foob>=$SIZE(nb) ) {
  324            barf("PDL::index: invalid y-index %d (valid range 0..%d)",
  325                 foob,$SIZE(nb)-1);
  326         }';
  327 my $index2d_init_bad =
  328        'register PDL_Indx fooa,foob;
  329         fooa = $inda();
  330         if( $ISBADVAR(fooa,inda) || fooa<0 || fooa>=$SIZE(na) ) {
  331            barf("PDL::index: invalid index 1");
  332         }
  333         foob = $indb();
  334         if( $ISBADVAR(foob,indb) || foob<0 || foob>=$SIZE(nb) ) {
  335            barf("PDL::index: invalid index 2");
  336         }';
  337 
  338 pp_def(
  339        'index2d',
  340        HandleBad => 1,
  341        DefaultFlow => 1,
  342        Reversible => 1,
  343        Pars => 'a(na,nb); indx inda(); indx indb(); [oca] c();',
  344        Code =>
  345        $index2d_init_good . ' $c() = $a(na => fooa, nb => foob);',
  346        BadCode =>
  347        $index2d_init_bad . '$c() = $a(na => fooa, nb => foob);',
  348        BackCode =>
  349        $index2d_init_good . ' $a(na => fooa, nb => foob) = $c();',
  350        BadBackCode =>
  351        $index2d_init_bad . '$a(na => fooa, nb => foob) = $c();',
  352        Doc => $doc,
  353        BadDoc =>
  354        'index2d barfs if either of the index values are bad.',
  355 );
  356 
  357 
  358 # indexND: CED 2-Aug-2002
  359 pp_add_exported('','indexND indexNDb');
  360 pp_addpm(<<'EOD-indexND');
  361 
  362 =head2 indexNDb
  363 
  364 =for ref
  365 
  366   Backwards-compatibility alias for indexND
  367 
  368 =head2 indexND
  369 
  370 =for ref
  371 
  372   Find selected elements in an N-D piddle, with optional boundary handling
  373 
  374 =for example
  375 
  376   $out = $source->indexND( $index, [$method] )
  377 
  378   $source = 10*xvals(10,10) + yvals(10,10);
  379   $index  = pdl([[2,3],[4,5]],[[6,7],[8,9]]);
  380   print $source->indexND( $index );
  381 
  382   [
  383    [23 45]
  384    [67 89]
  385   ]
  386 
  387 IndexND collapses C<$index> by lookup into C<$source>.  The
  388 0th dimension of C<$index> is treated as coordinates in C<$source>, and
  389 the return value has the same dimensions as the rest of C<$index>.
  390 The returned elements are looked up from C<$source>.  Dataflow
  391 works -- propagated assignment flows back into C<$source>.
  392 
  393 IndexND and IndexNDb were originally separate routines but they are both
  394 now implemented as a call to L<range|/range>, and have identical syntax to
  395 one another.
  396 
  397 =cut
  398 
  399 sub PDL::indexND {
  400         my($source,$index, $boundary) = @_;
  401         return PDL::range($source,$index,undef,$boundary);
  402 }
  403 
  404 *PDL::indexNDb = \&PDL::indexND;
  405 
  406 EOD-indexND
  407 
  408 pp_addpm(<<'EOD-range');
  409 
  410 sub PDL::range {
  411   my($source,$ind,$sz,$bound) = @_;
  412 
  413 # Convert to indx type up front (also handled in rangeb if necessary)
  414   my $index = (ref $ind && UNIVERSAL::isa($ind,'PDL') && $ind->type eq 'indx') ? $ind : indx($ind);
  415   my $size = defined($sz) ? PDL->pdl($sz) : undef;
  416 
  417 
  418   # Handle empty PDL case: return a properly constructed Empty.
  419   if($index->isempty) {
  420       my @sdims= $source->dims;
  421       splice(@sdims, 0, $index->dim(0) + ($index->dim(0)==0)); # added term is to treat Empty[0] like a single empty coordinate
  422       unshift(@sdims, $size->list) if(defined($size));
  423       return PDL->new_from_specification(0 x ($index->ndims-1), @sdims);
  424   }
  425 
  426 
  427   $index = $index->dummy(0,1) unless $index->ndims;
  428 
  429 
  430   # Pack boundary string if necessary
  431   if(defined $bound) {
  432     if(ref $bound eq 'ARRAY') {
  433       my ($s,$el);
  434       foreach $el(@$bound) {
  435         barf "Illegal boundary value '$el' in range"
  436           unless( $el =~ m/^([0123fFtTeEpPmM])/ );
  437         $s .= $1;
  438       }
  439       $bound = $s;
  440     }
  441     elsif($bound !~ m/^[0123ftepx]+$/  && $bound =~ m/^([0123ftepx])/i ) {
  442       $bound = $1;
  443     }
  444   }
  445 
  446   no warnings; # shut up about passing undef into rangeb
  447   $source->rangeb($index,$size,$bound);
  448 }
  449 EOD-range
  450 
  451 =head2 range
  452 
  453 =cut
  454 
  455 pp_def(
  456         'rangeb',
  457         OtherPars => 'SV *index; SV *size; SV *boundary',
  458         Doc => <<'EOD',
  459 =for ref
  460 
  461 Engine for L<range|/range>
  462 
  463 =for example
  464 
  465 Same calling convention as L<range|/range>, but you must supply all
  466 parameters.  C<rangeb> is marginally faster as it makes a direct PP call,
  467 avoiding the perl argument-parsing step.
  468 
  469 =cut
  470 
  471 =head2 range
  472 
  473 =for ref
  474 
  475 Extract selected chunks from a source piddle, with boundary conditions
  476 
  477 =for example
  478 
  479         $out = $source->range($index,[$size,[$boundary]])
  480 
  481 Returns elements or rectangular slices of the original piddle, indexed by
  482 the C<$index> piddle.  C<$source> is an N-dimensional piddle, and C<$index> is
  483 a piddle whose first dimension has size up to N.  Each row of C<$index> is
  484 treated as coordinates of a single value or chunk from C<$source>, specifying
  485 the location(s) to extract.
  486 
  487 If you specify a single index location, then range is essentially an expensive
  488 slice, with controllable boundary conditions.
  489 
  490 B<INPUTS>
  491 
  492 C<$index> and C<$size> can be piddles or array refs such as you would
  493 feed to L<zeroes|PDL::Core/zeroes> and its ilk.  If C<$index>'s 0th dimension
  494 has size higher than the number of dimensions in C<$source>, then
  495 C<$source> is treated as though it had trivial dummy dimensions of
  496 size 1, up to the required size to be indexed by C<$index> -- so if
  497 your source array is 1-D and your index array is a list of 3-vectors,
  498 you get two dummy dimensions of size 1 on the end of your source array.
  499 
  500 You can extract single elements or N-D rectangular ranges from C<$source>,
  501 by setting C<$size>.  If C<$size> is undef or zero, then you get a single
  502 sample for each row of C<$index>.  This behavior is similar to
  503 L<indexNDb|/indexNDb>, which is in fact implemented as a call to L<range|/range>.
  504 
  505 If C<$size> is positive then you get a range of values from C<$source> at
  506 each location, and the output has extra dimensions allocated for them.
  507 C<$size> can be a scalar, in which case it applies to all dimensions, or an
  508 N-vector, in which case each element is applied independently to the
  509 corresponding dimension in C<$source>.  See below for details.
  510 
  511 C<$boundary> is a number, string, or list ref indicating the type of
  512 boundary conditions to use when ranges reach the edge of C<$source>.  If you
  513 specify no boundary conditions the default is to forbid boundary violations
  514 on all axes.  If you specify exactly one boundary condition, it applies to
  515 all axes.  If you specify more (as elements of a list ref, or as a packed
  516 string, see below), then they apply to dimensions in the order in which they
  517 appear, and the last one applies to all subsequent dimensions.  (This is
  518 less difficult than it sounds; see the examples below).
  519 
  520 =over 3
  521 
  522 =item 0 (synonyms: 'f','forbid') B<(default)>
  523 
  524 Ranges are not allowed to cross the boundary of the original PDL.  Disallowed
  525 ranges throw an error.  The errors are thrown at evaluation time, not
  526 at the time of the range call (this is the same behavior as L<slice|/slice>).
  527 
  528 =item 1 (synonyms: 't','truncate')
  529 
  530 Values outside the original piddle get BAD if you've got bad value
  531 support compiled into your PDL and set the badflag for the source PDL;
  532 or 0 if you haven't (you must set the badflag if you want BADs for out
  533 of bound values, otherwise you get 0).  Reverse dataflow works OK for
  534 the portion of the child that is in-bounds.  The out-of-bounds part of
  535 the child is reset to (BAD|0) during each dataflow operation, but
  536 execution continues.
  537 
  538 =item 2 (synonyms: 'e','x','extend')
  539 
  540 Values that would be outside the original piddle point instead to the
  541 nearest allowed value within the piddle.  See the CAVEAT below on
  542 mappings that are not single valued.
  543 
  544 =item 3 (synonyms: 'p','periodic')
  545 
  546 Periodic boundary conditions apply: the numbers in $index are applied,
  547 strict-modulo the corresponding dimensions of $source.  This is equivalent to
  548 duplicating the $source piddle throughout N-D space.  See the CAVEAT below
  549 about mappings that are not single valued.
  550 
  551 =item 4 (synonyms: 'm','mirror')
  552 
  553 Mirror-reflection periodic boundary conditions apply.  See the CAVEAT
  554 below about mappings that are not single valued.
  555 
  556 =back
  557 
  558 The boundary condition identifiers all begin with unique characters, so
  559 you can feed in multiple boundary conditions as either a list ref or a
  560 packed string.  (The packed string is marginally faster to run).  For
  561 example, the four expressions [0,1], ['forbid','truncate'], ['f','t'],
  562 and 'ft' all specify that violating the boundary in the 0th dimension
  563 throws an error, and all other dimensions get truncated.
  564 
  565 If you feed in a single string, it is interpreted as a packed boundary
  566 array if all of its characters are valid boundary specifiers (e.g. 'pet'),
  567 but as a single word-style specifier if they are not (e.g. 'forbid').
  568 
  569 B<OUTPUT>
  570 
  571 The output threads over both C<$index> and C<$source>.  Because implicit
  572 threading can happen in a couple of ways, a little thought is needed.  The
  573 returned dimension list is stacked up like this:
  574 
  575    (index thread dims), (index dims (size)), (source thread dims)
  576 
  577 The first few dims of the output correspond to the extra dims of
  578 C<$index> (beyond the 0 dim). They allow you to pick out individual
  579 ranges from a large, threaded collection.
  580 
  581 The middle few dims of the output correspond to the size dims
  582 specified in C<$size>, and contain the range of values that is extracted
  583 at each location in C<$source>.  Every nonzero element of C<$size> is copied to
  584 the dimension list here, so that if you feed in (for example) C<$size
  585 = [2,0,1]> you get an index dim list of C<(2,1)>.
  586 
  587 The last few dims of the output correspond to extra dims of C<$source> beyond
  588 the number of dims indexed by C<$index>.  These dims act like ordinary
  589 thread dims, because adding more dims to C<$source> just tacks extra dims
  590 on the end of the output.  Each source thread dim ranges over the entire
  591 corresponding dim of C<$source>.
  592 
  593 B<Dataflow>: Dataflow is bidirectional.
  594 
  595 B<Examples>:
  596 Here are basic examples of C<range> operation, showing how to get
  597 ranges out of a small matrix.  The first few examples show extraction
  598 and selection of individual chunks.  The last example shows
  599 how to mark loci in the original matrix (using dataflow).
  600 
  601  pdl> $src = 10*xvals(10,5)+yvals(10,5)
  602  pdl> print $src->range([2,3])    # Cut out a single element
  603  23
  604  pdl> print $src->range([2,3],1)  # Cut out a single 1x1 block
  605  [
  606   [23]
  607  ]
  608  pdl> print $src->range([2,3], [2,1]) # Cut a 2x1 chunk
  609  [
  610   [23 33]
  611  ]
  612  pdl> print $src->range([[2,3]],[2,1]) # Trivial list of 1 chunk
  613  [
  614   [
  615    [23]
  616    [33]
  617   ]
  618  ]
  619  pdl> print $src->range([[2,3],[0,1]], [2,1])   # two 2x1 chunks
  620  [
  621   [
  622    [23  1]
  623    [33 11]
  624   ]
  625  ]
  626  pdl> # A 2x2 collection of 2x1 chunks
  627  pdl> print $src->range([[[1,1],[2,2]],[[2,3],[0,1]]],[2,1])
  628  [
  629   [
  630    [
  631     [11 22]
  632     [23  1]
  633    ]
  634    [
  635     [21 32]
  636     [33 11]
  637    ]
  638   ]
  639  ]
  640  pdl> $src = xvals(5,3)*10+yvals(5,3)
  641  pdl> print $src->range(3,1)  # Thread over y dimension in $src
  642  [
  643   [30]
  644   [31]
  645   [32]
  646  ]
  647 
  648  pdl> $src = zeroes(5,4);
  649  pdl> $src->range(pdl([2,3],[0,1]),pdl(2,1)) .= xvals(2,2,1) + 1
  650  pdl> print $src
  651  [
  652   [0 0 0 0 0]
  653   [2 2 0 0 0]
  654   [0 0 0 0 0]
  655   [0 0 1 1 0]
  656  ]
  657 
  658 B<CAVEAT>: It's quite possible to select multiple ranges that
  659 intersect.  In that case, modifying the ranges doesn't have a
  660 guaranteed result in the original PDL -- the result is an arbitrary
  661 choice among the valid values.  For some things that's OK; but for
  662 others it's not. In particular, this doesn't work:
  663 
  664     pdl> $photon_list = new PDL::RandVar->sample(500)->reshape(2,250)*10
  665     pdl> histogram = zeroes(10,10)
  666     pdl> histogram->range($photon_list,1)++;  #not what you wanted
  667 
  668 The reason is that if two photons land in the same bin, then that bin
  669 doesn't get incremented twice.  (That may get fixed in a later version...)
  670 
  671 B<PERMISSIVE RANGING>: If C<$index> has too many dimensions compared
  672 to C<$source>, then $source is treated as though it had dummy
  673 dimensions of size 1, up to the required number of dimensions.  These
  674 virtual dummy dimensions have the usual boundary conditions applied to
  675 them.
  676 
  677 If the 0 dimension of C<$index> is ludicrously large (if its size is
  678 more than 5 greater than the number of dims in the source PDL) then
  679 range will insist that you specify a size in every dimension, to make
  680 sure that you know what you're doing.  That catches a common error with
  681 range usage: confusing the initial dim (which is usually small) with another
  682 index dim (perhaps of size 1000).
  683 
  684 If the index variable is Empty, then range() always returns the Empty PDL.
  685 If the index variable is not Empty, indexing it always yields a boundary
  686 violation.  All non-barfing conditions are treated as truncation, since
  687 there are no actual data to return.
  688 
  689 B<EFFICIENCY>: Because C<range> isn't an affine transformation (it
  690 involves lookup into a list of N-D indices), it is somewhat
  691 memory-inefficient for long lists of ranges, and keeping dataflow open
  692 is much slower than for affine transformations (which don't have to copy
  693 data around).
  694 
  695 Doing operations on small subfields of a large range is inefficient
  696 because the engine must flow the entire range back into the original
  697 PDL with every atomic perl operation, even if you only touch a single element.
  698 One way to speed up such code is to sever your range, so that PDL
  699 doesn't have to copy the data with each operation, then copy the
  700 elements explicitly at the end of your loop.  Here's an example that
  701 labels each region in a range sequentially, using many small
  702 operations rather than a single xvals assignment:
  703 
  704   ### How to make a collection of small ops run fast with range...
  705   $x =  $data->range($index, $sizes, $bound)->sever;
  706   $aa = $data->range($index, $sizes, $bound);
  707   map { $x($_ - 1) .= $_; } (1..$x->nelem);    # Lots of little ops
  708   $aa .= $x;
  709 
  710 C<range> is a perl front-end to a PP function, C<rangeb>.  Calling
  711 C<rangeb> is marginally faster but requires that you include all arguments.
  712 
  713 DEVEL NOTES
  714 
  715 * index thread dimensions are effectively clumped internally.  This
  716 makes it easier to loop over the index array but a little more brain-bending
  717 to tease out the algorithm.
  718 
  719 =cut
  720 
  721 EOD
  722         HandleBad => 1,
  723         DefaultFlow => 1,
  724         Reversible => 1,
  725         P2Child => 1,
  726         NoPdlThread => 1,
  727 
  728 #
  729 # rdim: dimensionality of each range (0 dim of index PDL)
  730 #
  731 # ntsize: number of nonzero size dimensions
  732 # sizes:  array of range sizes, indexed (0..rdim-1).  A zero element means
  733 #         that the dimension is omitted from the child dim list.
  734 # corners: parent coordinates of each corner, running fastest over coord index.
  735 #       (indexed 0 .. (nitems-1)*(rdim)+rdim-1)
  736 # nitems: total number of list elements   (product of itdims)
  737 # itdim:  number of index thread dimensions
  738 # itdims: Size of each index thread dimension, indexed (0..itdim-1)
  739 #
  740 # bsize: Number of independently specified boundary conditions
  741 # nsizes: Number of independently specified range dim sizes
  742 # boundary: Array containing all the boundary condition specs
  743 # indord: Order/size of the indexing dim (0th dim of $index)
  744 
  745         Comp => 'PDL_Indx rdim;
  746                  PDL_Indx nitems;
  747                  PDL_Indx itdim;
  748                  PDL_Indx ntsize;
  749                  PDL_Indx bsize;
  750                  PDL_Indx nsizes;
  751                  PDL_Indx sizes[$COMP(rdim)];
  752                  PDL_Indx itdims[$COMP(itdim)];
  753                  PDL_Indx corners[$COMP(rdim) * $COMP(nitems)];
  754                  char boundary[$COMP(rdim)];
  755                  ',
  756         MakeComp => <<'EOD-MakeComp',
  757 pdl *ind_pdl;
  758 pdl *size_pdl;
  759 
  760 /***
  761  * Check and condition the index piddle.  Some of this is apparently
  762  * done by XS -- but XS doesn't check for existing SVs that are undef.
  763  */
  764 if ((index==NULL) || (index == &PL_sv_undef))
  765    { $CROAK("rangeb: index variable must be defined"); }
  766 
  767 if(!(ind_pdl = PDL->SvPDLV(index))) /* assignment */
  768    { $CROAK("rangeb: unable to convert index variable to a PDL"); }
  769 
  770 PDL->make_physdims(ind_pdl);
  771 
  772 /* Generalized empties are ok, but not in the special 0 dim (the index vector) */
  773 if(ind_pdl->dims[0] == 0)
  774     { $CROAK("rangeb: can't handle Empty indices -- call range instead"); }
  775 
  776 /***
  777  * Ensure that the index is a PDL_Indx.  If there's no loss of information,
  778  * just upgrade it -- otherwise, make a temporary copy.
  779  */
  780 switch(ind_pdl->datatype) {
  781  default:                              /* Most types: */
  782    ind_pdl = PDL->hard_copy(ind_pdl);  /*   copy and fall through */
  783  case PDL_B: case PDL_S: case PDL_US: case PDL_L: case PDL_LL:
  784    PDL->converttype(&ind_pdl,PDL_IND,1); /* convert in place. */
  785    break;
  786  case PDL_IND:
  787    PDL->make_physical(ind_pdl);
  788    break;
  789 }
  790 
  791 /***
  792  * Figure sizes of the COMP arrrays and allocate them.
  793  */
  794 {
  795   PDL_Indx i,nitems;
  796 
  797   $COMP(rdim) = ind_pdl->ndims ? ind_pdl->dims[0] : 1;
  798   for(i=nitems=1; i < ind_pdl->ndims; i++)  /* Accumulate item list size */
  799     nitems *= ind_pdl->dims[i];
  800   $COMP(nitems) = nitems;
  801   $COMP(itdim) = ind_pdl->ndims ? ind_pdl->ndims - 1 : 0;
  802   $DOCOMPDIMS();
  803 }
  804 
  805 /***
  806  * Fill in the boundary condition array
  807  */
  808 {
  809   char *bstr;
  810   STRLEN blen;
  811   bstr = SvPV(boundary,blen);
  812 
  813   if(blen == 0) {
  814     /* If no boundary is specified then every dim gets forbidden */
  815     int i;
  816     for (i=0;i<$COMP(rdim);i++)
  817       $COMP(boundary[i]) = 0;
  818   } else {
  819     int i;
  820     for(i=0;i<$COMP(rdim);i++) {
  821       switch(bstr[i < blen ? i : blen-1 ]) {
  822       case '0': case 'f': case 'F':               /* forbid */
  823         $COMP(boundary[i]) = 0;
  824         break;
  825       case '1': case 't': case 'T':               /* truncate */
  826         $COMP(boundary[i]) = 1;
  827         break;
  828       case '2': case 'e': case 'E':               /* extend */
  829         $COMP(boundary[i]) = 2;
  830         break;
  831       case '3': case 'p': case 'P':               /* periodic */
  832         $COMP(boundary[i]) = 3;
  833         break;
  834       case '4': case 'm': case 'M':               /* mirror */
  835         $COMP(boundary[i]) = 4;
  836         break;
  837       default:
  838         {
  839           /* No need to check if i < blen -- this will barf out the
  840            * first time it gets hit.  I didn't use $ CROAK 'coz that
  841            * macro doesn't let you pass in a string variable -- only a
  842            * constant.
  843            */
  844           barf("Error in rangeb: Unknown boundary condition '%c' in range",bstr[i]);
  845         }
  846         break;
  847       } // end of switch
  848     }
  849   }
  850 }
  851 /***
  852  * Store the sizes of the index-thread dims
  853  */
  854 {
  855   PDL_Indx i;
  856   PDL_Indx nd = ind_pdl->ndims - 1;
  857   for(i=0; i < nd ; i++)
  858     $COMP(itdims[i]) = ind_pdl->dims[i+1];
  859 }
  860 
  861 /***
  862  * Check and condition the size piddle, and store sizes of the ranges
  863  */
  864 {
  865   PDL_Indx i,ntsize;
  866 
  867   if( (size == NULL) || (size == &PL_sv_undef) ) {
  868     // NO size was passed in (not normally executed even if you passed in no size to range(),
  869     // as range() generates a size array...
  870     for(i=0;i<$COMP(rdim);i++)
  871           $COMP(sizes[i]) = 0;
  872 
  873   } else {
  874     /* Normal case with sizes present in a PDL */
  875 
  876     if(!(size_pdl = PDL->SvPDLV(size))) /* assignment */
  877       $CROAK("Unable to convert size to a PDL in range");
  878 
  879     if(size_pdl->nvals == 0) {
  880       // no values in the size_pdl - Empty or Null.  Just copy 0s to all the range dims
  881       for(i=0;i<$COMP(rdim);i++)
  882         $COMP(sizes[i]) = 0;
  883 
  884     } else {
  885 
  886       // Convert size PDL to PDL_IND to support indices
  887       switch(size_pdl->datatype) {
  888       default:                              /* Most types: */
  889         size_pdl = PDL->hard_copy(size_pdl);  /*   copy and fall through */
  890       case PDL_B: case PDL_S: case PDL_US: case PDL_L:  case PDL_LL:
  891         PDL->converttype(&size_pdl,PDL_IND,1); /* convert in place. */
  892         break;
  893       case PDL_IND:
  894         PDL->make_physical(size_pdl);
  895         break;
  896       }
  897 
  898       $COMP(nsizes) = size_pdl->nvals; /* Store for later permissiveness check */
  899 
  900       /* Copy the sizes, or die if they're the wrong shape */
  901       if(size_pdl->nvals == 1) {
  902         for(i=0;i<$COMP(rdim);i++) {
  903           $COMP(sizes[i]) = *((PDL_Indx *)(size_pdl->data));
  904         }
  905 
  906         /* Check for nonnegativity of sizes.  The rdim>0 mask ensures that */
  907         /* we don't barf on the Empty PDL (as an index). */
  908         if( $COMP(rdim) > 0 && $COMP(sizes[0]) < 0 ) {
  909           $CROAK("  Negative range size is not allowed in range\n");
  910         }
  911       }
  912       else if( size_pdl->nvals <= $COMP(rdim) && size_pdl->ndims == 1) {
  913         for(i=0;i<$COMP(rdim);i++) {
  914           $COMP(sizes[i]) = (   (i < size_pdl->nvals) ?
  915                                 ((PDL_Indx *)(size_pdl->data))[i] :
  916                                 0
  917                             );
  918           if($COMP(sizes[i]) < 0)
  919                 $CROAK("  Negative range sizes are not allowed in range\n");
  920         }
  921       }
  922       else {
  923         $CROAK(" Size must match index's 0th dim in range\n");
  924       }
  925 
  926     } /* end of nonempty size-piddle code */
  927   } /* end of defined-size-piddle code */
  928 
  929   /* Insert the number of nontrivial sizes (these get output dimensions) */
  930   for(i=ntsize=0;i<$COMP(rdim);i++)
  931     if($COMP(sizes[i]))
  932       ntsize++;
  933   $COMP(ntsize) = ntsize;
  934 }
  935 
  936 /***
  937  * Stash coordinates of the corners
  938  */
  939 
  940 {
  941   PDL_Indx i,j,k,ioff;
  942   PDL_Indx *cptr;
  943   PDL_Indx *iter = (PDL_Indx *)(PDL->smalloc((STRLEN) (sizeof(PDL_Indx) * ($COMP(itdim)))));
  944 
  945   /* initialize iterator to loop over index threads */
  946   cptr = iter;
  947   for(k=0;k<$COMP(itdim);k++)
  948     *(cptr++) = 0;
  949 
  950   cptr = $COMP(corners);
  951   do {
  952     /* accumulate offset into the index from the iterator */
  953     for(k=ioff=0;k<$COMP(itdim);k++)
  954       ioff += iter[k] * ind_pdl->dimincs[k+1];
  955 
  956     /* Loop over the 0th dim of index, copying coords. */
  957     /* This is the natural place to check for permissive ranging; too */
  958     /* bad we don't have access to the parent piddle here... */
  959 
  960     for(j=0;j<$COMP(rdim);j++)
  961         *(cptr++) = ((PDL_Indx *)(ind_pdl->data))[ioff + ind_pdl->dimincs[0] * j];
  962 
  963     /* Increment the iterator -- the test increments, the body carries. */
  964     for(k=0; k<$COMP(itdim) && (++(iter[k]))>=($COMP(itdims)[k]) ;k++)
  965       iter[k] = 0;
  966   } while(k<$COMP(itdim));
  967 
  968 
  969 
  970 }
  971 
  972 
  973 $SETREVERSIBLE(1);
  974 
  975 EOD-MakeComp
  976 
  977         RedoDims => <<'EOD-RedoDims' ,
  978 {
  979   PDL_Indx stdim = $PARENT(ndims) - $COMP(rdim);
  980   PDL_Indx dim,inc;
  981   PDL_Indx i,rdvalid;
  982 
  983     // Speed bump for ludicrous cases
  984     if( $COMP(rdim) > $PARENT(ndims)+5 && $COMP(nsizes) != $COMP(rdim)) {
  985       barf("Ludicrous number of extra dims in range index; leaving child null.\n    (%d implicit dims is > 5; index has %d dims; source has %d dim%s.)\n    This often means that your index PDL is incorrect.  To avoid this message,\n    allocate dummy dims in the source or use %d dims in range's size field.\n",$COMP(rdim)-$PARENT(ndims),$COMP(rdim),$PARENT(ndims),($PARENT(ndims))>1?"s":"",$COMP(rdim));
  986     }
  987 
  988     if(stdim < 0)
  989       stdim = 0;
  990 
  991     /* Set dimensionality of child */
  992     $CHILD(ndims) = $COMP(itdim) + $COMP(ntsize) + stdim;
  993     $SETNDIMS($COMP(itdim)+$COMP(ntsize)+stdim);
  994 
  995     inc = 1;
  996     /* Copy size dimensions to child, crunching as we go. */
  997     dim = $COMP(itdim);
  998     for(i=rdvalid=0;i<$COMP(rdim);i++) {
  999       if($COMP(sizes[i])) {
 1000         rdvalid++;
 1001         $CHILD(dimincs[dim]) = inc;
 1002         inc *= ($CHILD(dims[dim++]) = $COMP(sizes[i])); /* assignment */
 1003       }
 1004     }
 1005 
 1006     /* Copy index thread dimensions to child */
 1007     for(dim=0; dim<$COMP(itdim); dim++) {
 1008       $CHILD(dimincs[dim]) = inc;
 1009       inc *= ($CHILD(dims[dim]) = $COMP(itdims[dim])); /* assignment */
 1010     }
 1011 
 1012     /* Copy source thread dimensions to child */
 1013     dim = $COMP(itdim) + rdvalid;
 1014     for(i=0;i<stdim;i++) {
 1015       $CHILD(dimincs[dim]) = inc;
 1016       inc *= ($CHILD(dims[dim++]) = $PARENT(dims[i+$COMP(rdim)])); /* assignment */
 1017     }
 1018 
 1019     /* Cover bizarre case where the source PDL is empty - in that case, change */
 1020     /* all non-barfing boundary conditions to truncation, since we have no data */
 1021     /* to reflect, extend, or mirror. */
 1022     if($PARENT(dims[0])==0) {
 1023       for(dim=0; dim<$COMP(rdim); dim++) {
 1024         if($COMP(boundary[dim]))
 1025           $COMP(boundary[dim]) = 1; // force truncation
 1026       }
 1027     }
 1028 
 1029 
 1030   $CHILD(datatype) = $PARENT(datatype);
 1031 
 1032   $SETDIMS();
 1033 }
 1034 
 1035 EOD-RedoDims
 1036 
 1037         EquivCPOffsCode => <<'EOD-EquivCPOffsCode',
 1038 {
 1039   PDL_Indx *iter, *ip;  /* vector iterator */
 1040   PDL_Indx *sizes, *sp; /* size vector including stdims */
 1041   PDL_Indx *coords;     /* current coordinates */
 1042 
 1043   PDL_Indx k;           /* index */
 1044   PDL_Indx item;        /* index thread iterator */
 1045   PDL_Indx pdim = $PARENT_P(ndims);
 1046   PDL_Indx rdim = $COMP(rdim);
 1047   PDL_Indx prdim = (rdim < pdim) ? rdim : pdim;
 1048   PDL_Indx stdim = pdim - prdim;
 1049 
 1050   /* Allocate iterator and larger size vector -- do it all in one foop
 1051    * to avoid extra calls to smalloc.
 1052    */
 1053     if(!(iter = (PDL_Indx *)(PDL->smalloc((STRLEN) (sizeof(PDL_Indx) * ($PARENT_P(ndims) * 2 + rdim)))))) {
 1054     barf("couldn't get memory for range iterator");
 1055   }
 1056   sizes  = iter + $PARENT_P(ndims);
 1057   coords = sizes + $PARENT_P(ndims);
 1058 
 1059   /* Figure out size vector */
 1060   for(ip = $COMP(sizes), sp = sizes, k=0; k<rdim; k++)
 1061      *(sp++) = *(ip++);
 1062   for(; k < $PARENT_P(ndims); k++)
 1063      *(sp++) = $PARENT_P(dims[k]);
 1064 
 1065 
 1066   /* Loop over all the ranges in the index list */
 1067   for(item=0; item<$COMP(nitems); item++) {
 1068 
 1069     /* initialize in-range iterator to loop within each range */
 1070     for(ip = iter, k=0; k<$PARENT_P(ndims); k++)
 1071       *(ip++) = 0;
 1072 
 1073     do {
 1074       PDL_Indx poff = 0;
 1075       PDL_Indx coff;
 1076       PDL_Indx k2;
 1077       char trunc = 0;       /* Flag used to skip truncation case */
 1078 
 1079       /* Collect and boundary-check the current N-D coords */
 1080       for(k=0; k < prdim; k++){
 1081 
 1082         PDL_Indx ck = iter[k] + $COMP(corners[ item * rdim + k  ]) ;
 1083 
 1084         /* normal case */
 1085           if(ck < 0 || ck >= $PARENT_P(dims[k])) {
 1086             switch($COMP(boundary[k])) {
 1087             case 0: /* no boundary breakage allowed */
 1088           {
 1089         char barfstr[1024];
 1090         sprintf(barfstr,"index out-of-bounds in range (index vector #%ld)",item);
 1091         barf(barfstr);
 1092           }
 1093               break;
 1094             case 1: /* truncation */
 1095               trunc = 1;
 1096               break;
 1097             case 2: /* extension -- crop */
 1098               ck = (ck >= $PARENT_P(dims[k])) ? $PARENT_P(dims[k])-1 : 0;
 1099               break;
 1100             case 3: /* periodic -- mod it */
 1101               ck %= $PARENT_P(dims[k]);
 1102               if(ck < 0)   /* Fix mod breakage in C */
 1103                 ck += $PARENT_P(dims[k]);
 1104               break;
 1105             case 4: /* mirror -- reflect off the edges */
 1106               ck += $PARENT_P(dims[k]);
 1107               ck %= ($PARENT_P(dims[k]) * 2);
 1108               if(ck < 0) /* Fix mod breakage in C */
 1109                 ck += $PARENT_P(dims[k])*2;
 1110               ck -= $PARENT_P(dims[k]);
 1111               if(ck < 0) {
 1112                 ck *= -1;
 1113                 ck -= 1;
 1114               }
 1115               break;
 1116             default:
 1117               barf("Unknown boundary condition in range -- bug alert!");
 1118               break;
 1119             }
 1120           }
 1121 
 1122         coords[k] = ck;
 1123 
 1124       }
 1125 
 1126       /* Check extra dimensions -- pick up where k left off... */
 1127       for( ; k < rdim ; k++) {
 1128         /* Check for indexing off the end of the dimension list */
 1129 
 1130         PDL_Indx ck = iter[k] + $COMP(corners[ item * rdim + k  ]) ;
 1131 
 1132         switch($COMP(boundary[k])) {
 1133             case 0: /* No boundary breakage allowed -- nonzero corners cause barfage */
 1134               if(ck != 0)
 1135                  barf("Too many dims in range index (and you've forbidden boundary violations)");
 1136               break;
 1137             case 1: /* truncation - just truncate if the corner is nonzero */
 1138               trunc |= (ck != 0);
 1139               break;
 1140             case 2: /* extension -- ignore the corner (same as 3) */
 1141             case 3: /* periodic  -- ignore the corner */
 1142             case 4: /* mirror -- ignore the corner */
 1143               ck = 0;
 1144               break;
 1145             default:
 1146               barf("Unknown boundary condition in range -- bug alert!");
 1147               break;
 1148           }
 1149       }
 1150 
 1151       /* Find offsets into the child and parent arrays, from the N-D coords */
 1152       /* Note we only loop over real source dims (prdim) to accumulate -- */
 1153       /* because the offset is trivial and/or we're truncating for virtual */
 1154       /* dims caused by permissive ranging. */
 1155       coff = $CHILD_P(dimincs[0]) * item;
 1156       for(k2 = $COMP(itdim), poff = k = 0;
 1157           k < prdim;
 1158           k++) {
 1159         poff += coords[k]*$PARENT_P(dimincs[k]);
 1160         if($COMP(sizes[k]))
 1161           coff += iter[k] * $CHILD_P(dimincs[k2++]);
 1162       }
 1163 
 1164       /* Loop the copy over all the source thread dims (above rdim). */
 1165       do {
 1166         PDL_Indx poff1 = poff;
 1167         PDL_Indx coff1 = coff;
 1168 
 1169         /* Accumulate the offset due to source threading */
 1170         for(k2 = $COMP(itdim) + $COMP(ntsize), k = rdim;
 1171             k < pdim;
 1172             k++) {
 1173           poff1 += iter[k] * $PARENT_P(dimincs[k]);
 1174           coff1 += iter[k] * $CHILD_P(dimincs[k2++]);
 1175         }
 1176 
 1177         /* Finally -- make the copy
 1178          * EQUIVCPTRUNC works like EQUIVCPOFFS but with checking for
 1179          * out-of-bounds conditions.
 1180          */
 1181         $EQUIVCPTRUNC(coff1,poff1,trunc);
 1182 
 1183         /* Increment the source thread iterator */
 1184         for( k=$COMP(rdim);
 1185              k < $PARENT_P(ndims) && (++(iter[k]) >= $PARENT_P(dims[k]));
 1186              k++)
 1187           iter[k] = 0;
 1188       } while(k < $PARENT_P(ndims)); /* end of source-thread iteration */
 1189 
 1190       /* Increment the in-range iterator */
 1191       for(k = 0;
 1192           k < $COMP(rdim) && (++(iter[k]) >= $COMP(sizes[k]));
 1193           k++)
 1194         iter[k] = 0;
 1195     } while(k < $COMP(rdim)); /* end of main iteration */
 1196   } /* end of item do loop */
 1197 
 1198 }
 1199 
 1200 EOD-EquivCPOffsCode
 1201 
 1202 );
 1203 
 1204 
 1205 =head2 rld
 1206 =cut
 1207 
 1208 pp_def(
 1209         'rld',
 1210         Pars=>'indx a(n); b(n); [o]c(m);',
 1211         PMCode =><<'EOD',
 1212 sub PDL::rld {
 1213   my ($x,$y) = @_;
 1214   my ($c);
 1215   if ($#_ == 2) {
 1216     $c = $_[2];
 1217   } else {
 1218 # XXX Need to improve emulation of threading in auto-generating c
 1219     my ($size) = $x->sumover->max;
 1220     my (@dims) = $x->dims;
 1221     shift @dims;
 1222     $c = $y->zeroes($size,@dims);
 1223   }
 1224   &PDL::_rld_int($x,$y,$c);
 1225   $c;
 1226 }
 1227 EOD
 1228         Code=>'
 1229           PDL_Indx i,j=0,an;
 1230           $GENERIC(b) bv;
 1231           loop (n) %{
 1232             an = $a();
 1233             bv = $b();
 1234             for (i=0;i<an;i++) {
 1235               $c(m=>j) = bv;
 1236               j++;
 1237             }
 1238           %}',
 1239         Doc => <<'EOD'
 1240 =for ref
 1241 
 1242 Run-length decode a vector
 1243 
 1244 Given a vector C<$x> of the numbers of instances of values C<$y>, run-length
 1245 decode to C<$c>.
 1246 
 1247 =for example
 1248 
 1249  rld($x,$y,$c=null);
 1250 
 1251 =cut
 1252 
 1253 EOD
 1254 );
 1255 
 1256 =head2 rle
 1257 =cut
 1258 
 1259 pp_def(
 1260         'rle',
 1261         Pars=>'c(n); indx [o]a(m); [o]b(m);',
 1262 #this RedoDimsCode sets $SIZE(m)==$SIZE(n), but the slice in the PMCode below makes m<=n.
 1263         RedoDimsCode=>'$SIZE(m)=$PDL(c)->dims[0];',
 1264         PMCode=><<'EOC',
 1265 sub PDL::rle {
 1266   my $c = shift;
 1267   my ($x,$y) = @_==2 ? @_ : (null,null);
 1268   &PDL::_rle_int($c,$x,$y);
 1269   my $max_ind = ($c->ndims<2) ? ($x!=0)->sumover-1 :
 1270                                 ($x!=0)->clump(1..$x->ndims-1)->sumover->max-1;
 1271   return ($x->slice("0:$max_ind"),$y->slice("0:$max_ind"));
 1272 }
 1273 EOC
 1274         Code=>'
 1275           PDL_Indx j=0,sn=$SIZE(n);
 1276           $GENERIC(c) cv, clv;
 1277           clv = $c(n=>0);
 1278           $b(m=>0) = clv;
 1279           $a(m=>0) = 0;
 1280           loop (n) %{
 1281             cv = $c();
 1282             if (cv == clv) {
 1283               $a(m=>j)++;
 1284             } else {
 1285               j++;
 1286               $b(m=>j) = clv = cv;
 1287               $a(m=>j) = 1;
 1288             }
 1289           %}
 1290           for (j++;j<$SIZE(m);j++) {
 1291             $a(m=>j) = 0;
 1292             $b(m=>j) = 0;
 1293           }
 1294         ',
 1295         Doc => <<'EOD'
 1296 =for ref
 1297 
 1298 Run-length encode a vector
 1299 
 1300 Given vector C<$c>, generate a vector C<$x> with the number of each
 1301 element, and a vector C<$y> of the unique values.  New in PDL 2.017,
 1302 only the elements up to the first instance of C<0> in C<$x> are
 1303 returned, which makes the common use case of a 1-dimensional C<$c> simpler.
 1304 For threaded operation, C<$x> and C<$y> will be large enough
 1305 to hold the largest row of C<$y>, and only the elements up to the
 1306 first instance of C<0> in each row of C<$x> should be considered.
 1307 
 1308 =for example
 1309 
 1310  $c = floor(4*random(10));
 1311  rle($c,$x=null,$y=null);
 1312  #or
 1313  ($x,$y) = rle($c);
 1314 
 1315  #for $c of shape [10, 4]:
 1316  $c = floor(4*random(10,4));
 1317  ($x,$y) = rle($c);
 1318 
 1319  #to see the results of each row one at a time:
 1320  foreach (0..$c->dim(1)-1){
 1321   my ($as,$bs) = ($x(:,($_)),$y(:,($_)));
 1322   my ($ta,$tb) = where($as,$bs,$as!=0); #only the non-zero elements of $x
 1323   print $c(:,($_)) . " rle==> " , ($ta,$tb) , "\trld==> " . rld($ta,$tb) . "\n";
 1324  }
 1325 
 1326 =cut
 1327 
 1328 EOD
 1329 );
 1330 
 1331 # this one can convert vaffine piddles without(!) physicalising them
 1332 # maybe it can replace 'converttypei' in the future?
 1333 #
 1334 # XXX do not know whether the HandleBad stuff will work here
 1335 #
 1336 pp_def('flowconvert',
 1337        HandleBad => 1,
 1338        DefaultFlow => 1,
 1339        Reversible => 1,
 1340        Pars => 'PARENT(); [oca]CHILD()',
 1341        OtherPars => 'int totype;',
 1342        Reversible => 1,
 1343        # Forced types
 1344        FTypes => {CHILD => '$COMP(totype)'},
 1345        Code =>
 1346        '$CHILD() = $PARENT();',
 1347        BadCode =>
 1348        'if ( $ISBAD(PARENT()) ) {
 1349            $SETBAD(CHILD());
 1350         } else {
 1351            $CHILD() = $PARENT();
 1352         }',
 1353        BackCode => '$PARENT() = $CHILD();',
 1354        BadBackCode =>
 1355        'if ( $ISBAD(CHILD()) ) {
 1356            $SETBAD(PARENT());
 1357         } else {
 1358            $PARENT() = $CHILD();
 1359         }',
 1360        Doc => 'internal',
 1361 );
 1362 
 1363 
 1364 pp_def(
 1365         'converttypei',
 1366         HandleBad => 1,
 1367         DefaultFlow => 1,
 1368         GlobalNew => 'converttypei_new',
 1369         OtherPars => 'int totype;',
 1370         P2Child => 1,
 1371         NoPdlThread => 1,
 1372         Identity => 1,
 1373         Reversible => 1,
 1374 # Forced types
 1375         FTypes => {CHILD => '$COMP(totype)'},
 1376         Doc => 'internal',
 1377 );
 1378 
 1379 
 1380 
 1381 # the perl wrapper clump is now defined in Core.pm
 1382 # this is just the low level interface
 1383 pp_def(
 1384         '_clump_int',
 1385         DefaultFlow => 1,
 1386         OtherPars => 'int n',
 1387         P2Child => 1,
 1388     NoPdlThread=>1,
 1389         Priv => 'int nnew; int nrem;',
 1390         RedoDims => 'int i; PDL_Indx d1;
 1391 
 1392              /* truncate overly long clumps to just clump existing dimensions */
 1393          if($COMP(n) > $PARENT(ndims))
 1394                         $COMP(n) = $PARENT(ndims);
 1395 
 1396          if($COMP(n) < -1)
 1397                         $COMP(n) = $PARENT(ndims) + $COMP(n) + 1;
 1398 
 1399                  $PRIV(nrem) = ($COMP(n)==-1 ? $PARENT(threadids[0]) : $COMP(n));
 1400                  $PRIV(nnew) = $PARENT(ndims) - $PRIV(nrem) + 1;
 1401                  $SETNDIMS($PRIV(nnew));
 1402                  d1=1;
 1403                  for(i=0; i<$PRIV(nrem); i++) {
 1404                         d1 *= $PARENT(dims[i]);
 1405                  }
 1406                  $CHILD(dims[0]) = d1;
 1407                  for(; i<$PARENT(ndims); i++) {
 1408                         $CHILD(dims[i-$PRIV(nrem)+1]) = $PARENT(dims[i]);
 1409                  }
 1410                  $SETDIMS();
 1411                  $SETDELTATHREADIDS(1-$PRIV(nrem));
 1412                  ',
 1413         EquivCPOffsCode => '
 1414                 PDL_Indx i;
 1415                 for(i=0; i<$CHILD_P(nvals); i++) {
 1416                         $EQUIVCPOFFS(i,i);
 1417                 }
 1418                 ',
 1419         Reversible => 1,
 1420         Doc => 'internal',
 1421 );
 1422 
 1423 
 1424 
 1425 =head2 xchg
 1426 =cut
 1427 
 1428 pp_def(
 1429         'xchg',
 1430         OtherPars => 'int n1; int n2;',
 1431         DefaultFlow => 1,
 1432         Reversible => 1,
 1433         P2Child => 1,
 1434         NoPdlThread => 1,
 1435         XCHGOnly => 1,
 1436         EquivDimCheck => 'if ($COMP(n1) <0)
 1437                                 $COMP(n1) += $PARENT(threadids[0]);
 1438                           if ($COMP(n2) <0)
 1439                                 $COMP(n2) += $PARENT(threadids[0]);
 1440                           if ($COMP(n1) <0 ||$COMP(n2) <0 ||
 1441                              $COMP(n1) >= $PARENT(threadids[0]) ||
 1442                              $COMP(n2) >= $PARENT(threadids[0]))
 1443                 barf("One of dims %d, %d out of range: should be 0<=dim<%d",
 1444                         $COMP(n1),$COMP(n2),$PARENT(threadids[0]));',
 1445         EquivPDimExpr => '(($CDIM == $COMP(n1)) ? $COMP(n2) : ($CDIM == $COMP(n2)) ? $COMP(n1) : $CDIM)',
 1446         EquivCDimExpr => '(($PDIM == $COMP(n1)) ? $COMP(n2) : ($PDIM == $COMP(n2)) ? $COMP(n1) : $PDIM)',
 1447         Doc => <<'EOD',
 1448 =for ref
 1449 
 1450 exchange two dimensions
 1451 
 1452 Negative dimension indices count from the end.
 1453 
 1454 The command
 1455 
 1456 =for example
 1457 
 1458  $y = $x->xchg(2,3);
 1459 
 1460 creates C<$y> to be like C<$x> except that the dimensions 2 and 3
 1461 are exchanged with each other i.e.
 1462 
 1463  $y->at(5,3,2,8) == $x->at(5,3,8,2)
 1464 
 1465 =cut
 1466 
 1467 EOD
 1468 );
 1469 
 1470 pp_addpm(<< 'EOD');
 1471 
 1472 =head2 reorder
 1473 
 1474 =for ref
 1475 
 1476 Re-orders the dimensions of a PDL based on the supplied list.
 1477 
 1478 Similar to the L<xchg|/xchg> method, this method re-orders the dimensions
 1479 of a PDL. While the L<xchg|/xchg> method swaps the position of two dimensions,
 1480 the reorder method can change the positions of many dimensions at
 1481 once.
 1482 
 1483 =for usage
 1484 
 1485  # Completely reverse the dimension order of a 6-Dim array.
 1486  $reOrderedPDL = $pdl->reorder(5,4,3,2,1,0);
 1487 
 1488 The argument to reorder is an array representing where the current dimensions
 1489 should go in the new array. In the above usage, the argument to reorder
 1490 C<(5,4,3,2,1,0)>
 1491 indicates that the old dimensions (C<$pdl>'s dims) should be re-arranged to make the
 1492 new pdl (C<$reOrderPDL>) according to the following:
 1493 
 1494    Old Position   New Position
 1495    ------------   ------------
 1496    5              0
 1497    4              1
 1498    3              2
 1499    2              3
 1500    1              4
 1501    0              5
 1502 
 1503 You do not need to specify all dimensions, only a complete set
 1504 starting at position 0.  (Extra dimensions are left where they are).
 1505 This means, for example, that you can reorder() the X and Y dimensions of
 1506 an image, and not care whether it is an RGB image with a third dimension running
 1507 across color plane.
 1508 
 1509 =for example
 1510 
 1511 Example:
 1512 
 1513  pdl> $x = sequence(5,3,2);       # Create a 3-d Array
 1514  pdl> p $x
 1515  [
 1516   [
 1517    [ 0  1  2  3  4]
 1518    [ 5  6  7  8  9]
 1519    [10 11 12 13 14]
 1520   ]
 1521   [
 1522    [15 16 17 18 19]
 1523    [20 21 22 23 24]
 1524    [25 26 27 28 29]
 1525   ]
 1526  ]
 1527  pdl> p $x->reorder(2,1,0); # Reverse the order of the 3-D PDL
 1528  [
 1529   [
 1530    [ 0 15]
 1531    [ 5 20]
 1532    [10 25]
 1533   ]
 1534   [
 1535    [ 1 16]
 1536    [ 6 21]
 1537    [11 26]
 1538   ]
 1539   [
 1540    [ 2 17]
 1541    [ 7 22]
 1542    [12 27]
 1543   ]
 1544   [
 1545    [ 3 18]
 1546    [ 8 23]
 1547    [13 28]
 1548   ]
 1549   [
 1550    [ 4 19]
 1551    [ 9 24]
 1552    [14 29]
 1553   ]
 1554  ]
 1555 
 1556 The above is a simple example that could be duplicated by calling
 1557 C<$x-E<gt>xchg(0,2)>, but it demonstrates the basic functionality of reorder.
 1558 
 1559 As this is an index function, any modifications to the
 1560 result PDL will change the parent.
 1561 
 1562 =cut
 1563 
 1564 sub PDL::reorder {
 1565         my ($pdl,@newDimOrder) = @_;
 1566 
 1567         my $arrayMax = $#newDimOrder;
 1568 
 1569         #Error Checking:
 1570         if( $pdl->getndims < scalar(@newDimOrder) ){
 1571                 my $errString = "PDL::reorder: Number of elements (".scalar(@newDimOrder).") in newDimOrder array exceeds\n";
 1572                 $errString .= "the number of dims in the supplied PDL (".$pdl->getndims.")";
 1573                 barf($errString);
 1574         }
 1575 
 1576         # Check to make sure all the dims are within bounds
 1577         for my $i(0..$#newDimOrder) {
 1578           my $dim = $newDimOrder[$i];
 1579           if($dim < 0 || $dim > $#newDimOrder) {
 1580               my $errString = "PDL::reorder: Dim index $newDimOrder[$i] out of range in position $i\n(range is 0-$#newDimOrder)";
 1581               barf($errString);
 1582           }
 1583         }
 1584 
 1585         # Checking that they are all present and also not duplicated is done by thread() [I think]
 1586 
 1587         # a quicker way to do the reorder
 1588         return $pdl->thread(@newDimOrder)->unthread(0);
 1589 }
 1590 
 1591 EOD
 1592 
 1593 =head2 mv
 1594 =cut
 1595 
 1596 pp_def(
 1597         'mv',
 1598         OtherPars => 'int n1; int n2;',
 1599         DefaultFlow => 1,
 1600         Reversible => 1,
 1601         P2Child => 1,
 1602         NoPdlThread => 1,
 1603         XCHGOnly => 1,
 1604         EquivDimCheck => 'if ($COMP(n1) <0)
 1605                                 $COMP(n1) += $PARENT(threadids[0]);
 1606                           if ($COMP(n2) <0)
 1607                                 $COMP(n2) += $PARENT(threadids[0]);
 1608                           if ($COMP(n1) <0 ||$COMP(n2) <0 ||
 1609                              $COMP(n1) >= $PARENT(threadids[0]) ||
 1610                              $COMP(n2) >= $PARENT(threadids[0]))
 1611                 barf("One of dims %d, %d out of range: should be 0<=dim<%d",
 1612                         $COMP(n1),$COMP(n2),$PARENT(threadids[0]));',
 1613         EquivPDimExpr => '(($COMP(n1) < $COMP(n2)) ?
 1614         (($CDIM < $COMP(n1) || $CDIM > $COMP(n2)) ?
 1615                 $CDIM : (($CDIM == $COMP(n2)) ? $COMP(n1) : $CDIM+1))
 1616         : (($COMP(n2) < $COMP(n1)) ?
 1617                 (($CDIM > $COMP(n1) || $CDIM < $COMP(n2)) ?
 1618                         $CDIM : (($CDIM == $COMP(n2)) ? $COMP(n1) : $CDIM-1))
 1619                 : $CDIM))',
 1620         EquivCDimExpr => '(($COMP(n2) < $COMP(n1)) ?
 1621         (($PDIM < $COMP(n2) || $PDIM > $COMP(n1)) ?
 1622                 $PDIM : (($PDIM == $COMP(n1)) ? $COMP(n2) : $PDIM+1))
 1623         : (($COMP(n1) < $COMP(n2)) ?
 1624                 (($PDIM > $COMP(n2) || $PDIM < $COMP(n1)) ?
 1625                         $PDIM : (($PDIM == $COMP(n1)) ? $COMP(n2) : $PDIM-1))
 1626                 : $PDIM))',
 1627         Doc => << 'EOD',
 1628 =for ref
 1629 
 1630 move a dimension to another position
 1631 
 1632 The command
 1633 
 1634 =for example
 1635 
 1636  $y = $x->mv(4,1);
 1637 
 1638 creates C<$y> to be like C<$x> except that the dimension 4 is moved to the
 1639 place 1, so:
 1640 
 1641  $y->at(1,2,3,4,5,6) == $x->at(1,5,2,3,4,6);
 1642 
 1643 The other dimensions are moved accordingly.
 1644 Negative dimension indices count from the end.
 1645 =cut
 1646 EOD
 1647 );
 1648 
 1649 pp_addhdr << 'EOH';
 1650 #define sign(x) ( (x) < 0 ? -1 : 1)
 1651 EOH
 1652 
 1653 =head2 oslice
 1654 =cut
 1655 
 1656 # I think the quotes in the =item ":" lines
 1657 # confuse the perldoc stuff
 1658 #
 1659 pp_def(
 1660         'oslice',
 1661         Doc => << 'EOD',
 1662 =for ref
 1663 
 1664 DEPRECATED:  'oslice' is the original 'slice' routine in pre-2.006_006
 1665 versions of PDL.  It is left here for reference but will disappear in
 1666 PDL 3.000
 1667 
 1668 Extract a rectangular slice of a piddle, from a string specifier.
 1669 
 1670 C<slice> was the original Swiss-army-knife PDL indexing routine, but is
 1671 largely superseded by the L<NiceSlice|PDL::NiceSlice> source prefilter
 1672 and its associated L<nslice|PDL::Core/nslice> method.  It is still used as the
 1673 basic underlying slicing engine for L<nslice|PDL::Core/nslice>,
 1674 and is especially useful in particular niche applications.
 1675 
 1676 =for example
 1677 
 1678  $x->slice('1:3');  #  return the second to fourth elements of $x
 1679  $x->slice('3:1');  #  reverse the above
 1680  $x->slice('-2:1'); #  return last-but-one to second elements of $x
 1681 
 1682 The argument string is a comma-separated list of what to do
 1683 for each dimension. The current formats include
 1684 the following, where I<a>, I<b> and I<c> are integers and can
 1685 take legal array index values (including -1 etc):
 1686 
 1687 =over 8
 1688 
 1689 =item :
 1690 
 1691 takes the whole dimension intact.
 1692 
 1693 =item ''
 1694 
 1695 (nothing) is a synonym for ":"
 1696 (This means that C<$x-E<gt>slice(':,3')> is equal to C<$x-E<gt>slice(',3')>).
 1697 
 1698 =item a
 1699 
 1700 slices only this value out of the corresponding dimension.
 1701 
 1702 =item (a)
 1703 
 1704 means the same as "a" by itself except that the resulting
 1705 dimension of length one is deleted (so if C<$x> has dims C<(3,4,5)> then
 1706 C<$x-E<gt>slice(':,(2),:')> has dimensions C<(3,5)> whereas
 1707 C<$x-E<gt>slice(':,2,:')> has dimensions C<(3,1,5))>.
 1708 
 1709 =item a:b
 1710 
 1711 slices the range I<a> to I<b> inclusive out of the dimension.
 1712 
 1713 =item a:b:c
 1714 
 1715 slices the range I<a> to I<b>, with step I<c> (i.e. C<3:7:2> gives the indices
 1716 C<(3,5,7)>). This may be confusing to Matlab users but several other
 1717 packages already use this syntax.
 1718 
 1719 
 1720 =item '*'
 1721 
 1722 inserts an extra dimension of width 1 and
 1723 
 1724 =item '*a'
 1725 
 1726 inserts an extra (dummy) dimension of width I<a>.
 1727 
 1728 =back
 1729 
 1730 An extension is planned for a later stage allowing
 1731 C<$x-E<gt>slice('(=1),(=1|5:8),3:6(=1),4:6')>
 1732 to express a multidimensional diagonal of C<$x>.
 1733 
 1734 Trivial out-of-bounds slicing is allowed: if you slice a source
 1735 dimension that doesn't exist, but only index the 0th element, then
 1736 C<slice> treats the source as if there were a dummy dimension there.
 1737 The following are all equivalent:
 1738 
 1739         xvals(5)->dummy(1,1)->slice('(2),0')  # Add dummy dim, then slice
 1740         xvals(5)->slice('(2),0')              # Out-of-bounds slice adds dim.
 1741         xvals(5)->slice((2),0)                # NiceSlice syntax
 1742         xvals(5)->((2))->dummy(0,1)           # NiceSlice syntax
 1743 
 1744 This is an error:
 1745 
 1746         xvals(5)->slice('(2),1')        # nontrivial out-of-bounds slice dies
 1747 
 1748 Because slicing doesn't directly manipulate the source and destination
 1749 pdl -- it just sets up a transformation between them -- indexing errors
 1750 often aren't reported until later.  This is either a bug or a feature,
 1751 depending on whether you prefer error-reporting clarity or speed of execution.
 1752 
 1753 =cut
 1754 
 1755 EOD
 1756         P2Child => 1,
 1757         NoPdlThread => 1,
 1758         DefaultFlow => 1,
 1759         OtherPars => 'char* str',
 1760         Comp => 'int nnew; int nthintact; int intactnew; int ndum;
 1761                  int corresp[$COMP(intactnew)]; PDL_Indx start[$COMP(intactnew)];
 1762                  PDL_Indx inc[$COMP(intactnew)]; PDL_Indx end[$COMP(intactnew)];
 1763                  int nolddims;
 1764                  int whichold[$COMP(nolddims)]; int oldind[$COMP(nolddims)];
 1765                  ',
 1766         AffinePriv => 1,
 1767         MakeComp => q~
 1768                 int i;
 1769                 int nthnew; int nthold; int nthreal;
 1770                 PDL_Indx dumsize;
 1771                 char *s; char *ns;
 1772                 int nums[3]; int nthnum;
 1773                 $COMP(nnew)=0;
 1774                 $COMP(ndum)=0;
 1775                 $COMP(nolddims) = 0;
 1776                 if(str[0] == '(')
 1777                         $COMP(nolddims)++;
 1778                 else if (str[0] == '*')
 1779                         $COMP(ndum)++;
 1780                 else if (str[0] != '\0') /* handle empty string */
 1781                         $COMP(nnew)++;
 1782                 for(i=0; str[i]; i++)
 1783                         if(str[i] == ',') {
 1784                                 if(str[i+1] == '(')
 1785                                         $COMP(nolddims)++;
 1786                                 else if(str[i+1] == '*')
 1787                                         $COMP(ndum)++;
 1788                                 else
 1789                                         $COMP(nnew)++;
 1790                         }
 1791                 $COMP(nthintact) = $COMP(nolddims) + $COMP(nnew);
 1792                 $COMP(intactnew) = $COMP(nnew)+$COMP(ndum);
 1793                 $DOCOMPDIMS();
 1794                 nthnew=0; nthold=0; i=0; nthreal=0;
 1795                 s=str-1;
 1796                 do {
 1797                         s++;
 1798                         if(isdigit(*s) || *s == '-') {
 1799                                 nthnew++; nthreal++;
 1800                                 $COMP(inc[nthnew-1]) = 1;
 1801                                 $COMP(corresp[nthnew-1]) = nthreal-1;
 1802                                 $COMP(start[nthnew-1]) = strtoll(s,&s,10);
 1803                                 if(*s != ':') {
 1804                                         $COMP(end[nthnew-1]) =
 1805                                                 $COMP(start[nthnew-1]);
 1806                                         goto outlab;
 1807                                 }
 1808                                 s++;
 1809                                 if(!isdigit(*s) && !(*s == '-')) {
 1810                                         barf("Invalid slice str ind1 '%s': '%s'",str,s);
 1811                                 }
 1812                                 $COMP(end[nthnew-1]) = strtoll(s,&s,10);
 1813                                 if(*s != ':') {goto outlab;}
 1814                                 s++;
 1815                                 if(!isdigit(*s) && !(*s == '-')) {
 1816                                         barf("Invalid slice str ind2 '%s': '%s'",str,s);
 1817                                 }
 1818                                 $COMP(inc[nthnew-1]) = strtoll(s,&s,10);
 1819                         } else switch(*s) {
 1820                         case ':':
 1821                                 s++;
 1822                                 /* FALLTHRU */
 1823                         case ',': case '\0':  /* In these cases, no inc s */
 1824                                 if ($COMP(intactnew) > 0) {
 1825                                   $COMP(start[nthnew]) = 0;
 1826                                   $COMP(end[nthnew]) = -1;
 1827                                   $COMP(inc[nthnew]) = 1;
 1828                                   $COMP(corresp[nthnew]) = nthreal;
 1829                                   nthnew++; nthreal++;
 1830                                 }
 1831                                 break;
 1832                         case '(':
 1833                                 s++;
 1834                                 $COMP(oldind[nthold]) = strtoll(s,&s,10);
 1835                                 $COMP(whichold[nthold]) = nthreal;
 1836                                 nthold++; nthreal++;
 1837                                 if(*s != ')') {
 1838                                         barf("Sliceoblit must end with ')': '%s': '%s'",str,s);
 1839                                 }
 1840                                 s++;
 1841                                 break;
 1842                         case '*':
 1843                                 s++;
 1844                                 if(isdigit(*s)) {
 1845                                         dumsize = strtoll(s,&s,10);
 1846                                 } else {dumsize = 1;}
 1847                                 $COMP(corresp[nthnew]) = -1;
 1848                                 $COMP(start[nthnew]) = 0;
 1849                                 $COMP(end[nthnew]) = dumsize-1;
 1850                                 $COMP(inc[nthnew]) = 1;
 1851                                 nthnew++;
 1852                                 break;
 1853                         }
 1854                    outlab:
 1855                         if(*s != ',' && *s != '\0') {
 1856                                 barf("Invalid slice str '%s': '%s'",str,s);
 1857                         }
 1858                 } while(*s);
 1859                 $SETREVERSIBLE(1); /* XXX Only if incs>0, no dummies */
 1860         ~,
 1861         RedoDims => '
 1862                 int i; PDL_Indx start; PDL_Indx end; PDL_Indx inc;
 1863                 if ($COMP(nthintact) > $PARENT(ndims)) {
 1864 
 1865         /* Slice has more dims than parent.  Check that the extra dims are
 1866          * all zero, and if they are then give back What You Probably Wanted,
 1867          * which is a slice with dummy dimensions of order 1 in place of each excessive
 1868          * dimension.  (Note that there are two ways to indicate a zero index: "0" and "-<w>",
 1869          * where <w> happens to be the size of that dim in the original
 1870          * piddle.  The latter case still causes an error.  That is a feature.)
 1871          *    --CED 15-March-2002
 1872          */
 1873                         int ii,parentdim,ok;
 1874                         int n_xtra_dims=0, n_xtra_olddims=0;
 1875 
 1876                            /* Check index for each extra dim in the ordinary affine list */
 1877 
 1878                         for(ok=1, ii = 0; ok && ii < $COMP(intactnew) ; ii++) {
 1879                                 parentdim = $COMP(corresp[ii]);
 1880 /*                              fprintf(stderr,"ii=%d,parent=%d, ndum=%d, nnew=%d...",ii,parentdim,$COMP(ndum),$COMP(nnew));                            */
 1881                                 if(parentdim >= $PARENT(ndims)) {
 1882 
 1883                                         ok = ( ( $COMP(start[ii]) == 0 ) &&
 1884                                                 ( $COMP(end[ii]) == 0 || $COMP(end[ii])== -1 )
 1885                                         );
 1886                                         if(ok) {
 1887                                                 /* Change this into a dummy dimension, rank 1 */
 1888                                                 $COMP(corresp[ii]) = -1;
 1889                                                 $COMP(start[ii])   = 0;
 1890                                                 $COMP(end[ii])     = 0;
 1891                                                 $COMP(inc[ii])     = 1;
 1892                                                 $COMP(ndum)++;      /* One more dummy dimension... */
 1893                                                 $COMP(nnew)--;      /* ... one less real dimension */
 1894                                                 $COMP(nthintact)--; /* ... one less intact dim */
 1895 /*                                              fprintf(stderr,"ok, ndum=%d, nnew=%d\n",$COMP(ndum), $COMP(nnew));*/
 1896                                         }
 1897 /*                              fflush(stderr);*/
 1898                                 }
 1899                         }
 1900 
 1901                           /* Check index for each indexed parent dimension */
 1902                         for(ii=0; ok && ii < $COMP(nolddims); ii++) {
 1903                                 if($COMP(whichold[ii]) >= $PARENT(ndims)) {
 1904                                         ok = ( $COMP(whichold[ii]) < $PARENT(ndims) ) ||
 1905                                                 ( $COMP(oldind[ii]) == 0 ) ||
 1906                                                 ( $COMP(oldind[ii]) == -1) ;
 1907                                         if(ok) {
 1908                                           int ij;
 1909                                           /* crunch indexed dimensions -- slow but sure */
 1910                                           $COMP(nolddims)--;
 1911                                           for(ij=ii; ij<$COMP(nolddims); ij++) {
 1912                                                 $COMP(oldind[ij]) = $COMP(oldind[ij+1]);
 1913                                                 $COMP(whichold[ij]) = $COMP(whichold[ij+1]);
 1914                                           }
 1915                                           $COMP(nthintact)--;
 1916                                         }
 1917                                 }
 1918                         }
 1919 /*      fprintf(stderr,"ok=%d\n",ok);fflush(stderr);*/
 1920                         if(ok) {
 1921                            /* Valid slice: all extra dims are zero. Adjust indices accordingly. */
 1922 /*                        $COMP(intactnew) -= $COMP(nthintact) - $PARENT(ndims); */
 1923 /*                        $COMP(nthintact) = $PARENT(ndims);*/
 1924                         } else {
 1925 
 1926                            /* Invalid slice: nonzero extra dimension.  Clean up and die.  */
 1927 
 1928                          $SETNDIMS(0); /* dirty fix */
 1929                          $PRIV(offs) = 0;
 1930                          $SETDIMS();
 1931                          $CROAK("Too many dims in slice");
 1932                         }
 1933                 }
 1934                 $SETNDIMS($PARENT(ndims)-$COMP(nthintact)+$COMP(intactnew));
 1935                 $DOPRIVDIMS();
 1936                 $PRIV(offs) = 0;
 1937                 for(i=0; i<$COMP(intactnew); i++) {
 1938                         int parentdim = $COMP(corresp[i]);
 1939                         start = $COMP(start[i]); end = $COMP(end[i]);
 1940                         inc = $COMP(inc[i]);
 1941                         if(parentdim!=-1) {
 1942                                 if(-start > $PARENT(dims[parentdim]) ||
 1943                                    -end > $PARENT(dims[parentdim])) {
 1944                                         /* set a state flag to re-trip the RedoDims code later, in
 1945                                          * case this barf is caught in an eval. This slice will
 1946                                          * always croak, so it may be smarter to find a way to
 1947                                          * replace this whole piddle with a "barf" piddle, but this
 1948                                          * will work for now. */
 1949                                         PDL->changed($CHILD_PTR(), PDL_PARENTDIMSCHANGED, 0);
 1950                                         barf("Negative slice cannot start or end above limit");
 1951                                 }
 1952                                 if(start < 0)
 1953                                         start = $PARENT(dims[parentdim]) + start;
 1954                                 if(end < 0)
 1955                                         end = $PARENT(dims[parentdim]) + end;
 1956                                 if(start >= $PARENT(dims[parentdim]) ||
 1957                                    end >= $PARENT(dims[parentdim])) {
 1958                                         /* set a state flag to re-trip the RedoDims code later, in
 1959                                          * case this barf is caught in an eval. This slice will
 1960                                          * always croak, so it may be smarter to find a way to
 1961                                          * replace this whole piddle with a "barf" piddle, but this
 1962                                          * will work for now. */
 1963                                         PDL->changed($CHILD_PTR(), PDL_PARENTDIMSCHANGED, 0);
 1964                                         barf("Slice cannot start or end above limit");
 1965                                 }
 1966                                 if(sign(end-start)*sign(inc) < 0)
 1967                                         inc = -inc;
 1968                                 $PRIV(incs[i]) = $PARENT(dimincs[parentdim]) * inc;
 1969                                 $PRIV(offs) += start * $PARENT(dimincs[parentdim]);
 1970                         } else {
 1971                                 $PRIV(incs[i]) = 0;
 1972                         }
 1973                         $CHILD(dims[i]) = ((PDL_Indx)((end-start)/inc))+1;
 1974                         if ($CHILD(dims[i]) <= 0)
 1975                            barf("slice internal error: computed slice dimension must be positive");
 1976                 }
 1977                 for(i=$COMP(nthintact); i<$PARENT(ndims); i++) {
 1978                         int cdim = i - $COMP(nthintact) + $COMP(intactnew);
 1979                         $PRIV(incs[cdim]) = $PARENT(dimincs[i]);
 1980                         $CHILD(dims[cdim]) = $PARENT(dims[i]);
 1981                 }
 1982                 for(i=0; i<$COMP(nolddims); i++) {
 1983                         int oi = $COMP(oldind[i]);
 1984                         int wo = $COMP(whichold[i]);
 1985                         if(oi < 0)
 1986                                 oi += $PARENT(dims[wo]);
 1987                         if( oi >= $PARENT(dims[wo]) )
 1988                                 $CROAK("Cannot obliterate dimension after end");
 1989                         $PRIV(offs) += $PARENT(dimincs[wo])
 1990                                         * oi;
 1991                 }
 1992         /*
 1993                 for(i=0; i<$CHILD(ndims)-$PRIV(nnew); i++) {
 1994                         $CHILD(dims[i+$COMP(intactnew)]) =
 1995                                 $PARENT(dims[i+$COMP(nthintact)]);
 1996                         $PRIV(incs[i+$COMP(intactnew)]) =
 1997                                 $PARENT(dimincs[i+$COMP(nthintact)]);
 1998                 }
 1999         */
 2000                 $SETDIMS();
 2001         ',
 2002 );
 2003 
 2004 
 2005 pp_addpm(<<'EOD'
 2006 
 2007 =head2 using
 2008 
 2009 =for ref
 2010 
 2011 Returns array of column numbers requested
 2012 
 2013 =for usage
 2014 
 2015  line $pdl->using(1,2);
 2016 
 2017 Plot, as a line, column 1 of C<$pdl> vs. column 2
 2018 
 2019 =for example
 2020 
 2021  pdl> $pdl = rcols("file");
 2022  pdl> line $pdl->using(1,2);
 2023 
 2024 =cut
 2025 
 2026 *using = \&PDL::using;
 2027 sub PDL::using {
 2028   my ($x,@ind)=@_;
 2029   @ind = list $ind[0] if (blessed($ind[0]) && $ind[0]->isa('PDL'));
 2030   foreach (@ind) {
 2031     $_ = $x->slice("($_)");
 2032   }
 2033   @ind;
 2034 }
 2035 
 2036 EOD
 2037 );
 2038 
 2039 pp_add_exported('', 'using');
 2040 
 2041 pp_addhdr(<<END
 2042 static int cmp_pdll(const void *a_,const void *b_) {
 2043         int *a = (int *)a_; int *b=(int *)b_;
 2044         if(*a>*b) return 1;
 2045         else if(*a==*b) return 0;
 2046         else return -1;
 2047 }
 2048 END
 2049 );
 2050 
 2051 
 2052 pp_def( 'affine',
 2053         P2Child => 1,
 2054         NoPdlThread => 1,
 2055         DefaultFlow => 1,
 2056         Reversible => 1,
 2057         AffinePriv => 1,
 2058         GlobalNew => 'affine_new',
 2059         OtherPars => 'PDL_Indx offspar; SV *dimlist; SV *inclist;',
 2060         Comp => 'int nd; PDL_Indx offset; PDL_Indx sdims[$COMP(nd)];
 2061                 PDL_Indx sincs[$COMP(nd)];',
 2062         MakeComp => '
 2063                 int i,n2;
 2064                 PDL_Indx *tmpi;
 2065                 PDL_Indx *tmpd = PDL->packdims(dimlist,&($COMP(nd)));
 2066                 tmpi = PDL->packdims(inclist,&n2);
 2067                 if ($COMP(nd) < 0) {
 2068                       $CROAK("Affine: can not have negative no of dims");
 2069                 }
 2070                 if ($COMP(nd) != n2)
 2071                       $CROAK("Affine: number of incs does not match dims");
 2072                 $DOCOMPDIMS();
 2073                 $COMP(offset) = offspar;
 2074                 for (i=0; i<$COMP(nd); i++) {
 2075                         $COMP(sdims)[i] = tmpd[i];
 2076                         $COMP(sincs)[i] = tmpi[i];
 2077                 }
 2078                 ',
 2079         RedoDims => '
 2080                 PDL_Indx i;
 2081                 $SETNDIMS($COMP(nd));
 2082                 $DOPRIVDIMS();
 2083                 $PRIV(offs) = $COMP(offset);
 2084                 for (i=0;i<$CHILD(ndims);i++) {
 2085                         $PRIV(incs)[i] = $COMP(sincs)[i];
 2086                         $CHILD(dims)[i] = $COMP(sdims)[i];
 2087                 }
 2088                 $SETDIMS();
 2089                 ',
 2090         Doc => undef,
 2091 );
 2092 
 2093 =head2 diagonalI
 2094 =cut
 2095 
 2096 pp_def(
 2097         'diagonalI',
 2098         P2Child => 1,
 2099         NoPdlThread => 1,
 2100         DefaultFlow => 1,
 2101         Reversible => 1,
 2102         AffinePriv => 1,
 2103         OtherPars => 'SV *list',
 2104         Comp => 'int nwhichdims; int whichdims[$COMP(nwhichdims)];',
 2105         MakeComp => '
 2106                 int i,j;
 2107                 PDL_Indx *tmp= PDL->packdims(list,&($COMP(nwhichdims)));
 2108                 if($COMP(nwhichdims) < 1) {
 2109                         $CROAK("Diagonal: must have at least 1 dimension");
 2110                 }
 2111                 $DOCOMPDIMS();
 2112                 for(i=0; i<$COMP(nwhichdims); i++)
 2113                         $COMP(whichdims)[i] = tmp[i];
 2114                 qsort($COMP(whichdims), $COMP(nwhichdims), sizeof(int),
 2115                         cmp_pdll);
 2116         ',
 2117         RedoDims => '
 2118                 int nthp,nthc,nthd; int cd = $COMP(whichdims[0]);
 2119                 $SETNDIMS($PARENT(ndims)-$COMP(nwhichdims)+1);
 2120                 $DOPRIVDIMS();
 2121                 $PRIV(offs) = 0;
 2122                 if ($COMP(whichdims)[$COMP(nwhichdims)-1] >= $PARENT(ndims) ||
 2123                         $COMP(whichdims)[0] < 0)
 2124                         $CROAK("Diagonal: dim out of range");
 2125                 nthd=0; nthc=0;
 2126                 for(nthp=0; nthp<$PARENT(ndims); nthp++)
 2127                         if (nthd < $COMP(nwhichdims) &&
 2128                             nthp == $COMP(whichdims)[nthd]) {
 2129                                 if (!nthd) {
 2130                                         $CHILD(dims)[cd] = $PARENT(dims)[cd];
 2131                                         nthc++;
 2132                                         $PRIV(incs)[cd] = 0;
 2133                                 }
 2134                                 if (nthd && $COMP(whichdims)[nthd] ==
 2135                                     $COMP(whichdims)[nthd-1])
 2136                                        $CROAK("Diagonal: dims must be unique");
 2137                                 nthd++; /* advance pointer into whichdims */
 2138                                 if($CHILD(dims)[cd] !=
 2139                                     $PARENT(dims)[nthp]) {
 2140                                         $CROAK("Different dims %d and %d",
 2141                                                 $CHILD(dims)[cd],
 2142                                                 $PARENT(dims)[nthp]);
 2143                                 }
 2144                                 $PRIV(incs)[cd] += $PARENT(dimincs)[nthp];
 2145                         } else {
 2146                                 $PRIV(incs)[nthc] = $PARENT(dimincs)[nthp];
 2147                                 $CHILD(dims)[nthc] = $PARENT(dims)[nthp];
 2148                                 nthc++;
 2149                         }
 2150                 $SETDIMS();
 2151         ',
 2152         Doc => << 'EOD',
 2153 =for ref
 2154 
 2155 Returns the multidimensional diagonal over the specified dimensions.
 2156 
 2157 The diagonal is placed at the first (by number) dimension that is
 2158 diagonalized.
 2159 The other diagonalized dimensions are removed. So if C<$x> has dimensions
 2160 C<(5,3,5,4,6,5)> then after
 2161 
 2162 =for example
 2163 
 2164  $y = $x->diagonal(0,2,5);
 2165 
 2166 the piddle C<$y> has dimensions C<(5,3,4,6)> and
 2167 C<$y-E<gt>at(2,1,0,1)> refers
 2168 to C<$x-E<gt>at(2,1,2,0,1,2)>.
 2169 
 2170 NOTE: diagonal doesn't handle threadids correctly. XXX FIX
 2171 =cut
 2172 EOD
 2173 );
 2174 
 2175 =head2 lags
 2176 =cut
 2177 
 2178 pp_def(
 2179         'lags',
 2180         Doc => <<'EOD',
 2181 =for ref
 2182 
 2183 Returns a piddle of lags to parent.
 2184 
 2185 Usage:
 2186 
 2187 =for usage
 2188 
 2189   $lags = $x->lags($nthdim,$step,$nlags);
 2190 
 2191 I.e. if C<$x> contains
 2192 
 2193  [0,1,2,3,4,5,6,7]
 2194 
 2195 then
 2196 
 2197 =for example
 2198 
 2199  $y = $x->lags(0,2,2);
 2200 
 2201 is a (5,2) matrix
 2202 
 2203  [2,3,4,5,6,7]
 2204  [0,1,2,3,4,5]
 2205 
 2206 This order of returned indices is kept because the function is
 2207 called "lags" i.e. the nth lag is n steps behind the original.
 2208 
 2209 C<$step> and C<$nlags> must be positive. C<$nthdim> can be
 2210 negative and will then be counted from the last dim backwards
 2211 in the usual way (-1 = last dim).
 2212 =cut
 2213 EOD
 2214         P2Child => 1,
 2215         NoPdlThread => 1,
 2216         DefaultFlow => 1,
 2217         Reversible => 1, # XXX Not really
 2218         AffinePriv => 1,
 2219         OtherPars => 'int nthdim; int step; int n;',
 2220         RedoDims => '
 2221                 int i;
 2222                 if ($PRIV(nthdim) < 0)  /* the usual conventions */
 2223                    $PRIV(nthdim) = $PARENT(ndims) + $PRIV(nthdim);
 2224                 if ($PRIV(nthdim) < 0 || $PRIV(nthdim) >= $PARENT(ndims))
 2225                    $CROAK("lags: dim out of range");
 2226                 if ($COMP(n) < 1)
 2227                    $CROAK("lags: number of lags must be positive");
 2228                 if ($COMP(step) < 1)
 2229                    $CROAK("lags: step must be positive");
 2230                 $PRIV(offs) = 0;
 2231                 $SETNDIMS($PARENT(ndims)+1);
 2232                 $DOPRIVDIMS();
 2233                 for(i=0; i<$PRIV(nthdim); i++) {
 2234                         $CHILD(dims)[i] = $PARENT(dims)[i];
 2235                         $PRIV(incs)[i] = $PARENT(dimincs)[i];
 2236                 }
 2237                 $CHILD(dims)[i] = $PARENT(dims)[i] - $COMP(step) * ($COMP(n)-1);
 2238                 if ($CHILD(dims)[i] < 1)
 2239                   $CROAK("lags: product of step size and "
 2240                          "number of lags too large");
 2241                 $CHILD(dims)[i+1] = $COMP(n);
 2242                 $PRIV(incs)[i] = ($PARENT(dimincs)[i]);
 2243                 $PRIV(incs)[i+1] = - $PARENT(dimincs)[i] * $COMP(step);
 2244                 $PRIV(offs) += ($CHILD(dims)[i+1] - 1) * (-$PRIV(incs)[i+1]);
 2245                 i++;
 2246                 for(; i<$PARENT(ndims); i++) {
 2247                         $CHILD(dims)[i+1] = $PARENT(dims)[i];
 2248                         $PRIV(incs)[i+1] = $PARENT(dimincs)[i];
 2249                 }
 2250                 $SETDIMS();
 2251         '
 2252 );
 2253 
 2254 =head2 splitdim
 2255 =cut
 2256 
 2257 pp_def(
 2258         'splitdim',
 2259         Doc => <<'EOD',
 2260 =for ref
 2261 
 2262 Splits a dimension in the parent piddle (opposite of L<clump|PDL::Core/clump>)
 2263 
 2264 After
 2265 
 2266 =for example
 2267 
 2268  $y = $x->splitdim(2,3);
 2269 
 2270 the expression
 2271 
 2272  $y->at(6,4,m,n,3,6) == $x->at(6,4,m+3*n)
 2273 
 2274 is always true (C<m> has to be less than 3).
 2275 =cut
 2276 EOD
 2277         P2Child => 1,
 2278         NoPdlThread => 1,
 2279         DefaultFlow => 1,
 2280         Reversible => 1, # XXX Not really
 2281         OtherPars => 'int nthdim; int nsp;',
 2282         AffinePriv => 1,
 2283         RedoDims => '
 2284                 int i = $COMP(nthdim);
 2285                 int nsp = $COMP(nsp);
 2286                 if(nsp == 0) {die("Splitdim: Cannot split to 0\n");}
 2287                 if(i <0 || i >= $PARENT(ndims)) {
 2288                         die("Splitdim: nthdim (%d) must not be negative or greater or equal to number of dims (%d)\n",
 2289                                 i, $PARENT(ndims));
 2290                 }
 2291                 if(nsp > $PARENT(dims[i])) {
 2292                         die("Splitdim: nsp (%d) cannot be greater than dim (%"IND_FLAG")\n",
 2293                                 nsp, $PARENT(dims[i]));
 2294                 }
 2295                 $PRIV(offs) = 0;
 2296                 $SETNDIMS($PARENT(ndims)+1);
 2297                 $DOPRIVDIMS();
 2298                 for(i=0; i<$PRIV(nthdim); i++) {
 2299                         $CHILD(dims)[i] = $PARENT(dims)[i];
 2300                         $PRIV(incs)[i] = $PARENT(dimincs)[i];
 2301                 }
 2302                 $CHILD(dims)[i] = $COMP(nsp);
 2303                 $CHILD(dims)[i+1] = $PARENT(dims)[i] / $COMP(nsp);
 2304                 $PRIV(incs)[i] = $PARENT(dimincs)[i];
 2305                 $PRIV(incs)[i+1] = $PARENT(dimincs)[i] * $COMP(nsp);
 2306                 i++;
 2307                 for(; i<$PARENT(ndims); i++) {
 2308                         $CHILD(dims)[i+1] = $PARENT(dims)[i];
 2309                         $PRIV(incs)[i+1] = $PARENT(dimincs)[i];
 2310                 }
 2311                 $SETDIMS();
 2312         ',
 2313 );
 2314 
 2315 =head2 rotate
 2316 =cut
 2317 
 2318 pp_def('rotate',
 2319         Doc => <<'EOD',
 2320 =for ref
 2321 
 2322 Shift vector elements along with wrap. Flows data back&forth.
 2323 =cut
 2324 EOD
 2325         Pars=>'x(n); indx shift(); [oca]y(n)',
 2326         DefaultFlow => 1,
 2327         Reversible => 1,
 2328         Code=>'
 2329         PDL_Indx i,j;
 2330         PDL_Indx n_size = $SIZE(n);
 2331         if (n_size == 0)
 2332           barf("can not shift zero size piddle (n_size is zero)");
 2333         j = ($shift()) % n_size;
 2334         if (j < 0)
 2335                 j += n_size;
 2336         for(i=0; i<n_size; i++,j++) {
 2337             if (j == n_size)
 2338                j = 0;
 2339             $y(n=>j) = $x(n=>i);
 2340         }',
 2341         BackCode=>'
 2342         PDL_Indx i,j;
 2343         PDL_Indx n_size = $SIZE(n);
 2344         j = ($shift()) % n_size;
 2345         if (j < 0)
 2346                 j += n_size;
 2347         for(i=0; i<n_size; i++,j++) {
 2348             if (j == n_size)
 2349                j = 0;
 2350             $x(n=>i) = $y(n=>j);
 2351         }
 2352         '
 2353 );
 2354 
 2355 # This is a bit tricky. Hope I haven't missed any cases.
 2356 
 2357 =head2 threadI
 2358 =cut
 2359 
 2360 pp_def(
 2361         'threadI',
 2362         Doc => <<'EOD',
 2363 =for ref
 2364 
 2365 internal
 2366 
 2367 Put some dimensions to a threadid.
 2368 
 2369 =for example
 2370 
 2371  $y = $x->threadI(0,1,5); # thread over dims 1,5 in id 1
 2372 
 2373 =cut
 2374 
 2375 EOD
 2376         P2Child => 1,
 2377         NoPdlThread => 1,
 2378         DefaultFlow => 1,
 2379         Reversible => 1,
 2380         AffinePriv => 1,
 2381         CallCopy => 0,  # Don't CallCopy for subclassed objects because PDL::Copy calls ThreadI
 2382                         #  (Wouldn't cause recursive loop otherwise)
 2383         OtherPars => 'int id; SV *list',
 2384         Comp => 'int id; int nwhichdims; int whichdims[$COMP(nwhichdims)];
 2385                         int nrealwhichdims; ',
 2386         MakeComp => '
 2387                 int i,j;
 2388                 PDL_Indx *tmp= PDL->packdims(list,&($COMP(nwhichdims)));
 2389                 $DOCOMPDIMS();
 2390                 for(i=0; i<$COMP(nwhichdims); i++)
 2391                         $COMP(whichdims)[i] = tmp[i];
 2392                 $COMP(nrealwhichdims) = 0;
 2393                 for(i=0; i<$COMP(nwhichdims); i++) {
 2394                         for(j=i+1; j<$COMP(nwhichdims); j++)
 2395                                 if($COMP(whichdims[i]) == $COMP(whichdims[j]) &&
 2396                                    $COMP(whichdims[i]) != -1) {
 2397                                 $CROAK("Thread: duplicate arg %d %d %d",
 2398                                         i,j,$COMP(whichdims[i]));
 2399                         }
 2400                         if($COMP(whichdims)[i] != -1) {
 2401                                 $COMP(nrealwhichdims) ++;
 2402                         }
 2403                 }
 2404                 $COMP(id) = id;
 2405                 ',
 2406         RedoDims => '
 2407                 int nthc,i,j,flag;
 2408                 $SETNDIMS($PARENT(ndims));
 2409                 $DOPRIVDIMS();
 2410                 $PRIV(offs) = 0;
 2411                 nthc=0;
 2412                 for(i=0; i<$PARENT(ndims); i++) {
 2413                         flag=0;
 2414                         if($PARENT(nthreadids) > $COMP(id) && $COMP(id) >= 0 &&
 2415                            i == $PARENT(threadids[$COMP(id)])) {
 2416                            nthc += $COMP(nwhichdims);
 2417                         }
 2418                         for(j=0; j<$COMP(nwhichdims); j++) {
 2419                                 if($COMP(whichdims[j] == i)) {flag=1; break;}
 2420                         }
 2421                         if(flag) {
 2422                                 continue;
 2423                         }
 2424                         $CHILD(dims[nthc]) = $PARENT(dims[i]);
 2425                         $PRIV(incs[nthc]) = $PARENT(dimincs[i]);
 2426                         nthc++;
 2427                 }
 2428                 for(i=0; i<$COMP(nwhichdims); i++) {
 2429                         int cdim,pdim;
 2430                         cdim = i +
 2431                          ($PARENT(nthreadids) > $COMP(id) && $COMP(id) >= 0?
 2432                           $PARENT(threadids[$COMP(id)]) : $PARENT(ndims))
 2433                           - $COMP(nrealwhichdims);
 2434                         pdim = $COMP(whichdims[i]);
 2435                         if(pdim == -1) {
 2436                                 $CHILD(dims[cdim]) = 1;
 2437                                 $PRIV(incs[cdim]) = 0;
 2438                         } else {
 2439                                 $CHILD(dims[cdim]) = $PARENT(dims[pdim]);
 2440                                 $PRIV(incs[cdim]) = $PARENT(dimincs[pdim]);
 2441                         }
 2442                 }
 2443                 $SETDIMS();
 2444                 PDL->reallocthreadids($CHILD_PTR(),
 2445                         ($PARENT(nthreadids)<=$COMP(id) ?
 2446                                 $COMP(id)+1 : $PARENT(nthreadids)));
 2447                 for(i=0; i<$CHILD(nthreadids); i++) {
 2448                         $CHILD(threadids[i]) =
 2449                          ($PARENT(nthreadids) > i ?
 2450                           $PARENT(threadids[i]) : $PARENT(ndims)) +
 2451                          (i <= $COMP(id) ? - $COMP(nrealwhichdims) :
 2452                           $COMP(nwhichdims) - $COMP(nrealwhichdims));
 2453                 }
 2454                 $CHILD(threadids[$CHILD(nthreadids)]) = $CHILD(ndims);
 2455                 ',
 2456 );
 2457 
 2458 
 2459 =head2 identvaff
 2460 =cut
 2461 
 2462 # we don't really need this one since it can be achieved with
 2463 # a ->threadI(-1,[])
 2464 pp_def('identvaff',
 2465         P2Child => 1,
 2466         NoPdlThread => 1,
 2467         DefaultFlow => 1,
 2468         Reversible => 1,
 2469         AffinePriv => 1,
 2470         RedoDims => '
 2471                 int i;
 2472                 $SETNDIMS($PARENT(ndims));
 2473                 $DOPRIVDIMS();
 2474                 $PRIV(offs) = 0;
 2475                 for(i=0; i<$PARENT(ndims); i++) {
 2476                         $CHILD(dims[i]) = $PARENT(dims[i]);
 2477                         $PRIV(incs[i]) = $PARENT(dimincs[i]);
 2478                 }
 2479                 $SETDIMS();
 2480                 $SETDELTATHREADIDS(0);
 2481                 $CHILD(threadids[$CHILD(nthreadids)]) = $CHILD(ndims);
 2482                 ',
 2483         Doc => <<'EOD',
 2484 =for ref
 2485 
 2486 A vaffine identity transformation (includes thread_id copying).
 2487 
 2488 Mainly for internal use.
 2489 =cut
 2490 EOD
 2491 );
 2492 
 2493 
 2494 =head2 unthread
 2495 =cut
 2496 
 2497 pp_def(
 2498         'unthread',
 2499         Doc => <<'EOD',
 2500 =for ref
 2501 
 2502 All threaded dimensions are made real again.
 2503 
 2504 See [TBD Doc] for details and examples.
 2505 =cut
 2506 EOD
 2507         P2Child => 1,
 2508         NoPdlThread => 1,
 2509         DefaultFlow => 1,
 2510         Reversible => 1,
 2511         AffinePriv => 1,
 2512         OtherPars => 'int atind;',
 2513         RedoDims => '
 2514                 int i;
 2515                 $SETNDIMS($PARENT(ndims));
 2516                 $DOPRIVDIMS();
 2517                 $PRIV(offs) = 0;
 2518                 for(i=0; i<$PARENT(ndims); i++) {
 2519                         int corc;
 2520                         if(i<$COMP(atind)) {
 2521                                 corc = i;
 2522                         } else if(i < $PARENT(threadids[0])) {
 2523                                 corc = i + $PARENT(ndims)-$PARENT(threadids[0]);
 2524                         } else {
 2525                                 corc = i - $PARENT(threadids[0]) + $COMP(atind);
 2526                         }
 2527                         $CHILD(dims[corc]) = $PARENT(dims[i]);
 2528                         $PRIV(incs[corc]) = $PARENT(dimincs[i]);
 2529                 }
 2530                 $SETDIMS();
 2531         ',
 2532 );
 2533 
 2534 
 2535 pp_add_exported('', 'dice dice_axis');
 2536 pp_addpm(<<'EOD');
 2537 
 2538 =head2 dice
 2539 
 2540 =for ref
 2541 
 2542 Dice rows/columns/planes out of a PDL using indexes for
 2543 each dimension.
 2544 
 2545 This function can be used to extract irregular subsets
 2546 along many dimension of a PDL, e.g. only certain rows in an image,
 2547 or planes in a cube. This can of course be done with
 2548 the usual dimension tricks but this saves having to
 2549 figure it out each time!
 2550 
 2551 This method is similar in functionality to the L<slice|/slice>
 2552 method, but L<slice|/slice> requires that contiguous ranges or ranges
 2553 with constant offset be extracted. ( i.e. L<slice|/slice> requires
 2554 ranges of the form C<1,2,3,4,5> or C<2,4,6,8,10>). Because of this
 2555 restriction, L<slice|/slice> is more memory efficient and slightly faster
 2556 than dice
 2557 
 2558 =for usage
 2559 
 2560  $slice = $data->dice([0,2,6],[2,1,6]); # Dicing a 2-D array
 2561 
 2562 The arguments to dice are arrays (or 1D PDLs) for each dimension
 2563 in the PDL. These arrays are used as indexes to which rows/columns/cubes,etc
 2564 to dice-out (or extract) from the C<$data> PDL.
 2565 
 2566 Use C<X> to select all indices along a given dimension (compare also
 2567 L<mslice|PDL::Core/mslice>). As usual (in slicing methods) trailing
 2568 dimensions can be omitted implying C<X>'es for those.
 2569 
 2570 =for example
 2571 
 2572  pdl> $x = sequence(10,4)
 2573  pdl> p $x
 2574  [
 2575   [ 0  1  2  3  4  5  6  7  8  9]
 2576   [10 11 12 13 14 15 16 17 18 19]
 2577   [20 21 22 23 24 25 26 27 28 29]
 2578   [30 31 32 33 34 35 36 37 38 39]
 2579  ]
 2580  pdl> p $x->dice([1,2],[0,3]) # Select columns 1,2 and rows 0,3
 2581  [
 2582   [ 1  2]
 2583   [31 32]
 2584  ]
 2585  pdl> p $x->dice(X,[0,3])
 2586  [
 2587   [ 0  1  2  3  4  5  6  7  8  9]
 2588   [30 31 32 33 34 35 36 37 38 39]
 2589  ]
 2590  pdl> p $x->dice([0,2,5])
 2591  [
 2592   [ 0  2  5]
 2593   [10 12 15]
 2594   [20 22 25]
 2595   [30 32 35]
 2596  ]
 2597 
 2598 As this is an index function, any modifications to the
 2599 slice change the parent (use the C<.=> operator).
 2600 
 2601 =cut
 2602 
 2603 sub PDL::dice {
 2604 
 2605         my $self = shift;
 2606         my @dim_indexes = @_;  # array of dimension indexes
 2607 
 2608         # Check that the number of dim indexes <=
 2609         #    number of dimensions in the PDL
 2610         my $no_indexes = scalar(@dim_indexes);
 2611         my $noDims = $self->getndims;
 2612         barf("PDL::dice: Number of index arrays ($no_indexes) not equal to the dimensions of the PDL ($noDims")
 2613                          if $no_indexes > $noDims;
 2614         my $index;
 2615         my $pdlIndex;
 2616         my $outputPDL=$self;
 2617         my $indexNo = 0;
 2618 
 2619         # Go thru each index array and dice the input PDL:
 2620         foreach $index(@dim_indexes){
 2621                 $outputPDL = $outputPDL->dice_axis($indexNo,$index)
 2622                         unless !ref $index && $index eq 'X';
 2623 
 2624                 $indexNo++;
 2625         }
 2626 
 2627         return $outputPDL;
 2628 }
 2629 *dice = \&PDL::dice;
 2630 
 2631 
 2632 =head2 dice_axis
 2633 
 2634 =for ref
 2635 
 2636 Dice rows/columns/planes from a single PDL axis (dimension)
 2637 using index along a specified axis
 2638 
 2639 This function can be used to extract irregular subsets
 2640 along any dimension, e.g. only certain rows in an image,
 2641 or planes in a cube. This can of course be done with
 2642 the usual dimension tricks but this saves having to
 2643 figure it out each time!
 2644 
 2645 =for usage
 2646 
 2647  $slice = $data->dice_axis($axis,$index);
 2648 
 2649 =for example
 2650 
 2651  pdl> $x = sequence(10,4)
 2652  pdl> $idx = pdl(1,2)
 2653  pdl> p $x->dice_axis(0,$idx) # Select columns
 2654  [
 2655   [ 1  2]
 2656   [11 12]
 2657   [21 22]
 2658   [31 32]
 2659  ]
 2660  pdl> $t = $x->dice_axis(1,$idx) # Select rows
 2661  pdl> $t.=0
 2662  pdl> p $x
 2663  [
 2664   [ 0  1  2  3  4  5  6  7  8  9]
 2665   [ 0  0  0  0  0  0  0  0  0  0]
 2666   [ 0  0  0  0  0  0  0  0  0  0]
 2667   [30 31 32 33 34 35 36 37 38 39]
 2668  ]
 2669 
 2670 The trick to using this is that the index selects
 2671 elements along the dimensions specified, so if you
 2672 have a 2D image C<axis=0> will select certain C<X> values
 2673 - i.e. extract columns
 2674 
 2675 As this is an index function, any modifications to the
 2676 slice change the parent.
 2677 
 2678 =cut
 2679 
 2680 sub PDL::dice_axis {
 2681   my($self,$axis,$idx) = @_;
 2682 
 2683   # Convert to PDLs: array refs using new, otherwise use topdl:
 2684   my $ix = (ref($idx) eq 'ARRAY') ? ref($self)->new($idx) : ref($self)->topdl($idx);
 2685   my $n = $self->getndims;
 2686   my $x = $ix->getndims;
 2687   barf("index_axis: index must be <=1D") if $x>1;
 2688   return $self->mv($axis,0)->index1d($ix)->mv(0,$axis);
 2689 }
 2690 *dice_axis = \&PDL::dice_axis;
 2691 
 2692 
 2693 EOD
 2694 
 2695 pp_addpm(<<'EOD-slice');
 2696 
 2697 =head2 slice
 2698 
 2699 =for usage
 2700 
 2701   $slice = $data->slice([2,3],'x',[2,2,0],"-1:1:-1", "*3");
 2702 
 2703 =for ref
 2704 
 2705 Extract rectangular slices of a piddle, from a string specifier,
 2706 an array ref specifier, or a combination.
 2707 
 2708 C<slice> is the main method for extracting regions of PDLs and
 2709 manipulating their dimensionality.  You can call it directly or
 2710 via he L<NiceSlice|PDL::NiceSlice> source prefilter that extends
 2711 Perl syntax o include array slicing.
 2712 
 2713 C<slice> can extract regions along each dimension of a source PDL,
 2714 subsample or reverse those regions, dice each dimension by selecting a
 2715 list of locations along it, or basic PDL indexing routine.  The
 2716 selected subfield remains connected to the original PDL via dataflow.
 2717 In most cases this neither allocates more memory nor slows down
 2718 subsequent operations on either of the two connected PDLs.
 2719 
 2720 You pass in a list of arguments.  Each term in the list controls
 2721 the disposition of one axis of the source PDL and/or returned PDL.
 2722 Each term can be a string-format cut specifier, a list ref that
 2723 gives the same information without recourse to string manipulation,
 2724 or a PDL with up to 1 dimension giving indices along that axis that
 2725 should be selected.
 2726 
 2727 If you want to pass in a single string specifier for the entire
 2728 operation, you can pass in a comma-delimited list as the first
 2729 argument.  C<slice> detects this condition and splits the string
 2730 into a regular argument list.  This calling style is fully
 2731 backwards compatible with C<slice> calls from before PDL 2.006.
 2732 
 2733 B<STRING SYNTAX>
 2734 
 2735 If a particular argument to C<slice> is a string, it is parsed as a
 2736 selection, an affine slice, or a dummy dimension depending on the
 2737 form.  Leading or trailing whitespace in any part of each specifier is
 2738 ignored (though it is not ignored within numbers).
 2739 
 2740 =over 3
 2741 
 2742 =item C<< '' >>, C<< : >>, or C<< X >> -- keep
 2743 
 2744 The empty string, C<:>, or C<X> cause the entire corresponding
 2745 dimension to be kept unchanged.
 2746 
 2747 
 2748 =item C<< <n> >> -- selection
 2749 
 2750 A single number alone causes a single index to be selected from the
 2751 corresponding dimension.  The dimension is kept (and reduced to size
 2752 1) in the output.
 2753 
 2754 =item C<< (<n>) >> -- selection and collapse
 2755 
 2756 A single number in parenthesis causes a single index to be selected
 2757 from the corresponding dimension.  The dimension is discarded
 2758 (completely eliminated) in the output.
 2759 
 2760 =item C<< <n>:<m> >> -- select an inclusive range
 2761 
 2762 Two numbers separated by a colon selects a range of values from the
 2763 corresponding axis, e.g. C<< 3:4 >> selects elements 3 and 4 along the
 2764 corresponding axis, and reduces that axis to size 2 in the output.
 2765 Both numbers are regularized so that you can address the last element
 2766 of the axis with an index of C< -1 >.  If, after regularization, the
 2767 two numbers are the same, then exactly one element gets selected (just
 2768 like the C<< <n> >> case).  If, after regulariation, the second number
 2769 is lower than the first, then the resulting slice counts down rather
 2770 than up -- e.g. C<-1:0> will return the entire axis, in reversed
 2771 order.
 2772 
 2773 =item C<< <n>:<m>:<s> >> -- select a range with explicit step
 2774 
 2775 If you include a third parameter, it is the stride of the extracted
 2776 range.  For example, C<< 0:-1:2 >> will sample every other element
 2777 across the complete dimension.  Specifying a stride of 1 prevents
 2778 autoreversal -- so to ensure that your slice is *always* forward
 2779 you can specify, e.g., C<< 2:$n:1 >>.  In that case, an "impossible"
 2780 slice gets an Empty PDL (with 0 elements along the corresponding
 2781 dimension), so you can generate an Empty PDL with a slice of the
 2782 form C<< 2:1:1 >>.
 2783 
 2784 =item C<< *<n> >> -- insert a dummy dimension
 2785 
 2786 Dummy dimensions aren't present in the original source and are
 2787 "mocked up" to match dimensional slots, by repeating the data
 2788 in the original PDL some number of times.  An asterisk followed
 2789 by a number produces a dummy dimension in the output, for
 2790 example C<< *2 >> will generate a dimension of size 2 at
 2791 the corresponding location in the output dim list.  Omitting
 2792 the number (and using just an asterisk) inserts a dummy dimension
 2793 of size 1.
 2794 
 2795 =back
 2796 
 2797 B<ARRAY REF SYNTAX>
 2798 
 2799 If you feed in an ARRAY ref as a slice term, then it can have
 2800 0-3 elements.  The first element is the start of the slice along
 2801 the corresponding dim; the second is the end; and the third is
 2802 the stepsize.  Different combinations of inputs give the same
 2803 flexibility as the string syntax.
 2804 
 2805 =over 3
 2806 
 2807 =item C<< [] >> - keep dim intact
 2808 
 2809 An empty ARRAY ref keeps the entire corresponding dim
 2810 
 2811 =item C<< [ 'X' ] >> - keep dim intact
 2812 
 2813 =item C<< [ '*',$n ] >> - generate a dummy dim of size $n
 2814 
 2815 If $n is missing, you get a dummy dim of size 1.
 2816 
 2817 =item C<< [ $dex, , 0 ] >> - collapse and discard dim
 2818 
 2819 C<$dex> must be a single value.  It is used to index
 2820 the source, and the corresponding dimension is discarded.
 2821 
 2822 =item C<< [ $start, $end ] >> - collect inclusive slice
 2823 
 2824 In the simple two-number case, you get a slice that runs
 2825 up or down (as appropriate) to connect $start and $end.
 2826 
 2827 =item C<< [ $start, $end, $inc ] >> - collect inclusive slice
 2828 
 2829 The three-number case works exactly like the three-number
 2830 string case above.
 2831 
 2832 =back
 2833 
 2834 B<PDL args for dicing>
 2835 
 2836 If you pass in a 0- or 1-D PDL as a slicing argument, the
 2837 corresponding dimension is "diced" -- you get one position
 2838 along the corresponding dim, per element of the indexing PDL,
 2839 e.g. C<< $x->slice( pdl(3,4,9)) >> gives you elements 3, 4, and
 2840 9 along the 0 dim of C<< $x >>.
 2841 
 2842 Because dicing is not an affine transformation, it is slower than
 2843 direct slicing even though the syntax is convenient.
 2844 
 2845 
 2846 =for example
 2847 
 2848  $x->slice('1:3');  #  return the second to fourth elements of $x
 2849  $x->slice('3:1');  #  reverse the above
 2850  $x->slice('-2:1'); #  return last-but-one to second elements of $x
 2851 
 2852  $x->slice([1,3]);  # Same as above three calls, but using array ref syntax
 2853  $x->slice([3,1]);
 2854  $x->slice([-2,1]);
 2855 
 2856 =cut
 2857 
 2858 
 2859 ##############################
 2860 # 'slice' is now implemented as a small Perl wrapper around
 2861 # a PP call.  This permits unification of the former slice,
 2862 # dice, and nslice into a single call.  At the moment, dicing
 2863 # is implemented a bit kludgily (it is detected in the Perl
 2864 # front-end), but it is serviceable.
 2865 #  --CED 12-Sep-2013
 2866 
 2867 *slice = \&PDL::slice;
 2868 sub PDL::slice (;@) {
 2869     my ($source, @others) = @_;
 2870 
 2871     # Deal with dicing.  This is lame and slow compared to the
 2872     # faster slicing, but works okay.  We loop over each argument,
 2873     # and if it's a PDL we dispatch it in the most straightforward
 2874     # way.  Single-element and zero-element PDLs are trivial and get
 2875     # converted into slices for faster handling later.
 2876 
 2877     for my $i(0..$#others) {
 2878       if( blessed($others[$i]) && $others[$i]->isa('PDL') ) {
 2879         my $idx = $others[$i];
 2880         if($idx->ndims > 1) {
 2881           barf("slice: dicing parameters must be at most 1D (arg $i)\n");
 2882         }
 2883         my $nlm = $idx->nelem;
 2884 
 2885         if($nlm > 1) {
 2886 
 2887        #### More than one element - we have to dice (darn it).
 2888            my $n = $source->getndims;
 2889            $source = $source->mv($i,0)->index1d($idx)->mv(0,$i);
 2890            $others[$i] = '';
 2891 
 2892         } 
 2893     elsif($nlm) {
 2894 
 2895            #### One element - convert to a regular slice.
 2896            $others[$i] = $idx->flat->at(0);
 2897 
 2898         }
 2899     else {
 2900     
 2901            #### Zero elements -- force an extended empty.
 2902            $others[$i] = "1:0:1";
 2903         }
 2904       }
 2905     }
 2906 
 2907     PDL::sliceb($source,\@others);
 2908 }
 2909 
 2910 EOD-slice
 2911 
 2912 pp_add_exported('', 'slice');
 2913 
 2914 ##########
 2915 # This is a kludge to pull arbitrary data out of a single-element PDL, using the Types.pm stuff,
 2916 # to make it easier to slice using single-element PDL arguments inside a slice specifier.
 2917 # The string $sliceb_data_kludge generates some code that physicalizes a PDL, ensures it has
 2918 # only one element, and extracts that element in a type-independent manner.  It's a pain because
 2919 # we have to generate the switch statement using information in the Config typehash.  But it saves
 2920 # time compared to parsing out any passed-in PDLs on the Perl side.
 2921 #
 2922 use PDL::Types;
 2923 $sliceb_data_kludge = <<'KLUDGE';
 2924     { pdl *p = PDL->SvPDLV( *svp );
 2925       int i;
 2926       PDL->make_physical(p);
 2927       if(p->nvals==0)
 2928         barf("slice: empty PDL in slice specifier");
 2929       if(p->nvals > 1)
 2930         barf("slice: multi-element PDL in slice specifier");
 2931       if( !(p->data) ) {
 2932          barf("slice: no data in slice specifier PDL! I give up.");
 2933       }
 2934       switch( p->datatype ) {
 2935 KLUDGE
 2936 
 2937 for my $type( PDL::Types::typesrtkeys()) {
 2938     $sliceb_data_kludge .=
 2939 "        case $type: nn = *( ($PDL::Types::typehash{$type}->{realctype} *)(p->data) ); break;\n";
 2940 }
 2941 
 2942 $sliceb_data_kludge .= <<'KLUDGE';
 2943          default: barf("Unknown PDL type in slice specifier!  This should never happen."); break;
 2944       }
 2945     }
 2946 KLUDGE
 2947 
 2948 
 2949 ##############################
 2950 # sliceb is the core of slice.  The "foo/foob" nomenclature is used to avoid the argument
 2951 # counting inherent in a direct Code section call -- "slice" is a front-end that just rolls up a
 2952 # whole variable-length argument list into a single AV reference.  
 2953 #
 2954 # I'm also too lazy to figure out how to make a PMCode section work right on a dataflow PP operator.
 2955 # -- CED
 2956 
 2957 pp_def(
 2958         'sliceb',
 2959         P2Child => 1,
 2960         NoPdlThread => 1,
 2961         DefaultFlow => 1,
 2962         OtherPars => 'SV *args;',
 2963 #
 2964 # Comp stash definitions:
 2965 #  nargs - number of args in original call
 2966 #  odim[]   - maps argno to output dim (or -1 for squished dims)
 2967 #  idim[]   - maps argno to input dim  (or -1 for squished dims)
 2968 #  odim_top - one more than the highest odim encountered
 2969 #  idim_top - one more than the highest idim encountered
 2970 #  start[]  - maps argno to start index of slice range (inclusive)
 2971 #  inc[]    - maps argno to increment of slice range
 2972 #  end[]    - maps argno to end index of slice range (inclusive)
 2973 #
 2974         Comp => 'int nargs;
 2975                  int odim[$COMP(nargs)];
 2976                  int idim[$COMP(nargs)];
 2977                  int idim_top;
 2978                  int odim_top;
 2979                  PDL_Indx start[$COMP(nargs)];
 2980                  PDL_Indx inc[$COMP(nargs)];
 2981                  PDL_Indx end[$COMP(nargs)];
 2982                  ',
 2983         AffinePriv => 1,
 2984         MakeComp => <<'SLICE-MC'
 2985                 int i;
 2986                 int idim;
 2987                 int odim;
 2988                 int imax;
 2989                 int nargs;
 2990                 AV *arglist;
 2991 
 2992                 /*** Make sure we got an array ref as input and extract its corresponding AV ***/
 2993                 if(!(  args   &&
 2994                        SvROK(args)   &&
 2995                        SvTYPE(SvRV(args))==SVt_PVAV  )){
 2996                      barf("slice requires an ARRAY ref containing zero or more arguments");
 2997                 }
 2998 
 2999                 arglist = (AV *)(SvRV(args));
 3000 
 3001                 /* Detect special case of a single comma-delimited string; in that case, */
 3002                 /* split out our own arglist.                                            */
 3003 
 3004                 if( (av_len(arglist) == 0) ) {
 3005 
 3006                   /***   single-element list: pull first element ***/
 3007 
 3008                   SV **svp;
 3009                   svp = av_fetch(arglist, 0, 0);
 3010 
 3011                   if(svp && *svp && *svp != &PL_sv_undef && SvPOKp(*svp)) {
 3012 
 3013                     /*** The element exists and is not undef and has a cached string value ***/
 3014 
 3015                     char *s,*ss;
 3016                     s = ss = SvPVbyte_nolen(*svp);
 3017                     for(;  *ss && *ss != ',';  ss++) {}
 3018 
 3019                     if(*ss == ',') {
 3020               char *s1;
 3021 
 3022                       /* the string contains at least one comma.  ATTACK!      */
 3023                       /* We make a temporary array and populate it with        */
 3024                       /* SVs containing substrings -- basically split(/\,/)... */
 3025 
 3026                       AV *al = (AV *)sv_2mortal((SV *)(newAV()));
 3027 
 3028                       do {
 3029                         for(s1=s; *s1 && *s1 != ','; s1++);
 3030 
 3031                         av_push(al, newSVpvn(s, s1-s));
 3032                         if(*s1==',')
 3033                           s = ++s1;
 3034                         else
 3035                           s = s1;
 3036                       } while(*s);
 3037 
 3038                       arglist = al;  /* al is ephemeral and will evaporate at the next perl gc */
 3039 
 3040                     } /* end of contains-comma case */
 3041                   } /* end of nontrivial single-element detection */
 3042                 }/* end of single-element detection */
 3043 
 3044 
 3045                 nargs = $COMP(nargs) = av_len( arglist ) + 1;
 3046                 $DOCOMPDIMS();
 3047 
 3048 
 3049                 /**********************************************************************/
 3050                 /**** Loop over the elements of the AV input and parse into values ****/
 3051                 /**** in the start/inc/end array                                   ****/
 3052 
 3053                 for(odim=idim=i=0; i<nargs; i++) {
 3054                    SV *this;
 3055                    SV **thisp;
 3056                    SV *sv, **svp;
 3057                    char all_flag = 0;
 3058                    char squish_flag = 0;
 3059                    char dummy_flag = 0;
 3060                    char *str;
 3061                    PDL_Indx n0 = 0;
 3062                    PDL_Indx n1 = -1;
 3063                    PDL_Indx n2 = 0;  /* n2 is inc - defaults to 0 (i.e. calc in RedoDims) */
 3064 
 3065                    thisp = av_fetch( arglist, i, 0 );
 3066                    this = (thisp  ?  *thisp  : 0 );
 3067 
 3068                    /** Keep the whole dimension if the element is undefined or missing **/
 3069                    all_flag = (  (!this)   ||   (this==&PL_sv_undef)  );
 3070 
 3071                    if(!all_flag)   {
 3072                      /*
 3073                       * Main branch -- this element is not an empty string
 3074                       */
 3075 
 3076                      if(SvROK(this)) {
 3077 
 3078                        /*** It's a reference - it better be an array ref. ***/
 3079                        int nelem;
 3080                        AV *sublist;
 3081 
 3082                        if( SvTYPE(SvRV(this)) != SVt_PVAV ) barf("slice: non-ARRAY ref in the argument list!");
 3083 
 3084 
 3085 
 3086                        /*** It *is* an array ref!  Expand it into an AV so we can read it. ***/
 3087 
 3088                        sublist = (AV *)(SvRV(this));
 3089                        if(!sublist) {
 3090                          nelem = 0;
 3091                        } else {
 3092                          nelem = av_len(sublist) + 1;
 3093                        }
 3094 
 3095                        if(nelem > 3) barf("slice: array refs can have at most 3 elements!");
 3096 
 3097 
 3098                        if(nelem==0) {      /* No elements - keep it all */
 3099                          all_flag = 1;
 3100 
 3101                        } else /* Got at least one element */{
 3102 
 3103                          /* Load the first into n0 and check for dummy or all-clear */
 3104                          /* (if element is missing use the default value already in n0) */
 3105 
 3106                          svp = av_fetch(sublist, 0, 0);
 3107                          if(svp && *svp && *svp != &PL_sv_undef) {
 3108 
 3109                 /* There is a first element.  Check if it's a PDL, then a string, then an IV */
 3110 
 3111                             if(SvROK(*svp) && sv_isa(*svp, "PDL")) {
 3112                               PDL_Indx nn;
 3113 SLICE-MC
 3114 . $sliceb_data_kludge  # Quick and dirty single-element parser (from above)
 3115 . <<'SLICE-MC'
 3116                               n0 = nn;
 3117                             } else if( SvPOKp(*svp)) {
 3118 
 3119                                /* Not a PDL but has associated string */
 3120 
 3121                                char *str = SvPVbyte_nolen(*svp);
 3122                                switch(*str) {
 3123                                  case 'X':
 3124                                       all_flag = 1; break;
 3125                                  case '*':
 3126                                       dummy_flag = 1;
 3127                                       n0 = 1;         /* n0 is 1 so 2nd field is element count */
 3128                                       n1 = 1;         /* size defaults to 1 for dummy dims */
 3129                                       n2 = 1;         /* n2 is forced to 1 so ['*',0] gives an empty */
 3130                                       break;
 3131                                  default:             /* Doesn't start with '*' or 'X' */
 3132                                       n0 = SvIV(*svp);
 3133                                       n1 = n0;         /* n1 defaults to n0 if n0 is present */
 3134                                       break;
 3135                                }
 3136                            } else /* the element has no associated string - just parse */ {
 3137                                 n0 = SvIV(*svp);
 3138                                 n1 = n0;           /* n1 defaults to n0 if n0 is present */
 3139                            }
 3140                         } /* end of defined check.  if it's undef, leave the n's at their default value. */
 3141 
 3142 
 3143                         /* Read the second element into n1 and check for alternate squish syntax */
 3144                         if( (nelem > 1) && (!all_flag) ) {
 3145                           svp = av_fetch(sublist, 1, 0);
 3146 
 3147                           if( svp && *svp && *svp != &PL_sv_undef ) {
 3148 
 3149                 if( SvROK(*svp) && sv_isa(*svp, "PDL")) {
 3150                   PDL_Indx nn;
 3151 SLICE-MC
 3152 . $sliceb_data_kludge
 3153 . <<'SLICE-MC'
 3154                               n1 = nn;
 3155 
 3156                             } else if( SvPOKp(*svp) ) {
 3157                   /* Second element has a string - make sure it's not 'X'. */
 3158                               char *str = SvPVbyte_nolen(*svp);
 3159                               if(*str == 'X') {
 3160                                 squish_flag = 1;
 3161                                 n1 = n0;
 3162                               } else {
 3163                     n1 = SvIV(*svp);
 3164                   }
 3165                             } else {
 3166                    /* Not a PDL, no string -- just get the IV */
 3167                                n1 = SvIV(*svp);
 3168                             }
 3169                           } /* If not defined, leave at the default */
 3170                         } /* End of second-element check */
 3171 
 3172 
 3173                         /*** Now try to read the third element (inc).  ***/
 3174                         if( (nelem > 2) && !(all_flag) && !(squish_flag) && !(dummy_flag) ) {
 3175                           svp = av_fetch(sublist, 2, 0);
 3176                           if( svp && *svp && *svp != &PL_sv_undef ) {
 3177 
 3178                 if(SvROK(*svp) && sv_isa(*svp, "PDL")) {
 3179                   PDL_Indx nn;
 3180 SLICE-MC
 3181 . $sliceb_data_kludge
 3182 . << 'SLICE-MC'
 3183                               n2 = nn;
 3184                 } else {
 3185                               STRLEN len;
 3186                               SvPV( *svp, len );
 3187                               if(len>0) {           /* nonzero length -> actual value given */
 3188                                 n2 = SvIV(*svp);    /* if the step is passed in as 0, it is a squish */
 3189                                 if(n2==0) {
 3190                                   n1 = n0;
 3191                                   squish_flag = 1;
 3192                                 }
 3193                               }
 3194                 }
 3195                           } /* end of nontrivial third-element parsing */
 3196                         } /* end of third-element parsing  */
 3197                        } /* end of nontrivial sublist parsing */
 3198 
 3199                      } else /* this argument is not an ARRAY ref - parse as a scalar */ {
 3200 
 3201                         /****************************************************************/
 3202                         /*** String handling part of slice is here.  Parse out each   ***/
 3203                         /*** term:                                                    ***/
 3204                         /***   <n> (or NV) - extract one element <n> from this dim    ***/
 3205                         /***   <n>:<m>     - extract <n> to <m>; autoreverse if nec.  ***/
 3206                         /***   <n>:<m>:<s> - extract <n> to <m>, stepping by <s>      ***/
 3207                         /***  (<n>)        - extract element and discard this dim     ***/
 3208                         /***  *<n>         - insert a dummy dimension of size <n>     ***/
 3209                         /***  :            - keep this dim in its entirety            ***/
 3210                         /***  X            - keep this dim in its entirety            ***/
 3211                         /****************************************************************/
 3212 
 3213                         if(SvPOKp(this)) {
 3214                           /* this argument has a cached string */
 3215 
 3216                           char *s;
 3217                           STRLEN len;
 3218                           int subargno = 0;
 3219                           int flagged = 0;
 3220                           int squish_closed = 0;
 3221                           char buf[161];
 3222                           char ii;
 3223                           s = SvPVbyte(this, len);
 3224 
 3225                           /* Very stoopid parsing - should probably make some calls to perl string utilities... */
 3226                           while(*s) {
 3227                             if( isspace( *s ) ) {
 3228                               s++;  /* ignore and loop again */
 3229                             } else {
 3230                               /* not whitespace */
 3231 
 3232                               switch(*(s++)) {
 3233                                 case '*':
 3234                                   if(flagged || subargno)
 3235                                     barf("slice: Erroneous '*' (arg %d)",i);
 3236                                   dummy_flag = flagged = 1;
 3237                                   n0 = 1;  /* default this number to 1 (size 1); '*0' yields an empty */
 3238                                   n1 = 1;  /* no second arg allowed - default to 1 so n0 is element count */
 3239                                   n2 = -1; /* -1 so we count down to n1 from n0 */
 3240                                   break;
 3241 
 3242                                 case '(':
 3243                                   if(flagged || subargno)
 3244                                     barf("slice: Erroneous '(' (arg %d)",i);
 3245                                   squish_flag = flagged = 1;
 3246 
 3247                                   break;
 3248 
 3249                                 case 'X': case 'x':
 3250                                   if(flagged || subargno > 1)
 3251                                     barf("slice: Erroneous 'X' (arg %d)",i);
 3252                                     if(subargno==0) {
 3253                                       all_flag = flagged = 1;
 3254                                     }
 3255                                     else /* subargno is 1 - squish */ {
 3256                                       squish_flag = squish_closed = flagged = 1;
 3257                                     }
 3258                                   break;
 3259 
 3260                                 case '+': case '-':
 3261                                 case '0': case '1': case '2': case '3': case '4':
 3262                                 case '5': case '6': case '7': case '8': case '9':
 3263                                  switch(subargno) {
 3264 
 3265                                    case 0: /* first arg - change default to 1 element */
 3266                                            n0 = strtoll(--s, &s, 10);
 3267                                            n1 = n0;
 3268                                            if(dummy_flag) {
 3269                                              n0 = 1;
 3270                                            }
 3271                                            break;
 3272 
 3273                                    case 1: /* second arg - parse and keep end */
 3274                                            n1 = strtoll(--s, &s, 10);
 3275                                            break;
 3276 
 3277                                    case 2: /* third arg - parse and keep inc */
 3278                                            if( squish_flag || dummy_flag ) {
 3279                                             barf("slice: erroneous third field in slice specifier (arg %d)",i);
 3280                                            }
 3281                                            n2 = strtoll(--s, &s, 10);
 3282                                            break;
 3283 
 3284                                    default: /* oops */
 3285                                      barf("slice: too many subargs in scalar slice specifier %d",i);
 3286                                      break;
 3287                                  }
 3288                                  break;
 3289 
 3290                                 case ')':
 3291                                  if( squish_closed || !squish_flag || subargno > 0) {
 3292                                   barf("nslice: erroneous ')' (arg %d)",i);
 3293                                  }
 3294                                  squish_closed = 1;
 3295                                  break;
 3296 
 3297                                 case ':':
 3298                                  if(squish_flag && !squish_closed) {
 3299                                    barf("slice: must close squishing parens (arg %d)",i);
 3300                                  }
 3301                  if( subargno == 0 ) {
 3302                    n1 = -1;   /* Set "<n>:" default to get the rest of the range */
 3303                  }
 3304                                  if( subargno > 1 ) {
 3305                                    barf("slice: too many ':'s in scalar slice specifier %d",i);
 3306                                  }
 3307                                  subargno++;
 3308                                  break;
 3309 
 3310                                 case ',':
 3311                                  barf("slice: ','  not allowed (yet) in scalar slice specifiers!");
 3312                                  break;
 3313 
 3314                                 default:
 3315                                  barf("slice: unexpected '%c' in slice specifier (arg %d)",*s,i);
 3316                                  break;
 3317                               }
 3318                             } /* end of conditional in parse loop */
 3319 
 3320                           } /* end of parse loop */
 3321 
 3322                         } else /* end of string parsing */ {
 3323 
 3324                           /* Simplest case -- there's no cached string, so it   */
 3325                           /* must be a number.  In that case it's a simple      */
 3326                           /* extraction.  Treated as a separate case for speed. */
 3327                           n0 = SvNV(this);
 3328                           n1 = SvNV(this);
 3329                           n2 = 0;
 3330                         }
 3331 
 3332                      } /* end of scalar handling */
 3333 
 3334                   } /* end of defined-element handling (!all_flag) */
 3335 
 3336                   if( (!all_flag) + (!squish_flag) + (!dummy_flag) < 2 ) {
 3337                     barf("Looks like you triggered a bug in  slice.  two flags set in dim %d",i);
 3338                   }
 3339 
 3340                   /* Force all_flag case to be a "normal" slice */
 3341                   if(all_flag) {
 3342                     n0 = 0;
 3343                     n1 = -1;
 3344                     n2 = 1;
 3345                   }
 3346 
 3347                   /* Copy parsed values into the limits */
 3348                   $COMP(start[i]) = n0;
 3349                   $COMP(end[i])   = n1;
 3350                   $COMP(inc[i])   = n2;
 3351 
 3352                   /* Deal with dimensions */
 3353                   if(squish_flag) {
 3354                     $COMP(odim[i]) = -1;
 3355                   } else {
 3356                     $COMP(odim[i]) = odim++;
 3357                   }
 3358                   if(dummy_flag) {
 3359                     $COMP(idim[i]) = -1;
 3360                   } else {
 3361                     $COMP(idim[i]) = idim++;
 3362                   }
 3363 
 3364                 } /* end of arg-parsing loop */
 3365 
 3366                 $COMP(idim_top) = idim;
 3367                 $COMP(odim_top) = odim;
 3368 
 3369                 $SETREVERSIBLE(1);
 3370 
 3371              /*** End of MakeComp for slice       */
 3372              /****************************************/
 3373 SLICE-MC
 3374            ,
 3375            RedoDims => q{
 3376                   int o_ndims_extra = 0;
 3377                   PDL_Indx i;
 3378           PDL_Indx PDIMS;
 3379 
 3380                   if( $COMP(idim_top) < $PARENT(ndims) ) {
 3381                     o_ndims_extra = $PARENT(ndims) - $COMP(idim_top);
 3382                   }
 3383 
 3384                   /* slurped dims from the arg parsing, plus any extra thread dims */
 3385                   $SETNDIMS( $COMP(odim_top) + o_ndims_extra );
 3386                   $DOPRIVDIMS();
 3387                   $PRIV(offs) = 0;  /* Offset vector to start of slice */
 3388 
 3389                   for(i=0; i<$COMP(nargs); i++) {
 3390                       PDL_Indx start, end;
 3391 
 3392                       /** Belt-and-suspenders **/
 3393                       if( ($COMP(idim[i]) < 0)  && ($COMP(odim[i]) < 0)  ) {
 3394                         PDL->changed($CHILD_PTR(), PDL_PARENTDIMSCHANGED, 0);
 3395                         barf("slice: Hmmm, both dummy and squished -- this can never happen.  I quit.");
 3396                       }
 3397 
 3398                       /** First handle dummy dims since there's no input from the parent **/
 3399                       if( $COMP(idim[i]) < 0 ) {
 3400                          /* dummy dim - offset or diminc. */
 3401                          $CHILD( dims[ $COMP(odim[i]) ] ) = $COMP(end[i]) - $COMP(start[i]) + 1;
 3402                          $PRIV( incs[ $COMP(odim[i]) ] ) = 0;
 3403                       } else {
 3404                         PDL_Indx pdsize;
 3405 
 3406                         /** This is not a dummy dim -- deal with a regular slice along it.     **/
 3407                         /** Get parent dim size for this idim, and/or allow permissive slicing **/
 3408 
 3409                         if( $COMP(idim[i]) < $PARENT(ndims)) {
 3410                           pdsize = $PARENT( dims[$COMP(idim[i])] );
 3411                         } else {
 3412                           pdsize = 1;
 3413                         }
 3414 
 3415                         start = $COMP(start[i]);
 3416             end = $COMP(end[i]);
 3417 
 3418             if( pdsize==0 && start==0 && end==-1 && $COMP(inc[i])==0 ) {
 3419               /** Trap special case: full slices of an empty dim are empty **/
 3420               $CHILD( dims[ $COMP(odim[i]) ] ) = 0;
 3421                           $PRIV( incs[$COMP(odim[i]) ] ) = 0;
 3422             } else {
 3423               
 3424               /** Regularize and bounds-check the start location **/
 3425               if(start < 0)
 3426                 start += pdsize;
 3427               if( start < 0 || start >= pdsize ) {
 3428                 PDL->changed($CHILD_PTR(), PDL_PARENTDIMSCHANGED, 0);
 3429                 if(i >= $PARENT( ndims )) {
 3430                   barf("slice: slice has too many dims (indexes dim %d; highest is %d)",i,$PARENT( ndims )-1);
 3431                 } else {
 3432               barf("slice: slice starts out of bounds in pos %d (start is %d; source dim %d runs 0 to %d)",i,start,$COMP(idim[i]),pdsize-1);
 3433                 }
 3434               }
 3435               
 3436               if( $COMP(odim[i]) < 0) {
 3437                 
 3438                 /* squished dim - just update the offset. */
 3439                 /* start is always defined and regularized if we are here here, since */
 3440                 /* both idim[i] and odim[i] can't be <0 */
 3441                 
 3442                 $PRIV(offs) += start * $PARENT( dimincs[ $COMP(idim[i]) ] );
 3443                 
 3444               } else /* normal operation */ {
 3445                 PDL_Indx siz;
 3446                 PDL_Indx inc;
 3447                 
 3448                 /** Regularize and bounds-check the end location **/
 3449                 if(end<0)
 3450                             end += pdsize;
 3451                 if( end < 0 || end >= pdsize ) {
 3452                   PDL->changed($CHILD_PTR(), PDL_PARENTDIMSCHANGED, 0);
 3453                   barf("slice: slice ends out of bounds in pos %d (end is %d; source dim %d runs 0 to %d)",i,end,$COMP(idim[i]),pdsize-1);
 3454                 }
 3455                 
 3456                 /* regularize inc */
 3457                 
 3458                 inc = $COMP(inc[i]);
 3459                 if(!inc)
 3460                   inc = (start <= end) ? 1 : -1;
 3461                 
 3462                 siz = (end - start + inc) / inc ;
 3463                 if(siz<0)
 3464                   siz=0;
 3465                 $CHILD( dims[ $COMP(odim[i]) ] ) = siz;
 3466                 $PRIV(  incs[ $COMP(odim[i]) ] ) = $PARENT( dimincs[ $COMP(idim[i]) ] ) * inc;
 3467                 $PRIV(offs) += start * $PARENT( dimincs[ $COMP(idim[i]) ] );
 3468               } /* end of normal slice case */
 3469             } /* end of trapped strange slice case */
 3470                       } /* end of non-dummy slice case */
 3471                   } /* end of nargs loop */
 3472 
 3473                   /* Now fill in thread dimensions as needed.  idim and odim persist from the parsing loop */
 3474                   /* up above. */
 3475                   for(i=0; i<o_ndims_extra; i++) {
 3476                     $CHILD( dims   [ $COMP(odim_top) + i ] ) = $PARENT( dims   [ $COMP(idim_top) + i ] );
 3477                     $PRIV( incs[ $COMP(odim_top) + i ] ) = $PARENT( dimincs[ $COMP(idim_top) + i ] );
 3478                   }
 3479 
 3480                 $SETDIMS();
 3481 
 3482         } # end of RedoDims for slice
 3483 );
 3484 
 3485 
 3486 pp_addpm({At => 'Bot'},<< 'EOD');
 3487 
 3488 =head1 BUGS
 3489 
 3490 For the moment, you can't slice one of the zero-length dims of an
 3491 empty piddle.  It is not clear how to implement this in a way that makes
 3492 sense.
 3493 
 3494 Many types of index errors are reported far from the indexing
 3495 operation that caused them.  This is caused by the underlying architecture:
 3496 slice() sets up a mapping between variables, but that mapping isn't
 3497 tested for correctness until it is used (potentially much later).
 3498 
 3499 =head1 AUTHOR
 3500 
 3501 Copyright (C) 1997 Tuomas J. Lukka.  Contributions by
 3502 Craig DeForest, deforest@boulder.swri.edu.
 3503 Documentation contributions by David Mertens.
 3504 All rights reserved. There is no warranty. You are allowed
 3505 to redistribute this software / documentation under certain
 3506 conditions. For details, see the file COPYING in the PDL
 3507 distribution. If this file is separated from the PDL distribution,
 3508 the copyright notice should be included in the file.
 3509 
 3510 =cut
 3511 
 3512 EOD
 3513 
 3514 
 3515 pp_done();