pp_addpm({At => 'Top'},<< 'EOD');

=head1 NAME

PDL::Slices -- Stupid index tricks

=head1 SYNOPSIS

  use PDL;
  $a = ones(3,3);
  $b = $a->slice('-1:0,(1)');
  $c = $a->dummy(2);


=head1 DESCRIPTION

This package provides many of the powerful PerlDL core index
manipulation routines. These routines mostly contain two-way data flow 
so you can get a unit matrix by

 $a = zeroes(1000,1000);
 $a->diagonal(0,1) ++;

which is quite efficient. See L<PDL::Indexing> and L<PDL::Tips> for
more examples.

These functions are usually two-way:

 $b = $a->slice("1:3");
 $b += 5;   		# $a is changed!

If you want to force a copy and no "flow" backwards, you need

 $b = $a->slice("1:3")->copy;
 $b += 5;		# $a is not changed.

alternatively, you can use

 $b = $a->slice("1:3")->sever;

which does not copy the struct but beware that after

 $b = $a->slice("1:3");
 $c = $b->sever;

the variables C<$b> and C<$c> point to the same object but with 
C<-E<gt>copy> they do not.

The fact that there is this kind of flow makes PDL a very
powerful language in many ways: since you can alter the original
data by altering some easier-to-use representation of it, many things are
much easier to accomplish, just like making the above unit matrix.

Slicing is so central to the PDL language that a special compile-time
syntax has been introduced to handle it compactly; see L<PDL::NiceSlice> 
for details.

=head1 BUGS 

For the moment, you can't slice the empty piddle.  This should probably
change:  slices of the empty piddle should probably return the empty piddle.

Many types of index errors are reported far from the indexing
operation that caused them.  This is caused by the underlying architecture:
slice() sets up a mapping between variables, but that mapping isn't 
tested for correctness until it is used (potentially much later).

=cut

use PDL::Core ':Internal';

EOD

pp_addpm({At => 'Bot'},<< 'EOD');

=head1 AUTHOR

Copyright (C) 1997 Tuomas J. Lukka.  Contributions by
Craig DeForest, deforest@boulder.swri.edu.
All rights reserved. There is no warranty. You are allowed
to redistribute this software / documentation under certain
conditions. For details, see the file COPYING in the PDL
distribution. If this file is separated from the PDL distribution,
the copyright notice should be included in the file.

=cut

EOD

# $::PP_VERBOSE=1;

pp_add_boot(
"  PDL->readdata_affine = pdl_readdata_affineinternal;\n" . 
"     PDL->writebackdata_affine = pdl_writebackdata_affineinternal;\n"
);

## Several routines use the 'Dims' and 'ParentInds'
## rules - these currently do nothing

pp_def(
       'affineinternal',
       HandleBad => 1,
       AffinePriv => 1,
       DefaultFlow => 1,
       P2Child => 1,
       ReadDataFuncName => "pdl_readdata_affineinternal",
       WriteBackDataFuncName => "pdl_writebackdata_affineinternal",
       MakeComp => '$CROAK("AFMC MUSTNT BE CALLED");',
       RedoDims => '$CROAK("AFRD MUSTNT BE CALLED");',
       EquivCPOffsCode => '
		int i; int poffs=$PRIV(offs); int nd;
		for(i=0; i<$CHILD_P(nvals); i++) {
			$EQUIVCPOFFS(i,poffs);
			for(nd=0; nd<$CHILD_P(ndims); nd++) {
				poffs += $PRIV(incs[nd]);
				if(nd<$CHILD_P(ndims)-1 &&
				   (i+1)%$CHILD_P(dimincs[nd+1]) ||
				   nd == $CHILD_P(ndims)-1)
					break;
				poffs -= $PRIV(incs[nd]) *
					$CHILD_P(dims[nd]);
			}
                }',
       Doc => undef,    # 'internal',
);

pp_def(
	's_identity',
        HandleBad => 1,
	P2Child => 1,
	DefaultFlow => 1,
	OtherPars => '',
	Reversible => 1,
	Dims => '$COPYDIMS();',
	ParentInds => '$COPYINDS();',
	Identity => 1,
	Doc => 
'=for ref

Internal vaffine identity function.

=cut
',
);

$doc = <<'EOD';
=for ref

C<index> and C<index2d> provide rudimentary index indirection.

=for example

 $c = index($source,$ind);
 $c = index2d($source2,$ind1,$ind2);

use the C<$ind> variables as indices to look up values in C<$source>.  
C<index2d> uses separate piddles for X and Y coordinates.  For more
general N-dimensional indexing, see L<PDL::Slices> or the L<PDL::NiceSlice>
syntax.

These functions are two-way, i.e. after

 $c = $a->index(pdl[0,5,8]);
 $c .= pdl [0,2,4];

the changes in C<$c> will flow back to C<$a>.

C<index> provids simple threading:  multiple-dimensioned arrays are treated
as collections of 1-D arrays, so that

 $a = xvals(10,10)+10*yvals(10,10);
 $b = $a->index(3);
 $c = $a->index(9-xvals(10));

puts a single column from C<$a> into C<$b>, and puts a single element
from each column of C<$a> into C<$c>.  If you want to extract multiple
columns from an array in one operation, see L<dice|/dice> or
L<indexND|/indexND>.

=cut
EOD

my $index_init_good = 
       'register int foo = $ind(); 
        if( foo<0 || foo>=$SIZE(n) ) {
	   barf("PDL::index: invalid index %d (valid range 0..%d)",
                foo,$SIZE(n)-1);
        }';
my $index_init_bad = 
       'register int foo = $ind(); 
        if( $ISBADVAR(foo,ind) || foo<0 || foo>=$SIZE(n) ) {
	   barf("PDL::index: invalid index %d (valid range 0..%d)",
                foo,$SIZE(n)-1);
        }';

pp_def(
       'index',
       HandleBad => 1,
       DefaultFlow => 1,
       Reversible => 1,	
       Pars => 'a(n); int ind(); [oca] c();',
       Code => 
       $index_init_good . ' $c() = $a(n => foo);',
       BadCode => 
       $index_init_bad . ' $c() = $a(n => foo);',
       BackCode => 
       $index_init_good . ' $a(n => foo) = $c();',
       BadBackCode => 
       $index_init_bad . ' $a(n => foo) = $c();',
       Doc => $doc,
       BadDoc =>
       'index barfs if any of the index values are bad.',
       );

my $index2d_init_good =
       'register int fooa,foob;
	fooa = $inda(); 
        if( fooa<0 || fooa>=$SIZE(na) ) {
           barf("PDL::index: invalid x-index %d (valid range 0..%d)",
                fooa,$SIZE(na)-1);
        }
	foob = $indb(); 
        if( foob<0 || foob>=$SIZE(nb) ) {
	   barf("PDL::index: invalid y-index %d (valid range 0..%d)",
                foob,$SIZE(nb)-1);
        }';
my $index2d_init_bad =
       'register int fooa,foob;
	fooa = $inda(); 
        if( $ISBADVAR(fooa,inda) || fooa<0 || fooa>=$SIZE(na) ) {
           barf("PDL::index: invalid index 1");
        }
	foob = $indb(); 
        if( $ISBADVAR(foob,indb) || foob<0 || foob>=$SIZE(nb) ) {
	   barf("PDL::index: invalid index 2");
        }';

pp_def(
       'index2d',
       HandleBad => 1,
       DefaultFlow => 1,
       Reversible => 1,	
       Pars => 'a(na,nb); int inda(); int indb(); [oca] c();',
       Code => 
       $index2d_init_good . ' $c() = $a(na => fooa, nb => foob);',
       BadCode => 
       $index2d_init_bad . '$c() = $a(na => fooa, nb => foob);',
       BackCode => 
       $index2d_init_good . ' $a(na => fooa, nb => foob) = $c();',
       BadBackCode => 
       $index2d_init_bad . '$a(na => fooa, nb => foob) = $c();',
       Doc => $doc,
       BadDoc =>
       'index2d barfs if either of the index values are bad.',
);


# indexND: CED 2-Aug-2002
pp_add_exported('','indexND indexNDb');
pp_addpm(<<'EOD-indexND');

=head2 indexNDb

=for ref 

  Backwards-compatibility alias for indexND

=head2 indexND

=for ref

  Find selected elements in an N-D piddle, with optional boundary handling

=for example

  $out = $source->indexND( $index, [$method] )

  $source = 10*xvals(10,10) + yvals(10,10);
  $index  = pdl([[2,3],[4,5]],[[6,7],[8,9]]);
  print $source->indexND( $index );

  [
   [23 45]
   [67 89]
  ]

IndexND collapses C<$index> by lookup into C<$source>.  The 
0th dimension of C<$index> is treated as coordinates in C<$source>, and 
the return value has the same dimensions as the rest of C<$index>.
The returned elements are looked up from C<$source>.  Dataflow
works -- propagated assignment flows back into C<$source>.  

IndexND and IndexNDb were originally separate routines but they are both
now implemented as a call to L<range|/range>, and have identical syntax to
one another.

=cut

sub PDL::indexND {
	my($source,$index, $boundary) = @_;
	return PDL::range($source,$index,undef,$boundary);
}

*PDL::indexNDb = \&PDL::indexND;

EOD-indexND

pp_addpm(<<'EOD-range');

sub PDL::range {
  my($source,$ind,$sz,$bound) = @_;
  my $index = PDL->pdl($ind);
  $index = $index->dummy(0,1) unless $index->ndims;
  my $size = defined($sz) ? PDL->pdl($sz) : undef;

	
  # Pack boundary string if necessary
  if(defined $bound) {
    if(ref $bound eq 'ARRAY') {
      my ($s,$el);
      foreach $el(@$bound) {
	barf "Illegal boundary value '$el' in range"
  	  unless( $el =~ m/^([0123fFtTeEpPmM])/ );
	$s .= $1;
      }
      $bound = $s;
    }
    elsif($bound !~ m/^[0123ftepx]+$/  && $bound =~ m/^([0123ftepx])/i ) {
      $bound = $1;
    }
  }

  no warnings; # shut up about passing undef into rangeb
  $source->rangeb($index,$size,$bound);
}
EOD-range

pp_def(
	'rangeb',
	Doc => <<'EOD',
=for ref

Engine for L<range|/range>

=for example

Same calling convention as L<range|/range>, but you must supply all
parameters.  C<rangeb> is marginally faster as it makes a direct PP call,
avoiding the perl argument-parsing step. 

=cut

=head2 range

=for ref 

Extract selected chunks from a source piddle, with boundary conditions

=for example

	$out = $source->range($index,[$size,[$boundary]])

Returns elements or rectangular slices of the original piddle, indexed by 
the C<$index> piddle.  C<$source> is an N-dimensional piddle, and C<$index> is 
a piddle whose first dimension has order up to N.  Each row of C<$index> is 
treated as coordinates of a single value or chunk from C<$source>, specifying 
the location(s) to extract.  

B<INPUTS>

C<$index> and C<$size> can be piddles or array refs such as you would
feed to L<zeroes|PDL::Core/zeroes> and its ilk.  If C<$index>'s 0th dimension
has order higher than the number of dimensions in C<$source>, then
C<$source> is treated as though it had trivial dummy dimensions of
order 1, up to the required size to be indexed by C<$index> -- so if
your source array is 1-D and your index array is a list of 3-vectors, 
you get two dummy dimensions of order 1 on the end of your source array.

You can extract single elements or N-D rectangular ranges from C<$source>,
by setting C<$size>.  If C<$size> is undef or zero, then you get a single
sample for each row of C<$index>.  This behavior is similar to 
L<indexNDb|/indexNDb>, which is in fact implemented as a call to L<range|/range>.

If C<$size> is positive then you get a range of values from C<$source> at 
each location, and the output has extra dimensions allocated for them.
C<$size> can be a scalar, in which case it applies to all dimensions, or an 
N-vector, in which case each element is applied independently to the 
corresponding dimension in C<$source>.  See below for details.

C<$boundary> is a number, string, or list ref indicating the type of 
boundary conditions to use when ranges reach the edge of C<$source>.  If you
specify no boundary conditions the default is to forbid boundary violations
on all axes.  If you specify exactly one boundary condition, it applies to 
all axes.  If you specify more (as elements of a list ref, or as a packed
string, see below), then they apply to dimensions in the order in which they
appear, and the last one applies to all subsequent dimensions.  (This is 
less difficult than it sounds; see the examples below).

=over 3

=item 0 (synonyms: 'f','forbid') B<(default)>

Ranges are not allowed to cross the boundary of the original PDL.  Disallowed
ranges throw an error.  The errors are thrown at evaluation time, not
at the time of the range call (this is the same behavior as L<slice|/slice>).

=item 1 (synonyms: 't','truncate')

Values outside the original piddle get BAD if you've got bad value
support compiled into your PDL; or 0 if you haven't.  Reverse dataflow
works OK for the portion of the child that is in-bounds.  The
out-of-bounds part of the child is reset to (BAD|0) during each
dataflow operation, but execution continues.

=item 2 (synonyms: 'e','x','extend')

Values that would be outside the original piddle point instead to the 
nearest allowed value within the piddle.  See the CAVEAT below on 
mappings that are not single valued.

=item 3 (synonyms: 'p','periodic')

Periodic boundary conditions apply: the numbers in $index are applied,
strict-modulo the corresponding dimensions of $source.  This is equivalent to
duplicating the $source piddle throughout N-D space.  See the CAVEAT below 
about mappings that are not single valued.

=item 4 (synonyms: 'm','mirror')

Mirror-reflection periodic boundary conditions apply.  See the CAVEAT
below about mappings that are not single valued.

=back

The boundary condition identifiers all begin with unique characters, so 
you can feed in multiple boundary conditions as either a list ref or a 
packed string.  (The packed string is marginally faster to run).  For 
example, the four expressions [0,1], ['forbid','truncate'], ['f','t'], 
and 'ft' all specify that violating the boundary in the 0th dimension 
throws an error, and all other dimensions get truncated.  

If you feed in a single string, it is interpreted as a packed boundary
array if all of its characters are valid boundary specifiers (e.g. 'pet'),
but as a single word-style specifier if they are not (e.g. 'forbid').

B<OUTPUT>

The output threads over both C<$index> and C<$source>.  Because implicit 
threading can happen in a couple of ways, a little thought is needed.  The
returned dimension list is stacked up like this:

   (index thread dims), (index dims (size)), (source thread dims)

The first few dims of the output correspond to the extra dims of
C<$index> (beyond the 0 dim). They allow you to pick out individual
ranges from a large, threaded collection.

The middle few dims of the output correspond to the size dims
specified in C<$size>, and contain the range of values that is extracted
at each location in C<$source>.  Every nonzero element of C<$size> is copied to
the dimension list here, so that if you feed in (for example) C<$size
= [2,0,1]> you get an index dim list of C<(2,1)>.

The last few dims of the output correspond to extra dims of C<$source> beyond
the number of dims indexed by C<$index>.  These dims act like ordinary 
thread dims, because adding more dims to C<$source> just tacks extra dims
on the end of the output.  Each source thread dim ranges over the entire
corresponding dim of C<$source>.

B<Dataflow>: Dataflow is bidirectional.

B<Examples>:
Here are basic examples of C<range> operation, showing how to get 
ranges out of a small matrix.  The first few examples show extraction
and selection of individual chunks.  The last example shows
how to mark loci in the original matrix (using dataflow).

 perldl> $src = 10*xvals(10,5)+yvals(10,5)
 perldl> print $src->range([2,3])    # Cut out a single element
 23
 perldl> print $src->range([2,3],1)  # Cut out a single 1x1 block
 [
  [23]
 ]
 perldl> print $src->range([2,3], [2,1]) # Cut a 2x1 chunk
 [
  [23 33]
 ]
 perldl> print $src->range([[2,3]],[2,1]) # Trivial list of 1 chunk
 [
  [
   [23]
   [33]
  ]
 ]
 perldl> print $src->range([[2,3],[0,1]], [2,1])   # two 2x1 chunks
 [
  [
   [23  1]
   [33 11]
  ]
 ]
 perldl> # A 2x2 collection of 2x1 chunks
 perldl> print $src->range([[[1,1],[2,2]],[[2,3],[0,1]]],[2,1])
 [
  [
   [
    [11 22]
    [23  1]
   ]
   [
    [21 32]
    [33 11]
   ]
  ]
 ]
 perldl> $src = xvals(5,3)*10+yvals(5,3)
 perldl> print $src->range(3,1)  # Thread over y dimension in $src 
 [
  [30]
  [31]
  [32]
 ]

 perldl> $src = zeroes(5,4);
 perldl> $src->range(pdl([2,3],[0,1]),pdl(2,1)) .= xvals(2,2,1) + 1
 perldl> print $src
 [
  [0 0 0 0 0]
  [2 2 0 0 0]
  [0 0 0 0 0]
  [0 0 1 1 0]
 ]

B<CAVEAT>: It's quite possible to select multiple ranges that
intersect.  In that case, modifying the ranges doesn't have a
guaranteed result in the original PDL -- the result is an arbitrary
choice among the valid values.  For some things that's OK; but for
others it's not. In particular, this doesn't work:

    perldl> $photon_list = new PDL::RandVar->sample(500)->reshape(2,250)*10
    perldl> histogram = zeroes(10,10)
    perldl> histogram->range($photon_list,1)++;  #not what you wanted

The reason is that if two photons land in the same bin, then that bin 
doesn't get incremented twice.  (That may get fixed in a later version...)

B<PERMISSIVE RANGING>: If C<$index> has too many dimensions compared
to C<$source>, then $source is treated as though it had dummy
dimensions of order 1, up to the required number of dimensions.  These
virtual dummy dimensions have the usual boundary conditions applied to
them.

B<Efficiency>: Because C<range> isn't an affine transformation (it
involves lookup into a list of N-D indices), it is somewhat
memory-inefficient for long lists of ranges, and keeping dataflow open
is somewhat slower than for affine transformations (which don't have
to copy data around).  Thus, if you have a few very large (order of
10^4 element) chunks, you might find that perl lists of ordinary
slices are faster.

C<range> is a perl front-end to a PP function, C<rangeb>.  Calling
C<rangeb> is marginally faster but requires that you include all arguments.

DEVEL NOTES

* index thread dimensions are effectively clumped internally.  This
makes it easier to loop over the index array but a little more brain-bending 
to tease out the algorithm.

* Currently the index threads really do run fastest in memory; this is 
probably the wrong direction to thread, for fastest behavior -- modifying 
the appropriate dimincs in RedoDims ought to take care of it.

=cut 

EOD
	HandleBad => 1,
	DefaultFlow => 1,
	Reversible => 1,
	P2Child => 1,
	OtherPars => 'SV *index; SV *size; SV *boundary',


#
# rdim: dimensionality of each range.
#
# ntsize: number of nonzero size dimensions
# sizes:  array of range sizes, indexed (0..rdim-1).  A zero element means
#	  that the dimension is omitted from the child dim list.
# corners: parent coordinates of each corner, running fastest over coord index.
#	(indexed 0 .. (nitems-1)*(rdim)+rdim-1)
# nitems: total number of list elements   (product of itdims)
# itdim:  number of index thread dimensions
# itdims: Size of each index thread dimension, indexed (0..itdim-1)
# 
# bsize: Number of independently specified boundary conditions
# boundary: Array containing all the boundary condition specs

	Comp => 'long rdim;
                 long nitems;
                 long itdim;
                 long ntsize;
                 long bsize;
		 long sizes[$COMP(rdim)];
		 long itdims[$COMP(itdim)];
		 long corners[$COMP(rdim) * $COMP(nitems)];
	         char boundary[$COMP(rdim)];
		 ',
	MakeComp => <<'EOD-MakeComp',
pdl *ind_pdl;
pdl *size_pdl;

/***
 * Check and condition the index piddle.  Some of this is apparently
 * done by XS -- but XS doesn't check for existing SVs that are undef.
 */
if ((index==NULL) || (index == &PL_sv_undef))
   { $CROAK("index variable must be defined in range"); }

if(!(ind_pdl = PDL->SvPDLV(index))) /* assignment */
   { $CROAK("unable to convert index to a PDL in range"); }
  
PDL->make_physdims(ind_pdl);


/*
 * if(ind_pdl->nvals == 0)
 *  { $CROAK("Cannot use the empty PDL as an index in range"); }
 */

/***
 * Ensure that the index is a long.  If there's no loss of information,
 * just upgrade it -- otherwise, make a temporary copy.
 */
switch(ind_pdl->datatype) {
 default:                              /* Most types: */
   ind_pdl = PDL->hard_copy(ind_pdl);  /*   copy and fall through */
 case PDL_B: case PDL_S: case PDL_US:  
   PDL->converttype(&ind_pdl,PDL_L,1); /* convert in place. */
   break;
 case PDL_L:                           
   break;
}

/***
 * Figure sizes of the COMP arrrays and allocate them.
 */
{
  long i,nitems;

  $COMP(rdim) = ind_pdl->ndims ? ind_pdl->dims[0] : 1;
  for(i=nitems=1; i < ind_pdl->ndims; i++)  /* Accumulate item list size */
    nitems *= ind_pdl->dims[i];
  $COMP(nitems) = nitems;
  $COMP(itdim) = ind_pdl->ndims ? ind_pdl->ndims - 1 : 0;
  $DOCOMPDIMS();
}

/***
 * Fill in the boundary condition array
 */
{
  char *bstr;
  STRLEN blen;
  bstr = SvPV(boundary,blen);
  
  if(blen == 0) {
    /* If no boundary is specified then every dim gets truncated */
    int i;
    for (i=0;i<$COMP(rdim);i++) 
      $COMP(boundary[i]) = 0;
  } else {
    int i;
    for(i=0;i<$COMP(rdim);i++) {
      switch(bstr[i < blen ? i : blen-1 ]) {
      case '0': case 'f': case 'F':               /* forbid */
	$COMP(boundary[i]) = 0;
	break;
      case '1': case 't': case 'T':               /* truncate */
	$COMP(boundary[i]) = 1;
	break;
      case '2': case 'e': case 'E':               /* extend */
	$COMP(boundary[i]) = 2;
	break;
      case '3': case 'p': case 'P':               /* periodic */
	$COMP(boundary[i]) = 3;
	break;
      case '4': case 'm': case 'M':               /* mirror */
	$COMP(boundary[i]) = 4;
	break;
      default:
	{
	  char buf[BUFSIZ];
	  /* No need to check if i < blen -- this will barf out the
	   * first time it gets hit.  I didn't use $ CROAK 'coz that
	   * macro doesn't let you pass in a string variable -- only a 
	   * constant.
	   */
	  sprintf(buf,"Error in rangeb: Unknown boundary condition '%c' in range",bstr[i]);
	  barf(buf);
	}
	break;
      }
    }
  }
}  
/***
 * Store the sizes of the index-thread dims
 */
{
  long i;
  long nd = ind_pdl->ndims - 1;
  for(i=0; i < nd ; i++)
    $COMP(itdims[i]) = ind_pdl->dims[i+1];
}

/*** 
 * Check and condition the size piddle, and store sizes of the ranges
 */
{
  long i,ntsize;

  if( (size == NULL) || (size == &PL_sv_undef) ) 
    for(i=0;i<$COMP(rdim);i++)
	  $COMP(sizes[i]) = 0;
  else {
    if(!(size_pdl = PDL->SvPDLV(size))) /* assignment */
      $CROAK("Unable to convert size to a PDL in range");
    if(size_pdl->nvals == 0)
      for(i=0;i<$COMP(rdim);i++)
	$COMP(sizes[i]) = 0;
    else {
      switch(size_pdl->datatype) {
      default:                              /* Most types: */
	size_pdl = PDL->hard_copy(size_pdl);  /*   copy and fall through */
      case PDL_B: case PDL_S: case PDL_US:  
	PDL->converttype(&size_pdl,PDL_L,1); /* convert in place. */
	break;
      case PDL_L:                           
	break;
      }
      
      /* Copy the sizes, or die if they're the wrong shape */
      if(size_pdl->nvals == 1) 
	for(i=0;i<$COMP(rdim);i++) 
	  $COMP(sizes[i]) = *((long *)(size_pdl->data));

      else if( size_pdl->nvals <= $COMP(rdim) &&size_pdl->ndims == 1) 
	for(i=0;i<$COMP(rdim);i++)
	  $COMP(sizes[i]) = ( 	(i < size_pdl->nvals) ? 
	                    	((long *)(size_pdl->data))[i] :
		                0
			    );
      else {
	$CROAK("Size must match index's 0th dim in range");
      }
	

    } /* end of nonempty size-piddle code */
  } /* end of defined-size-piddle code */

  /* Insert the number of nontrivial sizes (these get output dimensions) */
  for(i=ntsize=0;i<$COMP(rdim);i++)
    if($COMP(sizes[i]))
      ntsize++;
  $COMP(ntsize) = ntsize;
}

/***
 * Stash coordinates of the corners
 */

{
  long i,j,k,ioff;
  long *cptr;
  long *iter = (long *)(PDL->smalloc(sizeof(long) * ($COMP(itdim))));
  
  /* initialize iterator to loop over index threads */
  cptr = iter;
  for(k=0;k<$COMP(itdim);k++)
    *(cptr++) = 0;
  
  cptr = $COMP(corners);  
  do {
    /* accumulate offset into the index from the iterator */
    for(k=ioff=0;k<$COMP(itdim);k++)
      ioff += iter[k] * ind_pdl->dimincs[k+1];
    
    /* Loop over the 0th dim of index, copying coords. */ 
    /* This is the natural place to check for permissive ranging; too */
    /* bad we don't have access to the parent piddle here... */
	
    for(j=0;j<$COMP(rdim);j++) 
        *(cptr++) = ((long *)(ind_pdl->data))[ioff + ind_pdl->dimincs[0] * j];
    
    /* Increment the iterator -- the test increments, the body carries. */
    for(k=0; k<$COMP(itdim) && (++(iter[k]))>=($COMP(itdims)[k]) ;k++)
      iter[k] = 0;
  } while(k<$COMP(itdim));
}

$SETREVERSIBLE(1);

EOD-MakeComp

	RedoDims => <<'EOD-RedoDims' ,
{		      
  long stdim = $PARENT(ndims) - $COMP(rdim);
  long dim,inc;
  long i;

  if(stdim < 0) 
	stdim = 0;

  /* Set dimensionality of child */
  $CHILD(ndims) = $COMP(itdim) + $COMP(ntsize) + stdim;
  $SETNDIMS($COMP(itdim)+$COMP(ntsize)+stdim);

  /* Copy index thread dimensions to child */
  inc = 1;
  for(dim=0; dim<$COMP(itdim); dim++) {
    $CHILD(dimincs[dim]) = inc;
    inc *= ($CHILD(dims[dim]) = $COMP(itdims[dim])); /* assignment */

  }
  
  /* Copy size dimensions to child, crunching as we go. */
  for(i=0;i<$COMP(rdim);i++) {
    if($COMP(sizes[i])) {
      $CHILD(dimincs[dim]) = inc;
      inc *= ($CHILD(dims[dim++]) = $COMP(sizes[i])); /* assignment */
    }
  }

  /* Copy source thread dimensions to child */
  for(i=0;i<stdim;i++) {
    $CHILD(dimincs[dim]) = inc;
    inc *= ($CHILD(dims[dim++]) = $PARENT(dims[i+$COMP(rdim)])); /* assignment */
  }


  $CHILD(datatype) = $PARENT(datatype);
  $SETDIMS();
}

EOD-RedoDims

	EquivCPOffsCode => <<'EOD-EquivCPOffsCode',
{
  long *iter, *ip;  /* vector iterator */
  long *sizes, *sp; /* size vector including stdims */
  long *coords;     /* current coordinates */

  long k;           /* index */
  long item;        /* index thread iterator */
  long pdim = $PARENT_P(ndims);
  long rdim = $COMP(rdim);
  long prdim = (rdim < pdim) ? rdim : pdim;
  long stdim = pdim - prdim;
  	
  /* Allocate iterator and larger size vector -- do it all in one foop 
   * to avoid extra calls to smalloc.
   */
    if(!(iter = (long *)(PDL->smalloc(sizeof(long) * ($PARENT_P(ndims) * 2 + rdim))))) {
    barf("couldn't get memory for range iterator");
  }
  sizes  = iter + $PARENT_P(ndims);
  coords = sizes + $PARENT_P(ndims);

  /* Figure out size vector */
  for(ip = $COMP(sizes), sp = sizes, k=0; k<rdim; k++)
     *(sp++) = *(ip++);
  for(; k < $PARENT_P(ndims); k++) 
     *(sp++) = $PARENT_P(dims[k]);


  /* Loop over all the ranges in the index list */
  for(item=0; item<$COMP(nitems); item++) {

    /* initialize in-range iterator to loop within each range */
    for(ip = iter, k=0; k<$PARENT_P(ndims); k++) 
      *(ip++) = 0;
    
    do {
      long poff = 0;
      long coff;
      long k2;
      char trunc = 0;       /* Flag used to skip truncation case */

      /* Collect and boundary-check the current N-D coords */
      for(k=0; k < prdim; k++){

	long ck = iter[k] + $COMP(corners[ item * rdim + k  ]) ;

     	/* normal case */
    	  if(ck < 0 || ck >= $PARENT_P(dims[k])) {
 	    switch($COMP(boundary[k])) {
 	    case 0: /* no boundary breakage allowed */
 	      barf("index out-of-bounds in range");
 	      break;
 	    case 1: /* truncation */
	      trunc = 1;
 	      break;
 	    case 2: /* extension -- crop */
 	      ck = (ck >= $PARENT_P(dims[k])) ? $PARENT_P(dims[k])-1 : 0;
 	      break;
 	    case 3: /* periodic -- mod it */
 	      ck %= $PARENT_P(dims[k]);
 	      if(ck < 0)   /* Fix mod breakage in C */
 	        ck += $PARENT_P(dims[k]);
 	      break;
	    case 4: /* mirror -- reflect off the edges */
	      ck += $PARENT_P(dims[k]);
	      ck %= ($PARENT_P(dims[k]) * 2);
	      if(ck < 0) /* Fix mod breakage in C */
		ck += $PARENT_P(dims[k])*2;
	      ck -= $PARENT_P(dims[k]);
	      if(ck < 0) {
		ck *= -1;
	        ck -= 1;
	      }
	      break;
 	    default:
 	      barf("Unknown boundary condition in range -- bug alert!");
 	      break;
 	    }
 	  }

	coords[k] = ck;

      }

      /* Check extra dimensions -- pick up where k left off... */
      for( ; k < rdim ; k++) {
        /* Check for indexing off the end of the dimension list */

	long ck = iter[k] + $COMP(corners[ item * rdim + k  ]) ;

  	switch($COMP(boundary[k])) {
	    case 0: /* No boundary breakage allowed -- nonzero corners cause barfage */
	      if(ck != 0)
		 barf("too many dims in range index");
              break;
            case 1: /* truncation - just truncate if the corner is nonzero */
	      trunc |= (ck != 0);
	      break;
	    case 2: /* extension -- ignore the corner (same as 3) */
	    case 3: /* periodic  -- ignore the corner */
              ck = 0;
              break;
	    default:
	      barf("Unknown boudnary condition in range -- bug alert!");
	      /* Note clever misspelling of boundary to distinguish from other case */
	      break;
	  }
      }

      /* Find offsets into the child and parent arrays, from the N-D coords */
      /* Note we only loop over real source dims (prdim) to accumulate -- */
      /* because the offset is trivial and/or we're truncating for virtual */
      /* dims caused by permissive ranging. */
      coff = $CHILD_P(dimincs[0]) * item;
      for(k2 = $COMP(itdim), poff = k = 0; 
	  k < prdim;
	  k++) { 
	poff += coords[k]*$PARENT_P(dimincs[k]);
	if($COMP(sizes[k])) 
	  coff += iter[k] * $CHILD_P(dimincs[k2++]);
      }

      /* Loop the copy over all the source thread dims (above rdim). */
      do {
	long poff1 = poff;
	long coff1 = coff;
	
	/* Accumulate the offset due to source threading */
	for(k2 = $COMP(itdim) + $COMP(ntsize), k = rdim; 
	    k < pdim;
	    k++) {
	  poff1 += iter[k] * $PARENT_P(dimincs[k]);
	  coff1 += iter[k] * $CHILD_P(dimincs[k2++]);
	}

	/* Finally -- make the copy 
	 * EQUIVCPTRUNC works like EQUIVCPOFFS but with checking for 
	 * out-of-bounds conditions.  
	 */
	$EQUIVCPTRUNC(coff1,poff1,trunc);		
	
	/* Increment the source thread iterator */
	for( k=$COMP(rdim); 
	     k < $PARENT_P(ndims) && (++(iter[k]) >= $PARENT_P(dims[k])); 
	     k++)
	  iter[k] = 0;
      } while(k < $PARENT_P(ndims)); /* end of source-thread iteration */

      /* Increment the in-range iterator */
      for(k = 0;
	  k < $COMP(rdim) && (++(iter[k]) >= $COMP(sizes[k]));
	  k++)
	iter[k] = 0;
    } while(k < $COMP(rdim)); /* end of main iteration */
  } /* end of item do loop */

}

EOD-EquivCPOffsCode

);


pp_def(
	'rld',
	Pars=>'int a(n); b(n); [o]c(m);',
	PMCode =><<'EOD',
sub PDL::rld {
  my ($a,$b) = @_;
  my ($c);
  if ($#_ == 2) {
    $c = $_[2];
  } else {
# XXX Need to improve emulation of threading in auto-generating c
    my ($size) = $a->sumover->max;
    my (@dims) = $a->dims;
    shift @dims;
    $c = $b->zeroes($size,@dims);
  }
  &PDL::_rld_int($a,$b,$c);
  $c;
}
EOD
	Code=>'
	  int i,j=0,an;
	  $GENERIC(b) bv;
	  loop (n) %{
	    an = $a();
	    bv = $b();
	    for (i=0;i<an;i++) {
	      $c(m=>j) = bv;
	      j++;
	    }
	  %}',
	Doc => <<'EOD'
=for ref

Run-length decode a vector

Given a vector C<$a> of the numbers of instances of values C<$b>, run-length
decode to C<$c>.

=for example

 rld($a,$b,$c=null);

=cut
EOD
);
pp_def(
	'rle',
	Pars=>'c(n); int [o]a(n); [o]b(n);',
	Code=>'
	  int j=0,sn=$SIZE(n);
	  $GENERIC(c) cv, clv;
	  clv = $c(n=>0);
	  $b(n=>0) = clv;
	  $a(n=>0) = 0;
	  loop (n) %{
	    cv = $c();
	    if (cv == clv) {
	      $a(n=>j)++;
	    } else {
	      j++;
	      $b(n=>j) = clv = cv;
	      $a(n=>j) = 1;
	    }
	  %}
	  for (j++;j<sn;j++) {
	    $a(n=>j) = 0;
	    $b(n=>j) = 0;
	  }
	',
	Doc => <<'EOD'
=for ref

Run-length encode a vector

Given vector C<$c>, generate a vector C<$a> with the number of each element,
and a vector C<$b> of the unique values.  Only the elements up to the
first instance of C<0> in C<$a> should be considered.

=for example

 rle($c,$a=null,$b=null);

=cut
EOD
);

# this one can convert vaffine piddles without(!) physicalising them
# maybe it can replace 'converttypei' in the future?
#
# XXX do not know whether the HandleBad stuff will work here
#
pp_def('flowconvert',
       HandleBad => 1,
       DefaultFlow => 1,
       Reversible => 1,	
       Pars => 'PARENT(); [oca]CHILD()',
       OtherPars => 'int totype;',
       Reversible => 1,
       # Forced types
       FTypes => {CHILD => '$COMP(totype)'},
       Code => 
       '$CHILD() = $PARENT();',
       BadCode => 
       'if ( $ISBAD(PARENT()) ) { 
           $SETBAD(CHILD());
        } else {
           $CHILD() = $PARENT();
        }',
       BackCode => '$PARENT() = $CHILD();',
       BadBackCode => 
       'if ( $ISBAD(CHILD()) ) { 
           $SETBAD(PARENT());
        } else {
           $PARENT() = $CHILD();
        }',
       Doc => 'internal',
);


pp_def(
	'converttypei',
        HandleBad => 1,
	DefaultFlow => 1,
	GlobalNew => 'converttypei_new',
	OtherPars => 'int totype;',
	P2Child => 1,
	Identity => 1,
	Reversible => 1,
# Forced types
	FTypes => {CHILD => '$COMP(totype)'},
	Doc => 'internal',
);



# XXX Make clump work with optional parameter!
if(0) {
# Special-case
pp_def(
	'clump',
	DefaultFlow => 1,
	OtherPars => 'int n',
	P2Child => 1,
	Priv => 'int nnew; int nrem;',
	RedoDims => 'int i; int d1;
		if($COMP(n) > $PARENT(ndims)) 
			/* Now with more flavor: truncate overly long clumps to 
			   just clump existing dimensions...  (CED 17-Mar-2002) */

			$COMP(n) = $PARENT(ndims);

			/* Old croaking code: */
			/*$CROAK("Too many dimensions %d to clump from %d", */
			/*		$COMP(n),$PARENT(ndims)); */

		 $COMP(nrem) = ($COMP(n)==-1 ? $PARENT(threadids[0]) : $COMP(n));
		 $PRIV(nnew) = $PARENT(ndims) - $COMP(nrem) + 1;
		 $SETNDIMS($PRIV(nnew));
		 d1=1;
		 for(i=0; i<$PRIV(nrem); i++) {
		 	d1 *= $PARENT(dims[i]);
		 }
		 $CHILD(dims[0]) = d1;
		 for(; i<$PARENT(ndims); i++) {
		 	$CHILD(dims[i-$PRIV(nrem)+1]) = $PARENT(dims[i]);
		 }
		 $SETDIMS();
		 $SETDELTATHREADIDS(1-$COMP(nrem));
		 ',
	EquivCPOffsCode => '
		int i;
		for(i=0; i<$CHILD_P(nvals); i++) {
			$EQUIVCPOFFS(i,i);
		}
		',
	Reversible => 1,
);
} else {

# Affine! Make sure vaffine chaining understands to stop in the right
# place.
# the perl wrapper clump is now defined in Core.pm
# this is just the low level interface
pp_def(
	'_clump_int',
	P2Child => 1,
	DefaultFlow => 1,
	Reversible => 1,
	AffinePriv => 1,
	OtherPars => 'int n',
	RedoDims => 'int i; int d1;
		int nrem; int nnew;
		if($COMP(n) > $PARENT(ndims)) {
			/* Now with more flavor:  truncate clumping in this case to 
			 * the total number of dimensions that actually exist...
			 *  --CED 17-Mar-2002
			 */
 			$COMP(n) = -1;
			
#ifdef older_croaking_code
			$SETNDIMS(0);  /* fix to make sure we do not get problems later */
			$PRIV(offs) = 0;
			$SETDIMS();
			$CROAK("Too many dimensions %d to clump from %d",
				$COMP(n),$PARENT(ndims));
#endif
		}
		 nrem = ($COMP(n)< 0 ? $PARENT(threadids[0]) + 1 + ($COMP(n)) : $COMP(n));
	         if(nrem < 0) {
			$CROAK("Too many dimensions %d to leave behind when clumping from %d",-$COMP(n),$PARENT(ndims));
		 }

		 nnew = $PARENT(ndims) - nrem + 1;
		 $SETNDIMS(nnew);
		 $DOPRIVDIMS();
		 $PRIV(offs) = 0;
		 d1=1;
		 for(i=0; i<nrem; i++) {
		 	d1 *= $PARENT(dims[i]);
		 }
		 $CHILD(dims[0]) = d1;
		 $PRIV(incs[0]) = 1;
		 for(; i<$PARENT(ndims); i++) {
		 	$CHILD(dims[i-nrem+1]) = $PARENT(dims[i]);
			$PRIV(incs[i-nrem+1]) = $PARENT(dimincs[i]);
		 }
		 $SETDIMS();
		 $SETDELTATHREADIDS(1-nrem);
		 ',
	Doc => 'internal',
);
}

pp_def(
	'xchg',
	OtherPars => 'int n1; int n2;',
	DefaultFlow => 1,
	Reversible => 1,
	P2Child => 1,
	XCHGOnly => 1,
	EquivDimCheck => 'if ($COMP(n1) <0)
				$COMP(n1) += $PARENT(threadids[0]);
			  if ($COMP(n2) <0)
				$COMP(n2) += $PARENT(threadids[0]);
			  if ($COMP(n1) <0 ||$COMP(n2) <0 ||
			     $COMP(n1) >= $PARENT(threadids[0]) ||
			     $COMP(n2) >= $PARENT(threadids[0]))
		barf("One of dims %d, %d out of range: should be 0<=dim<%d",
			$COMP(n1),$COMP(n2),$PARENT(threadids[0]));',
	EquivPDimExpr => '(($CDIM == $COMP(n1)) ? $COMP(n2) : ($CDIM == $COMP(n2)) ? $COMP(n1) : $CDIM)',
	EquivCDimExpr => '(($PDIM == $COMP(n1)) ? $COMP(n2) : ($PDIM == $COMP(n2)) ? $COMP(n1) : $PDIM)',
	Doc => <<'EOD',
=for ref

exchange two dimensions

Negative dimension indices count from the end.

The command

=for example

 $b = $a->xchg(2,3);

creates C<$b> to be like C<$a> except that the dimensions 2 and 3
are exchanged with each other i.e.

 $b->at(5,3,2,8) == $a->at(5,3,8,2)

=cut
EOD
);

pp_addpm(<< 'EOD');

=head2 reorder

=for ref

Re-orders the dimensions of a PDL based on the supplied list.

Similar to the L<xchg|/xchg> method, this method re-orders the dimensions
of a PDL. While the L<xchg|/xchg> method swaps the position of two dimensions,
the reorder method can change the positions of many dimensions at
once.

=for usage

 # Completely reverse the dimension order of a 6-Dim array.
 $reOrderedPDL = $pdl->reorder(5,4,3,2,1,0); 

The argument to reorder is an array representing where the current dimensions
should go in the new array. In the above usage, the argument to reorder 
C<(5,4,3,2,1,0)>
indicates that the old dimensions (C<$pdl>'s dims) should be re-arranged to make the
new pdl (C<$reOrderPDL>) according to the following:

   Old Position   New Position
   ------------   ------------
   5              0
   4              1
   3 		  2
   2		  3
   1		  4
   0		  5

=for example 

Example:

 perldl> $a = sequence(5,3,2);	  # Create a 3-d Array
 perldl> p $a
 [
  [
   [ 0  1  2  3  4]
   [ 5  6  7  8  9]
   [10 11 12 13 14]
  ]
  [
   [15 16 17 18 19]
   [20 21 22 23 24]
   [25 26 27 28 29]
  ]
 ]
 perldl> p $a->reorder(2,1,0); # Reverse the order of the 3-D PDL
 [
  [
   [ 0 15]
   [ 5 20]
   [10 25]
  ]
  [
   [ 1 16]
   [ 6 21]
   [11 26]
  ]
  [
   [ 2 17]
   [ 7 22]
   [12 27]
  ]
  [
   [ 3 18]
   [ 8 23]
   [13 28]
  ]
  [
   [ 4 19]
   [ 9 24]
   [14 29]
  ]
 ]

The above is a simple example that could be duplicated by calling
C<$a-E<gt>xchg(0,2)>, but it demonstrates the basic functionality of reorder.

As this is an index function, any modifications to the
result PDL will change the parent.

=cut

sub PDL::reorder {
	my ($pdl,@newDimOrder) = @_;
	
	my $arrayMax = $#newDimOrder;
  
	#Error Checking:
	if( $pdl->getndims != scalar(@newDimOrder) ){
		my $errString = "PDL::reoderDims: Number of elements (".scalar(@newDimOrder).") in newDimOrder array doesn't\n";
		$errString .= "match the number of dims in the supplied PDL ($ndims)";
		barf($errString);
	}
	# a quicker way to do the reorder
	return $pdl->thread(@newDimOrder)->unthread(0);
}

EOD

pp_def(
	'mv',
	OtherPars => 'int n1; int n2;',
	DefaultFlow => 1,
	Reversible => 1,
	P2Child => 1,
	XCHGOnly => 1,
	EquivDimCheck => 'if ($COMP(n1) <0)
				$COMP(n1) += $PARENT(threadids[0]);
			  if ($COMP(n2) <0)
				$COMP(n2) += $PARENT(threadids[0]);
			  if ($COMP(n1) <0 ||$COMP(n2) <0 ||
			     $COMP(n1) >= $PARENT(threadids[0]) ||
			     $COMP(n2) >= $PARENT(threadids[0]))
		barf("One of dims %d, %d out of range: should be 0<=dim<%d",
			$COMP(n1),$COMP(n2),$PARENT(threadids[0]));',
	EquivPDimExpr => '(($COMP(n1) < $COMP(n2)) ?
	(($CDIM < $COMP(n1) || $CDIM > $COMP(n2)) ?
		$CDIM : (($CDIM == $COMP(n2)) ? $COMP(n1) : $CDIM+1))
	: (($COMP(n2) < $COMP(n1)) ?
		(($CDIM > $COMP(n1) || $CDIM < $COMP(n2)) ?
			$CDIM : (($CDIM == $COMP(n2)) ? $COMP(n1) : $CDIM-1))
		: $CDIM))',
	EquivCDimExpr => '(($COMP(n2) < $COMP(n1)) ?
	(($PDIM < $COMP(n2) || $PDIM > $COMP(n1)) ?
		$PDIM : (($PDIM == $COMP(n1)) ? $COMP(n2) : $PDIM+1))
	: (($COMP(n1) < $COMP(n2)) ?
		(($PDIM > $COMP(n2) || $PDIM < $COMP(n1)) ?
			$PDIM : (($PDIM == $COMP(n1)) ? $COMP(n2) : $PDIM-1))
		: $PDIM))',
	Doc => << 'EOD',
=for ref

move a dimension to another position

The command

=for example

 $b = $a->mv(4,1);

creates C<$b> to be like C<$a> except that the dimension 4 is moved to the
place 1, so:

 $b->at(1,2,3,4,5,6) == $a->at(1,5,2,3,4,6);

The other dimensions are moved accordingly.
Negative dimension indices count from the end.
EOD
);

pp_def(
	'oneslice',
	Doc => <<'EOD',
=for ref

experimental function - not for public use

=for example

 $a = oneslice();

This is not for public use currently. See the source if you have to.
This function can be used to accomplish run-time changing of
transformations i.e. changing the size of some piddle at run-time.

However, the mechanism is not yet finalized and this is just a demonstration.

EOD
	P2Child => 1,
	DefaultFlow => 1,
	Reversible => 1,
	OtherPars => 'int nth; int from; int step; int nsteps;',
	AffinePriv => 1,
	RedoDims => '
		int nth = $PRIV(nth);
		int from = $PRIV(from);
		int step = $PRIV(step);
		int nsteps = $PRIV(nsteps);
		int i;
		printf("ONESLICE_REDODIMS %d %d %d %d\n",nth,from,step,nsteps);
		if(nth >= $PARENT(ndims)) {
			die("Oneslice: too large nthdim");
		}
		if(from + step * (nsteps-1) >= $PARENT(dims[nth])) {
			die("Oneslice: too many, too large steps");
		}
		if(from < 0 || step < 0) {
			die("Oneslice: can only support positive from & step");
		}
		$PRIV(offs) = 0;
		$SETNDIMS($PARENT(ndims));
		$DOPRIVDIMS();
		for(i=0; i<$PARENT(ndims); i++) {
			$CHILD(dims)[i] = $PARENT(dims)[i];
			$PRIV(incs)[i] = $PARENT(dimincs)[i];
		}
		$CHILD(dims)[nth] = nsteps;
		$PRIV(incs)[nth] *= step;
		$PRIV(offs) += from * $PARENT(dimincs)[nth];
		$SETDELTATHREADIDS(0);
		$SETDIMS();
	',
	FooCode => # This is why we have this stupid function
	'	$COMP(from) = i1;
		$COMP(step) = i2;
		$COMP(nsteps) = i3;
		printf("ONESLICE_FOOFUNC %d %d %d %d\n",
		    $COMP(nth),$COMP(from),$COMP(step),$COMP(nsteps));
	',
);


pp_addhdr << 'EOH';
#define sign(x) ( (x) < 0 ? -1 : 1)
EOH

# I think the quotes in the =item ":" lines
# confuse the perldoc stuff
#
pp_def(
	'slice',
	Doc => << 'EOD',
=for ref

Extract a rectangular slice of a piddle, from a string specifier.

C<slice> was the original Swiss-army-knife PDL indexing routine, but is
largely superseded by the L<NiceSlice|PDL::NiceSlice> source prefilter
and its associated L<nslice|PDL::Core/nslice> method.  It is still used as the
basic underlying slicing engine for L<nslice|PDL::Core/nslice>,
and is especially useful in particular niche applications.

=for example

 $a->slice('1:3');  #  return the second to fourth elements of $a
 $a->slice('3:1');  #  reverse the above
 $a->slice('-2:1'); #  return last-but-one to second elements of $a

The argument string is a comma-separated list of what to do
for each dimension. The current formats include
the following, where I<a>, I<b> and I<c> are integers and can
take legal array index values (including -1 etc):

=over 8

=item :

takes the whole dimension intact.

=item ''

(nothing) is a synonym for ":"
(This means that C<$a-E<gt>slice(':,3')> is equal to C<$a-E<gt>slice(',3')>).

=item a

slices only this value out of the corresponding dimension.

=item (a)

means the same as "a" by itself except that the resulting
dimension of length one is deleted (so if C<$a> has dims C<(3,4,5)> then
C<$a-E<gt>slice(':,(2),:')> has dimensions C<(3,5)> whereas
C<$a-E<gt>slice(':,2,:')> has dimensions C<(3,1,5))>.

=item a:b

slices the range I<a> to I<b> inclusive out of the dimension.

=item a:b:c

slices the range I<a> to I<b>, with step I<c> (i.e. C<3:7:2> gives the indices
C<(3,5,7)>). This may be confusing to Matlab users but several other
packages already use this syntax.

=item '*'

inserts an extra dimension of width 1 and

=item '*a'

inserts an extra (dummy) dimension of width I<a>.

=back

An extension is planned for a later stage allowing
C<$a-E<gt>slice('(=1),(=1|5:8),3:6(=1),4:6')>
to express a multidimensional diagonal of C<$a>.

Trivial out-of-bounds slicing is allowed: if you slice a source
dimension that doesn't exist, but only index the 0th element, then
C<slice> treats the source as if there were a dummy dimension there.
The following are all equivalent:

	xvals(5)->dummy(1,1)->slice('(2),0')  # Add dummy dim, then slice
	xvals(5)->slice('(2),0')              # Out-of-bounds slice adds dim.
	xvals(5)->slice((2),0)                # NiceSlice syntax
	xvals(5)->((2))->dummy(0,1) 	      # NiceSlice syntax

This is an error:

	xvals(5)->slice('(2),1')	# nontrivial out-of-bounds slice dies

Because slicing doesn't directly manipulate the source and destination 
pdl -- it just sets up a transformation between them -- indexing errors 
often aren't reported until later.  This is either a bug or a feature,
depending on whether you prefer error-reporting clarity or speed of execution.

=cut
EOD
	P2Child => 1,
	DefaultFlow => 1,
	OtherPars => 'char* str',
	Comp => 'int nnew; int nthintact; int intactnew; int ndum;
	         int corresp[$COMP(intactnew)]; int start[$COMP(intactnew)];
		 int inc[$COMP(intactnew)]; int end[$COMP(intactnew)];
		 int nolddims;
		 int whichold[$COMP(nolddims)]; int oldind[$COMP(nolddims)];
		 ',
	AffinePriv => 1,
	MakeComp => q~
		int i;
		int nthnew; int nthold; int nthreal;
		int dumsize;
		char *s; char *ns;
		int nums[3]; int nthnum;
		$COMP(nnew)=0;
		$COMP(ndum)=0;
		$COMP(nolddims) = 0;
		if(str[0] == '(')
			$COMP(nolddims)++;
		else if (str[0] == '*')
			$COMP(ndum)++;
		else if (str[0] != '\0') /* handle empty string */
			$COMP(nnew)++;
		for(i=0; str[i]; i++)
			if(str[i] == ',') {
				if(str[i+1] == '(')
					$COMP(nolddims)++;
				else if(str[i+1] == '*')
					$COMP(ndum)++;
				else
					$COMP(nnew)++;
			}
		$COMP(nthintact) = $COMP(nolddims) + $COMP(nnew);
		$COMP(intactnew) = $COMP(nnew)+$COMP(ndum);
		$DOCOMPDIMS();
		nthnew=0; nthold=0; i=0; nthreal=0;
		s=str-1;
		do {
			s++;
			if(isdigit(*s) || *s == '-') {
				nthnew++; nthreal++;
				$COMP(inc[nthnew-1]) = 1;
				$COMP(corresp[nthnew-1]) = nthreal-1;
				$COMP(start[nthnew-1]) = strtol(s,&s,10);
				if(*s != ':') {
					$COMP(end[nthnew-1]) =
						$COMP(start[nthnew-1]);
					goto outlab;
				}
				s++;
				if(!isdigit(*s) && !(*s == '-')) {
					barf("Invalid slice str ind1 '%s': '%s'",str,s);
				}
				$COMP(end[nthnew-1]) = strtol(s,&s,10);
				if(*s != ':') {goto outlab;}
				s++;
				if(!isdigit(*s) && !(*s == '-')) {
					barf("Invalid slice str ind2 '%s': '%s'",str,s);
				}
				$COMP(inc[nthnew-1]) = strtol(s,&s,10);
			} else switch(*s) {
			case ':':
				s++;
				/* FALLTHRU */
			case ',': case '\0':  /* In these cases, no inc s */
				if ($COMP(intactnew) > 0) {
				  $COMP(start[nthnew]) = 0;
				  $COMP(end[nthnew]) = -1;
				  $COMP(inc[nthnew]) = 1;
				  $COMP(corresp[nthnew]) = nthreal;
				  nthnew++; nthreal++;
				}
				break;
			case '(':
				s++;
				$COMP(oldind[nthold]) = strtol(s,&s,10);
				$COMP(whichold[nthold]) = nthreal;
				nthold++; nthreal++;
				if(*s != ')') {
					barf("Sliceoblit must end with ')': '%s': '%s'",str,s);
				}
				s++;
				break;
			case '*':
				s++;
				if(isdigit(*s)) {
					dumsize = strtol(s,&s,10);
				} else {dumsize = 1;}
				$COMP(corresp[nthnew]) = -1;
				$COMP(start[nthnew]) = 0;
				$COMP(end[nthnew]) = dumsize-1;
				$COMP(inc[nthnew]) = 1;
				nthnew++;
				break;
			}
		   outlab:
			if(*s != ',' && *s != '\0') {
				barf("Invalid slice str '%s': '%s'",str,s);
			}
		} while(*s);
		$SETREVERSIBLE(1); /* XXX Only if incs>0, no dummies */
	~,
	RedoDims => '
		int i; int start; int end; int inc;
		if ($COMP(nthintact) > $PARENT(ndims)) {

	/* Slice has more dims than parent.  Check that the extra dims are 		
         * all zero, and if they are then give back What You Probably Wanted,
	 * which is a slice with dummy dimensions of order 1 in place of each excessive 
	 * dimension.  (Note that there are two ways to indicate a zero index: "0" and "-<w>",
	 * where <w> happens to be the order of that dimension in the original
	 * piddle.  The latter case still causes an error.  That is a feature.)
	 *    --CED 15-March-2002
	 */
			int ii,parentdim,ok;
			int n_xtra_dims=0, n_xtra_olddims=0;

			   /* Check index for each extra dim in the ordinary affine list */

			for(ok=1, ii = 0; ok && ii < $COMP(intactnew) ; ii++) {
				parentdim = $COMP(corresp[ii]);
/*				fprintf(stderr,"ii=%d,parent=%d, ndum=%d, nnew=%d...",ii,parentdim,$COMP(ndum),$COMP(nnew));				*/
				if(parentdim >= $PARENT(ndims)) {					

					ok = ( ( $COMP(start[ii]) == 0 ) && 
					  	( $COMP(end[ii]) == 0 || $COMP(end[ii])== -1 )
					);
					if(ok) {
						/* Change this into a dummy dimension, rank 1 */
						$COMP(corresp[ii]) = -1;
						$COMP(start[ii])   = 0;
						$COMP(end[ii])     = 0;
						$COMP(inc[ii])     = 1;
						$COMP(ndum)++;      /* One more dummy dimension... */
						$COMP(nnew)--;      /* ... one less real dimension */
						$COMP(nthintact)--; /* ... one less intact dim */
/*						fprintf(stderr,"ok, ndum=%d, nnew=%d\n",$COMP(ndum), $COMP(nnew));*/
					}
/*				fflush(stderr);*/
				}
			}	
			
			  /* Check index for each indexed parent dimension */
			for(ii=0; ok && ii < $COMP(nolddims); ii++) {
				if($COMP(whichold[ii]) >= $PARENT(ndims)) {
					ok = ( $COMP(whichold[ii]) < $PARENT(ndims) ) ||
						( $COMP(oldind[ii]) == 0 ) ||
						( $COMP(oldind[ii]) == -1) ;
					if(ok) {
					  int ij;
					  /* crunch indexed dimensions -- slow but sure */
					  $COMP(nolddims)--;
					  for(ij=ii; ij<$COMP(nolddims); ij++) {
						$COMP(oldind[ij]) = $COMP(oldind[ij+1]);
					 	$COMP(whichold[ij]) = $COMP(whichold[ij+1]);
					  }
					  $COMP(nthintact)--;
					}
				}
			}	
/*	fprintf(stderr,"ok=%d\n",ok);fflush(stderr);*/
			if(ok) {
			   /* Valid slice: all extra dims are zero. Adjust indices accordingly. */
/*			  $COMP(intactnew) -= $COMP(nthintact) - $PARENT(ndims); */
/*			  $COMP(nthintact) = $PARENT(ndims);*/
			} else { 	

			   /* Invalid slice: nonzero extra dimension.  Clean up and die.  */

  			 $SETNDIMS(0); /* dirty fix */
			 $PRIV(offs) = 0;	
			 $SETDIMS();
			 $CROAK("Too many dims in slice");
			}
 	        }
		$SETNDIMS($PARENT(ndims)-$COMP(nthintact)+$COMP(intactnew));
		$DOPRIVDIMS();
		$PRIV(offs) = 0;
		for(i=0; i<$COMP(intactnew); i++) {
			int parentdim = $COMP(corresp[i]);
			start = $COMP(start[i]); end = $COMP(end[i]);
			inc = $COMP(inc[i]);
			if(parentdim!=-1) {
				if(-start > $PARENT(dims[parentdim]) ||
				   -end > $PARENT(dims[parentdim])) {
					barf("Negative slice cannot start or end above limit");
                                }
				if(start < 0)
					start = $PARENT(dims[parentdim]) + start;
				if(end < 0)
					end = $PARENT(dims[parentdim]) + end;
				if(start >= $PARENT(dims[parentdim]) ||
				   end >= $PARENT(dims[parentdim])) {
					barf("Slice cannot start or end above limit");
				}
				if(sign(end-start)*sign(inc) < 0)
					inc = -inc;
				$PRIV(incs[i]) = $PARENT(dimincs[parentdim]) * inc;
				$PRIV(offs) += start * $PARENT(dimincs[parentdim]);
			} else {
				$PRIV(incs[i]) = 0;
			}
			$CHILD(dims[i]) = ((int)((end-start)/inc))+1;
                        if ($CHILD(dims[i]) <= 0)
                           barf("slice internal error: computed slice dimension must be positive");
		}
		for(i=$COMP(nthintact); i<$PARENT(ndims); i++) {
			int cdim = i - $COMP(nthintact) + $COMP(intactnew);
			$PRIV(incs[cdim]) = $PARENT(dimincs[i]);
			$CHILD(dims[cdim]) = $PARENT(dims[i]);
		}
		for(i=0; i<$COMP(nolddims); i++) {
			int oi = $COMP(oldind[i]);
			int wo = $COMP(whichold[i]);
			if(oi < 0)
				oi += $PARENT(dims[wo]);
			if( oi >= $PARENT(dims[wo]) )
				$CROAK("Cannot obliterate dimension after end");
			$PRIV(offs) += $PARENT(dimincs[wo])
					* oi;
		}
	/*
		for(i=0; i<$CHILD(ndims)-$PRIV(nnew); i++) {
			$CHILD(dims[i+$COMP(intactnew)]) =
				$PARENT(dims[i+$COMP(nthintact)]);
			$PRIV(incs[i+$COMP(intactnew)]) =
				$PARENT(dimincs[i+$COMP(nthintact)]);
		}
	*/
		$SETDIMS();
	',
);

pp_addpm(<<'EOD'

=head2 using

=for ref

Returns array of column numbers requested

=for usage

 line $pdl->using(1,2);

Plot, as a line, column 1 of C<$pdl> vs. column 2

=for example

 perldl> $pdl = rcols("file");
 perldl> line $pdl->using(1,2);

=cut

*using = \&PDL::using;
sub PDL::using {
  my ($x,@ind)=@_;
  @ind = list $ind[0] if (ref $ind[0] eq 'PDL');
  foreach (@ind) {
    $_ = $x->slice("($_)");
  }
  @ind;
}

EOD
);

pp_add_exported('', 'using');

pp_addhdr(<<END
static int cmp_pdll(const void *a_,const void *b_) {
	PDL_Long *a = (PDL_Long *)a_; PDL_Long *b=(PDL_Long *)b_;
	if(*a>*b) return 1;
	else if(*a==*b) return 0;
	else return -1;
}
END
);
	

pp_def( 'affine',
	P2Child => 1,
	DefaultFlow => 1,
	Reversible => 1,
	AffinePriv => 1,
	GlobalNew => 'affine_new',
	OtherPars => 'int offspar; SV *dimlist; SV *inclist;',
	Comp => 'int nd; PDL_Long offset; PDL_Long sdims[$COMP(nd)];
		PDL_Long sincs[$COMP(nd)];',
	MakeComp => '
		int i,n2;
		PDL_Long *tmpi;
		PDL_Long *tmpd = PDL->packdims(dimlist,&($COMP(nd)));
		tmpi = PDL->packdims(inclist,&n2);		
		if ($COMP(nd) < 0) {
		      $CROAK("Affine: can not have negative no of dims");
		}
		if ($COMP(nd) != n2)
		      $CROAK("Affine: number of incs does not match dims");
		$DOCOMPDIMS();
		$COMP(offset) = offspar;
		for (i=0; i<$COMP(nd); i++) {
			$COMP(sdims)[i] = tmpd[i];
			$COMP(sincs)[i] = tmpi[i];
		}
		',
	RedoDims => '
		int i;
		$SETNDIMS($COMP(nd));
		$DOPRIVDIMS();
		$PRIV(offs) = $COMP(offset);
		for (i=0;i<$CHILD(ndims);i++) {
			$PRIV(incs)[i] = $COMP(sincs)[i];
			$CHILD(dims)[i] = $COMP(sdims)[i];
		}
		$SETDIMS();
		',
	Doc => undef,
);

pp_def(
	'diagonalI',
	P2Child => 1,
	DefaultFlow => 1,
	Reversible => 1,
	AffinePriv => 1,
	OtherPars => 'SV *list',
	Comp => 'int nwhichdims; PDL_Long whichdims[$COMP(nwhichdims)];',
	MakeComp => '
		int i,j;
		PDL_Long *tmp= PDL->packdims(list,&($COMP(nwhichdims)));
		if($COMP(nwhichdims) < 1) {
			$CROAK("Diagonal: must have at least 1 dimension");
		}
		$DOCOMPDIMS();
		for(i=0; i<$COMP(nwhichdims); i++)
			$COMP(whichdims)[i] = tmp[i];
		qsort($COMP(whichdims), $COMP(nwhichdims), sizeof(PDL_Long),
			cmp_pdll);
	',
	RedoDims => '
		int nthp,nthc,nthd; int cd = $COMP(whichdims[0]);
		$SETNDIMS($PARENT(ndims)-$COMP(nwhichdims)+1);
		$DOPRIVDIMS();
		$PRIV(offs) = 0;
		if ($COMP(whichdims)[$COMP(nwhichdims)-1] >= $PARENT(ndims) ||
			$COMP(whichdims)[0] < 0)
			$CROAK("Diagonal: dim out of range");
		nthd=0; nthc=0;
		for(nthp=0; nthp<$PARENT(ndims); nthp++)
			if (nthd < $COMP(nwhichdims) &&
			    nthp == $COMP(whichdims)[nthd]) {
				if (!nthd) {
					$CHILD(dims)[cd] = $PARENT(dims)[cd];
					nthc++;
					$PRIV(incs)[cd] = 0;
				}
				if (nthd && $COMP(whichdims)[nthd] ==
				    $COMP(whichdims)[nthd-1])
				       $CROAK("Diagonal: dims must be unique");
				nthd++; /* advance pointer into whichdims */
				if($CHILD(dims)[cd] !=
				    $PARENT(dims)[nthp]) {
					$CROAK("Different dims %d and %d",
						$CHILD(dims)[cd],
						$PARENT(dims)[nthp]);
				}
				$PRIV(incs)[cd] += $PARENT(dimincs)[nthp];
			} else {
				$PRIV(incs)[nthc] = $PARENT(dimincs)[nthp];
				$CHILD(dims)[nthc] = $PARENT(dims)[nthp];
				nthc++;
			}
		$SETDIMS();
	',
	Doc => << 'EOD',
=for ref

Returns the multidimensional diagonal over the specified dimensions.

The diagonal is placed at the first (by number) dimension that is
diagonalized.
The other diagonalized dimensions are removed. So if C<$a> has dimensions
C<(5,3,5,4,6,5)> then after

=for example

 $b = $a->diagonal(0,2,5);

the piddle C<$b> has dimensions C<(5,3,4,6)> and 
C<$b-E<gt>at(2,1,0,1)> refers
to C<$a-E<gt>at(2,1,2,0,1,2)>.

NOTE: diagonal doesn't handle threadids correctly. XXX FIX

EOD
);

pp_def(
	'lags',
	Doc => <<'EOD',
=for ref

Returns a piddle of lags to parent.

Usage:

=for usage

  $lags = $a->lags($nthdim,$step,$nlags);

I.e. if C<$a> contains

 [0,1,2,3,4,5,6,7]

then

=for example

 $b = $a->lags(0,2,2);

is a (5,2) matrix

 [2,3,4,5,6,7]
 [0,1,2,3,4,5]

This order of returned indices is kept because the function is
called "lags" i.e. the nth lag is n steps behind the original.

C<$step> and C<$nlags> must be positive. C<$nthdim> can be
negative and will then be counted from the last dim backwards
in the usual way (-1 = last dim).

EOD
	P2Child => 1,
	DefaultFlow => 1,
	Reversible => 1, # XXX Not really
	AffinePriv => 1,
	OtherPars => 'int nthdim; int step; int n;',
	RedoDims => '
		int i;
		if ($PRIV(nthdim) < 0)  /* the usual conventions */
	           $PRIV(nthdim) = $PARENT(ndims) + $PRIV(nthdim);
		if ($PRIV(nthdim) < 0 || $PRIV(nthdim) >= $PARENT(ndims))
		   $CROAK("lags: dim out of range");
		if ($COMP(n) < 1)
		   $CROAK("lags: number of lags must be positive");
		if ($COMP(step) < 1)
		   $CROAK("lags: step must be positive");
		$PRIV(offs) = 0;
		$SETNDIMS($PARENT(ndims)+1);
		$DOPRIVDIMS();
		for(i=0; i<$PRIV(nthdim); i++) {
			$CHILD(dims)[i] = $PARENT(dims)[i];
			$PRIV(incs)[i] = $PARENT(dimincs)[i];
		}
		$CHILD(dims)[i] = $PARENT(dims)[i] - $COMP(step) * ($COMP(n)-1);
		if ($CHILD(dims)[i] < 1)
		  $CROAK("lags: product of step size and "
			 "number of lags too large");
		$CHILD(dims)[i+1] = $COMP(n);
		$PRIV(incs)[i] = ($PARENT(dimincs)[i]);
		$PRIV(incs)[i+1] = - $PARENT(dimincs)[i] * $COMP(step);
                $PRIV(offs) += ($CHILD(dims)[i+1] - 1) * (-$PRIV(incs)[i+1]);
		i++;
		for(; i<$PARENT(ndims); i++) {
			$CHILD(dims)[i+1] = $PARENT(dims)[i];
			$PRIV(incs)[i+1] = $PARENT(dimincs)[i];
		}
		$SETDIMS();
	'
);

pp_def(
	'splitdim',
	Doc => <<'EOD',
=for ref

Splits a dimension in the parent piddle (opposite of L<clump|PDL::Core/clump>)

After

=for example

 $b = $a->splitdim(2,3);

the expression

 $b->at(6,4,x,y,3,6) == $a->at(6,4,x+3*y)

is always true (C<x> has to be less than 3).

EOD
	P2Child => 1,
	DefaultFlow => 1,
	Reversible => 1, # XXX Not really
	OtherPars => 'int nthdim; int nsp;',
	AffinePriv => 1,
	RedoDims => '
		int i = $COMP(nthdim);
		int nsp = $COMP(nsp);
		if(nsp == 0) {die("Splitdim: Cannot split to 0\n");}
		if(i <0 || i >= $PARENT(ndims)) {
			die("Splitdim: nthdim (%d) must not be negative or greater or equal to number of dims (%d)\n",
				i, $PARENT(ndims));
		}
		if(nsp > $PARENT(dims[i])) {
			die("Splitdim: nsp (%d) cannot be greater than dim (%d)\n",
				nsp, $PARENT(dims[i]));
		}
		$PRIV(offs) = 0;
		$SETNDIMS($PARENT(ndims)+1);
		$DOPRIVDIMS();
		for(i=0; i<$PRIV(nthdim); i++) {
			$CHILD(dims)[i] = $PARENT(dims)[i];
			$PRIV(incs)[i] = $PARENT(dimincs)[i];
		}
		$CHILD(dims)[i] = $COMP(nsp);
		$CHILD(dims)[i+1] = $PARENT(dims)[i] / $COMP(nsp);
		$PRIV(incs)[i] = $PARENT(dimincs)[i];
		$PRIV(incs)[i+1] = $PARENT(dimincs)[i] * $COMP(nsp);
		i++;
		for(; i<$PARENT(ndims); i++) {
			$CHILD(dims)[i+1] = $PARENT(dims)[i];
			$PRIV(incs)[i+1] = $PARENT(dimincs)[i];
		}
		$SETDIMS();
	',
);

pp_def('rotate',
	Doc => <<'EOD',
=for ref

Shift vector elements along with wrap. Flows data back&forth.

EOD
 	Pars=>'x(n); int shift(); [oca]y(n)',
        DefaultFlow => 1,
        Reversible => 1,
        Code=>'
        int i,j;
        int n_size = $SIZE(n);
	if (n_size == 0)
          barf("can not shift zero size piddle (n_size is zero)");
        j = ($shift()) % n_size;
        if (j < 0)
                j += n_size;
        for(i=0; i<n_size; i++,j++) {
            if (j == n_size)
               j = 0;
            $y(n=>j) = $x(n=>i);
        }',
        BackCode=>'
        int i,j;
        int n_size = $SIZE(n);
        j = ($shift()) % n_size;
        if (j < 0)
                j += n_size;
        for(i=0; i<n_size; i++,j++) {
            if (j == n_size)
               j = 0;
            $x(n=>i) = $y(n=>j);
        }
        '
);

# This is a bit tricky. Hope I haven't missed any cases.

pp_def(
	'threadI',
	Doc => <<'EOD',
=for ref

internal

Put some dimensions to a threadid.

=for example

 $b = $a->threadI(0,1,5); # thread over dims 1,5 in id 1

EOD
	P2Child => 1,
	DefaultFlow => 1,
	Reversible => 1,
	AffinePriv => 1,
	CallCopy => 0,  # Don't CallCopy for subclassed objects because PDL::Copy calls ThreadI
			#  (Wouldn't cause recursive loop otherwise)
	OtherPars => 'int id; SV *list',
	Comp => 'int id; int nwhichdims; PDL_Long whichdims[$COMP(nwhichdims)];
			int nrealwhichdims; ',
	MakeComp => '
		int i,j;
		PDL_Long *tmp= PDL->packdims(list,&($COMP(nwhichdims)));
		$DOCOMPDIMS();
		for(i=0; i<$COMP(nwhichdims); i++)
			$COMP(whichdims)[i] = tmp[i];
		$COMP(nrealwhichdims) = 0;
		for(i=0; i<$COMP(nwhichdims); i++) {
			for(j=i+1; j<$COMP(nwhichdims); j++)
				if($COMP(whichdims[i]) == $COMP(whichdims[j]) &&
				   $COMP(whichdims[i]) != -1) {
				$CROAK("Thread: duplicate arg %d %d %d",
					i,j,$COMP(whichdims[i]));
			}
			if($COMP(whichdims)[i] != -1) {
				$COMP(nrealwhichdims) ++;
			}
		}
		$COMP(id) = id;
		',
	RedoDims => '
		int nthc,i,j,flag;
		$SETNDIMS($PARENT(ndims));
		$DOPRIVDIMS();
		$PRIV(offs) = 0;
		nthc=0;
		for(i=0; i<$PARENT(ndims); i++) {
			flag=0;
			if($PARENT(nthreadids) > $COMP(id) &&
			   i == $PARENT(threadids[$COMP(id)])) {
			   nthc += $COMP(nwhichdims);
			}
			for(j=0; j<$COMP(nwhichdims); j++) {
				if($COMP(whichdims[j] == i)) {flag=1; break;}
			}
			if(flag) {
				continue;
			}
			$CHILD(dims[nthc]) = $PARENT(dims[i]);
			$PRIV(incs[nthc]) = $PARENT(dimincs[i]);
			nthc++;
		}
		for(i=0; i<$COMP(nwhichdims); i++) {
			int cdim,pdim;
			cdim = i +
			 ($PARENT(nthreadids) > $COMP(id) ?
			  $PARENT(threadids[$COMP(id)]) : $PARENT(ndims))
			  - $COMP(nrealwhichdims);
			pdim = $COMP(whichdims[i]);
			if(pdim == -1) {
				$CHILD(dims[cdim]) = 1;
				$PRIV(incs[cdim]) = 0;
			} else {
				$CHILD(dims[cdim]) = $PARENT(dims[pdim]);
				$PRIV(incs[cdim]) = $PARENT(dimincs[pdim]);
			}
		}
		$SETDIMS();
		PDL->reallocthreadids($CHILD_PTR(),
			($PARENT(nthreadids)<=$COMP(id) ?
				$COMP(id)+1 : $PARENT(nthreadids)));
		for(i=0; i<$CHILD(nthreadids); i++) {
			$CHILD(threadids[i]) =
			 ($PARENT(nthreadids) > i ?
			  $PARENT(threadids[i]) : $PARENT(ndims)) +
			 (i <= $COMP(id) ? - $COMP(nrealwhichdims) :
			  $COMP(nwhichdims) - $COMP(nrealwhichdims));
		}
		$CHILD(threadids[$CHILD(nthreadids)]) = $CHILD(ndims);
		',
);


# we don't really need this one since it can be achieved with
# a ->threadI(-1,[])
pp_def('identvaff',
	P2Child => 1,
	DefaultFlow => 1,
	Reversible => 1,
	AffinePriv => 1,
	RedoDims => '
		int i;
		$SETNDIMS($PARENT(ndims));
		$DOPRIVDIMS();
		$PRIV(offs) = 0;
		for(i=0; i<$PARENT(ndims); i++) {
			$CHILD(dims[i]) = $PARENT(dims[i]);
			$PRIV(incs[i]) = $PARENT(dimincs[i]);
		}
		$SETDIMS();
		$SETDELTATHREADIDS(0);
		$CHILD(threadids[$CHILD(nthreadids)]) = $CHILD(ndims);
		',
	Doc => <<'EOD',
=for ref

A vaffine identity transformation (includes thread_id copying).

Mainly for internal use.

EOD
);


pp_def(
	'unthread',
	Doc => <<'EOD',
=for ref

All threaded dimensions are made real again.

See [TBD Doc] for details and examples.

EOD
	P2Child => 1,
	DefaultFlow => 1,
	Reversible => 1,
	AffinePriv => 1,
	OtherPars => 'int atind;',
	RedoDims => '
		int i;
		$SETNDIMS($PARENT(ndims));
		$DOPRIVDIMS();
		$PRIV(offs) = 0;
		for(i=0; i<$PARENT(ndims); i++) {
			int corc;
			if(i<$COMP(atind)) {
				corc = i;
			} else if(i < $PARENT(threadids[0])) {
				corc = i + $PARENT(ndims)-$PARENT(threadids[0]);
			} else {
				corc = i - $PARENT(threadids[0]) + $COMP(atind);
			}
			$CHILD(dims[corc]) = $PARENT(dims[i]);
			$PRIV(incs[corc]) = $PARENT(dimincs[i]);
		}
		$SETDIMS();
	',
);


pp_add_exported('', 'dice dice_axis');
pp_addpm(<<'EOD');

=head2 dice

=for ref

Dice rows/columns/planes out of a PDL using indexes for
each dimension.

This function can be used to extract irregular subsets
along many dimension of a PDL, e.g. only certain rows in an image,
or planes in a cube. This can of course be done with
the usual dimension tricks but this saves having to
figure it out each time!

This method is similar in functionality to the L<slice|/slice>
method, but L<slice|/slice> requires that contiguous ranges or ranges
with constant offset be extracted. ( i.e. L<slice|/slice> requires 
ranges of the form C<1,2,3,4,5> or C<2,4,6,8,10>). Because of this
restriction, L<slice|/slice> is more memory efficient and slightly faster
than dice

=for usage

 $slice = $data->dice([0,2,6],[2,1,6]); # Dicing a 2-D array

The arguments to dice are arrays (or 1D PDLs) for each dimension
in the PDL. These arrays are used as indexes to which rows/columns/cubes,etc
to dice-out (or extract) from the C<$data> PDL. 

Use C<X> to select all indices along a given dimension (compare also
L<mslice|PDL::Core/mslice>). As usual (in slicing methods) trailing
dimensions can be omitted implying C<X>'es for those.

=for example 

 perldl> $a = sequence(10,4)
 perldl> p $a
 [
  [ 0  1  2  3  4  5  6  7  8  9]
  [10 11 12 13 14 15 16 17 18 19]
  [20 21 22 23 24 25 26 27 28 29]
  [30 31 32 33 34 35 36 37 38 39]
 ]
 perldl> p $a->dice([1,2],[0,3]) # Select columns 1,2 and rows 0,3
 [
  [ 1  2]
  [31 32]
 ]
 perldl> p $a->dice(X,[0,3])
 [
  [ 0  1  2  3  4  5  6  7  8  9]
  [30 31 32 33 34 35 36 37 38 39]
 ]
 perldl> p $a->dice([0,2,5])
 [
  [ 0  2  5]
  [10 12 15]
  [20 22 25]
  [30 32 35]
 ]

As this is an index function, any modifications to the
slice change the parent (use the C<.=> operator).

=cut

sub PDL::dice { 

	my $self = shift;
	my @dim_indexes = @_;  # array of dimension indexes
	
	# Check that the number of dim indexes <=
	#    number of dimensions in the PDL
	my $no_indexes = scalar(@dim_indexes);
	my $noDims = $self->getndims;
	barf("PDL::dice: Number of index arrays ($no_indexes) not equal to the dimensions of the PDL ($noDims")
			 if $no_indexes > $noDims;
	my $index;
	my $pdlIndex;
	my $outputPDL=$self;
	my $indexNo = 0;

	# Go thru each index array and dice the input PDL:
	foreach $index(@dim_indexes){
		$outputPDL = $outputPDL->dice_axis($indexNo,$index)
			unless !ref $index && $index eq 'X';

		$indexNo++;
	}

	return $outputPDL;
}  
*dice = \&PDL::dice;


=head2 dice_axis

=for ref

Dice rows/columns/planes from a single PDL axis (dimension)
using index along a specified axis

This function can be used to extract irregular subsets
along any dimension, e.g. only certain rows in an image,
or planes in a cube. This can of course be done with
the usual dimension tricks but this saves having to
figure it out each time!

=for usage

 $slice = $data->dice_axis($axis,$index);

=for example

 perldl> $a = sequence(10,4)
 perldl> $idx = pdl(1,2)
 perldl> p $a->dice_axis(0,$idx) # Select columns
 [
  [ 1  2]
  [11 12]
  [21 22]
  [31 32]
 ]
 perldl> $t = $a->dice_axis(1,$idx) # Select rows
 perldl> $t.=0
 perldl> p $a
 [
  [ 0  1  2  3  4  5  6  7  8  9]
  [ 0  0  0  0  0  0  0  0  0  0]
  [ 0  0  0  0  0  0  0  0  0  0]
  [30 31 32 33 34 35 36 37 38 39]
 ]

The trick to using this is that the index selects
elements along the dimensions specified, so if you
have a 2D image C<axis=0> will select certain C<X> values
- i.e. extract columns

As this is an index function, any modifications to the
slice change the parent.

=cut

sub PDL::dice_axis { 
  my($self,$axis,$idx) = @_;
  
  # Convert to PDLs: array refs using new, otherwise use topdl:
  my $ix = (ref($idx) eq 'ARRAY') ? ref($self)->new($idx) : ref($self)->topdl($idx);
  my $n = $self->getndims;
  my $a = $ix->getndims;
  barf("index_axis: index must be <=1D") if $a>1;
  for ($a..$n-1) {
     $ix = $ix->dummy(0);
  }
  	
  return $self->mv($axis,0)->index($ix)->mv($n-1,$axis);
}  
*dice_axis = \&PDL::dice_axis;


EOD

pp_done();
__DATA__

# A very useful transformation for e.g. axis values: hexagonal
# arrays can be made like this.

if(0) {
deftrans(
	Name => 'repeat',
	Pars => 'int whichind, int howmany',
	MakeComp => '
		$COMP(howmany) = howmany;
		$COMP(whichind) = whichind;
		$SETREVERSIBLE($COMP(howmany)==1);
		',
	Dims => '
		$SETNDIMS($PARENT(ndims));
		LOOPDIMS %{
			$CHILD(dims[$DIM]) = $PARENT(dims[$DIM]);
		%}
		$CHILD(dims[$COMP(whichind)]) *= $PRIV(howmany);
		$SETDIMS();
		',
	ParentInds =>
		'$COPYINDS();
		 $PARENTINDS($COMP(whichind)) %= $PARENT(dims[$PRIV(whichind)]);',
	Print => 'printf("REPEAT: %d, %d\n",
			$COMP(whichind), $COMP(howmany));'
);
}

# Parent's first index is value of indices.

if(0) {
deftrans(
	Name => 'indexed',
	Pars => 'pdl* indices',
	Dims => '
		$SETNDIMS($COMP(indices)->ndims);
		LOOPDIMS %{
			$CHILD(dims[$DIM]) = $COMP(indices)->dims[$DIM];
		%}
		',
	ParentInds =>
		'$COPYINDS();
		 $PARENTINDS(0) = PDL->get($COMP(indices),&($MYINDS(0)));'
);
}