"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Basic/Slices/slices.pd" between
PDL-2.082.tar.gz and PDL-2.083.tar.gz

About: PDL (Perl Data Language) aims to turn perl into an efficient numerical language for scientific computing (similar to IDL and MatLab).

slices.pd  (PDL-2.082):slices.pd  (PDL-2.083)
use PDL::Types qw(ppdefs_all); use PDL::Types qw(ppdefs_all);
use strict; use strict;
use warnings; use warnings;
pp_addpm({At => 'Top'},<< 'EOD'); pp_addpm({At => 'Top'},<<'EOD');
=head1 NAME =head1 NAME
PDL::Slices -- Indexing, slicing, and dicing PDL::Slices -- Indexing, slicing, and dicing
=head1 SYNOPSIS =head1 SYNOPSIS
use PDL; use PDL;
$x = ones(3,3); $x = ones(3,3);
$y = $x->slice('-1:0,(1)'); $y = $x->slice('-1:0,(1)');
skipping to change at line 178 skipping to change at line 178
'index barfs if any of the index values are bad.', 'index barfs if any of the index values are bad.',
); );
pp_def( pp_def(
'index1d', 'index1d',
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
HandleBad => 1, HandleBad => 1,
DefaultFlow => 1, DefaultFlow => 1,
TwoWay => 1, TwoWay => 1,
Pars => 'a(n); indx ind(m); [oca] c(m);', Pars => 'a(n); indx ind(m); [oca] c(m);',
Code => q{ Code => pp_line_numbers(__LINE__, <<'EOF'),
PDL_Indx i; PDL_Indx i;
for(i=0;i<$SIZE(m);i++) { for(i=0;i<$SIZE(m);i++) {
PDL_Indx this_ind = $ind(m=>i); PDL_Indx this_ind = $ind(m=>i);
PDL_IF_BAD(if( $ISBADVAR(this_ind, ind) ) { PDL_IF_BAD(if( $ISBADVAR(this_ind, ind) ) {
$SETBAD(c(m=>i)); $SETBAD(c(m=>i));
} else,) { } else,) {
if( this_ind<0 || this_ind >= $SIZE(n) ) { if( this_ind<0 || this_ind >= $SIZE(n) ) {
$CROAK("invalid index %"IND_FLAG" at pos %"IND_FLAG" (valid range 0..%"IND_FLAG")", $CROAK("invalid index %"IND_FLAG" at pos %"IND_FLAG" (valid range 0..%"IND_FLAG")",
this_ind, i, $SIZE(n)-1); this_ind, i, $SIZE(n)-1);
} }
$c(m=>i) = $a(n=>this_ind); $c(m=>i) = $a(n=>this_ind);
} }
} }
}, EOF
BackCode => q{ BackCode => pp_line_numbers(__LINE__, <<'EOF'),
PDL_Indx i; PDL_Indx i;
for(i=0;i<$SIZE(m);i++) { for(i=0;i<$SIZE(m);i++) {
PDL_Indx this_ind = $ind(m=>i); PDL_Indx this_ind = $ind(m=>i);
PDL_IF_BAD(if( $ISBADVAR(this_ind, ind) ) { PDL_IF_BAD(if( $ISBADVAR(this_ind, ind) ) {
/* do nothing */ /* do nothing */
} else,) { } else,) {
if( this_ind<0 || this_ind >= $SIZE(n) ) { if( this_ind<0 || this_ind >= $SIZE(n) ) {
$CROAK("invalid index %"IND_FLAG" at pos %"IND_FLAG" (valid range 0..%"IND_FLAG")", $CROAK("invalid index %"IND_FLAG" at pos %"IND_FLAG" (valid range 0..%"IND_FLAG")",
this_ind, i, $SIZE(n)-1); this_ind, i, $SIZE(n)-1);
} }
$a(n=>this_ind) = $c(m=>i); $a(n=>this_ind) = $c(m=>i);
} }
} }
}, EOF
Doc => $doc, Doc => $doc,
BadDoc => BadDoc =>
'index1d propagates BAD index elements to the output variable.' 'index1d propagates BAD index elements to the output variable.'
); );
my $index2d_init = my $index2d_init = pp_line_numbers(__LINE__, <<'EOF');
'register PDL_Indx this_ind_a = $inda(),this_ind_b = $indb(); register PDL_Indx this_ind_a = $inda(),this_ind_b = $indb();
if( PDL_IF_BAD($ISBADVAR(this_ind_a,inda) ||,) this_ind_a<0 || this_ind_ a>=$SIZE(na) ) if( PDL_IF_BAD($ISBADVAR(this_ind_a,inda) ||,) this_ind_a<0 || this_ind_ a>=$SIZE(na) )
$CROAK("invalid x-index %"IND_FLAG" (valid range 0..%"IND_FLAG")", $CROAK("invalid x-index %"IND_FLAG" (valid range 0..%"IND_FLAG")",
this_ind_a,$SIZE(na)-1); this_ind_a,$SIZE(na)-1);
if( PDL_IF_BAD($ISBADVAR(this_ind_b,indb) ||,) this_ind_b<0 || this_ind_ b>=$SIZE(nb) ) if( PDL_IF_BAD($ISBADVAR(this_ind_b,indb) ||,) this_ind_b<0 || this_ind_ b>=$SIZE(nb) )
$CROAK("invalid y-index %"IND_FLAG" (valid range 0..%"IND_FLAG")", $CROAK("invalid y-index %"IND_FLAG" (valid range 0..%"IND_FLAG")",
this_ind_b,$SIZE(nb)-1); this_ind_b,$SIZE(nb)-1);
'; EOF
pp_def( pp_def(
'index2d', 'index2d',
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
HandleBad => 1, HandleBad => 1,
DefaultFlow => 1, DefaultFlow => 1,
TwoWay => 1, TwoWay => 1,
Pars => 'a(na,nb); indx inda(); indx indb(); [oca] c();', Pars => 'a(na,nb); indx inda(); indx indb(); [oca] c();',
Code => Code =>
$index2d_init . ' $c() = $a(na => this_ind_a, nb => this_ind_b);', $index2d_init . pp_line_numbers(__LINE__-1, ' $c() = $a(na => this_ind_a, nb => this_ind_b);'),
BackCode => BackCode =>
$index2d_init . ' $a(na => this_ind_a, nb => this_ind_b) = $c();', $index2d_init . pp_line_numbers(__LINE__-1, ' $a(na => this_ind_a, nb => this_ind_b) = $c();'),
Doc => $doc, Doc => $doc,
BadDoc => BadDoc =>
'index2d barfs if either of the index values are bad.', 'index2d barfs if either of the index values are bad.',
); );
# indexND: CED 2-Aug-2002 # indexND: CED 2-Aug-2002
pp_add_exported('','indexND indexNDb'); pp_add_exported('','indexND indexNDb');
pp_addpm(<<'EOD-indexND'); pp_addpm(<<'EOD-indexND');
=head2 indexNDb =head2 indexNDb
skipping to change at line 293 skipping to change at line 293
sub PDL::indexND { sub PDL::indexND {
my($source,$index, $boundary) = @_; my($source,$index, $boundary) = @_;
return PDL::range($source,$index,undef,$boundary); return PDL::range($source,$index,undef,$boundary);
} }
*PDL::indexNDb = \&PDL::indexND; *PDL::indexNDb = \&PDL::indexND;
EOD-indexND EOD-indexND
pp_addpm(<<'EOD-range'); pp_addpm(<<'EOD-range');
sub PDL::range { sub PDL::range {
my($source,$ind,$sz,$bound) = @_; my($source,$ind,$sz,$bound) = @_;
# Convert to indx type up front (also handled in rangeb if necessary) # Convert to indx type up front (also handled in rangeb if necessary)
my $index = (ref $ind && UNIVERSAL::isa($ind,'PDL') && $ind->type eq 'indx') ? $ind : indx($ind); my $index = (ref $ind && UNIVERSAL::isa($ind,'PDL') && $ind->type eq 'indx') ? $ind : indx($ind);
my $size = defined($sz) ? PDL->pdl($sz) : undef; my $size = defined($sz) ? PDL->pdl($sz) : undef;
# Handle empty PDL case: return a properly constructed Empty. # Handle empty PDL case: return a properly constructed Empty.
if($index->isempty) { if($index->isempty) {
my @sdims= $source->dims; my @sdims= $source->dims;
skipping to change at line 334 skipping to change at line 333
} }
} }
no warnings; # shut up about passing undef into rangeb no warnings; # shut up about passing undef into rangeb
$source->rangeb($index,$size,$bound); $source->rangeb($index,$size,$bound);
} }
EOD-range EOD-range
pp_def( pp_def(
'rangeb', 'rangeb',
OtherPars => 'pdl *ind_pdl; SV *size; SV *boundary_sv', OtherPars => 'pdl *ind_pdl; SV *size_sv; SV *boundary_sv',
Doc => <<'EOD', Doc => <<'EOD',
=for ref =for ref
Engine for L</range> Engine for L</range>
=for example =for example
Same calling convention as L</range>, but you must supply all Same calling convention as L</range>, but you must supply all
parameters. C<rangeb> is marginally faster as it makes a direct PP call, parameters. C<rangeb> is marginally faster as it makes a direct PP call,
avoiding the perl argument-parsing step. avoiding the perl argument-parsing step.
skipping to change at line 539 skipping to change at line 538
[0 0 0 0 0] [0 0 0 0 0]
[0 0 1 1 0] [0 0 1 1 0]
] ]
B<CAVEAT>: It's quite possible to select multiple ranges that B<CAVEAT>: It's quite possible to select multiple ranges that
intersect. In that case, modifying the ranges doesn't have a intersect. In that case, modifying the ranges doesn't have a
guaranteed result in the original PDL -- the result is an arbitrary guaranteed result in the original PDL -- the result is an arbitrary
choice among the valid values. For some things that's OK; but for choice among the valid values. For some things that's OK; but for
others it's not. In particular, this doesn't work: others it's not. In particular, this doesn't work:
pdl> $photon_list = new PDL::RandVar->sample(500)->reshape(2,250)*10 pdl> $photon_list = PDL::RandVar->new->sample(500)->reshape(2,250)*10
pdl> $histogram = zeroes(10,10) pdl> $histogram = zeroes(10,10)
pdl> $histogram->range($photon_list,1)++; #not what you wanted pdl> $histogram->range($photon_list,1)++; #not what you wanted
The reason is that if two photons land in the same bin, then that bin 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...) doesn't get incremented twice. (That may get fixed in a later version...)
B<PERMISSIVE RANGING>: If C<$index> has too many dimensions compared B<PERMISSIVE RANGING>: If C<$index> has too many dimensions compared
to C<$source>, then $source is treated as though it had dummy to C<$source>, then $source is treated as though it had dummy
dimensions of size 1, up to the required number of dimensions. These dimensions of size 1, up to the required number of dimensions. These
virtual dummy dimensions have the usual boundary conditions applied to virtual dummy dimensions have the usual boundary conditions applied to
skipping to change at line 601 skipping to change at line 600
makes it easier to loop over the index array but a little more brain-bending makes it easier to loop over the index array but a little more brain-bending
to tease out the algorithm. to tease out the algorithm.
=cut =cut
EOD EOD
HandleBad => 1, HandleBad => 1,
TwoWay => 1, TwoWay => 1,
P2Child => 1, P2Child => 1,
# rdim: dimensionality of each range (0 dim of index PDL) # rdim: dimensionality of each range (0 dim of index PDL)
# # nitems: total number of list elements (product of itdims)
# itdim: number of index broadcast dimensions
# ntsize: number of nonzero size dimensions # ntsize: number of nonzero size dimensions
# bsize: Number of independently specified boundary conditions
# nsizes: Number of independently specified range dim sizes
# sizes: array of range sizes, indexed (0..rdim-1). A zero element means # sizes: array of range sizes, indexed (0..rdim-1). A zero element means
# that the dimension is omitted from the child dim list. # that the dimension is omitted from the child dim list.
# itdims: Size of each index broadcast dimension, indexed (0..itdim-1)
# corners: parent coordinates of each corner, running fastest over coord index. # corners: parent coordinates of each corner, running fastest over coord index.
# (indexed 0 .. (nitems-1)*(rdim)+rdim-1) # (indexed 0 .. (nitems-1)*(rdim)+rdim-1)
# nitems: total number of list elements (product of itdims)
# itdim: number of index broadcast dimensions
# itdims: Size of each index broadcast dimension, indexed (0..itdim-1)
#
# bsize: Number of independently specified boundary conditions
# nsizes: Number of independently specified range dim sizes
# boundary: Array containing all the boundary condition specs # boundary: Array containing all the boundary condition specs
# indord: Order/size of the indexing dim (0th dim of $index)
Comp => 'PDL_Indx rdim; Comp => 'PDL_Indx rdim;
PDL_Indx nitems; PDL_Indx nitems;
PDL_Indx itdim; PDL_Indx itdim;
PDL_Indx ntsize; PDL_Indx ntsize;
PDL_Indx bsize; PDL_Indx bsize;
PDL_Indx nsizes; PDL_Indx nsizes;
PDL_Indx sizes[$COMP(rdim)]; PDL_Indx sizes[$COMP(rdim)];
PDL_Indx itdims[$COMP(itdim)]; PDL_Indx itdims[$COMP(itdim)];
PDL_Indx corners[$COMP(rdim) * $COMP(nitems)]; PDL_Indx corners[$COMP(rdim) * $COMP(nitems)];
char boundary[$COMP(rdim)]; char boundary[$COMP(rdim)];
char size_pdl_destroy; char size_pdl_destroy;
char ind_pdl_destroy; char ind_pdl_destroy;
', ',
MakeComp => <<'EOD-MakeComp', MakeComp => pp_line_numbers(__LINE__, <<'EOD-MakeComp'),
pdl *size_pdl; pdl *size_pdl;
PDL_RETERROR(PDL_err, PDL->make_physdims(ind_pdl)); PDL_RETERROR(PDL_err, PDL->make_physdims(ind_pdl));
$COMP(size_pdl_destroy) = $COMP(ind_pdl_destroy) = 0; $COMP(size_pdl_destroy) = $COMP(ind_pdl_destroy) = 0;
/* Generalized empties are ok, but not in the special 0 dim (the index vector) * / /* Generalized empties are ok, but not in the special 0 dim (the index vector) * /
if(ind_pdl->dims[0] == 0) if(ind_pdl->dims[0] == 0)
{ $CROAK("can't handle Empty indices -- call range instead"); } { $CROAK("can't handle Empty indices -- call range instead"); }
/*** /***
* Ensure that the index is a PDL_Indx. If there's no loss of information, * Ensure that the index is a PDL_Indx. If there's no loss of information,
skipping to change at line 656 skipping to change at line 652
$COMP(ind_pdl_destroy) = 1; $COMP(ind_pdl_destroy) = 1;
case PDL_SB: case PDL_B: case PDL_S: case PDL_US: case PDL_L: case PDL_UL: case PDL_LL: case PDL_ULL: case PDL_SB: case PDL_B: case PDL_S: case PDL_US: case PDL_L: case PDL_UL: case PDL_LL: case PDL_ULL:
PDL_RETERROR(PDL_err, PDL->converttype(ind_pdl,PDL_IND)); /* convert in place . */ PDL_RETERROR(PDL_err, PDL->converttype(ind_pdl,PDL_IND)); /* convert in place . */
break; break;
case PDL_IND: case PDL_IND:
PDL_RETERROR(PDL_err, PDL->make_physical(ind_pdl)); PDL_RETERROR(PDL_err, PDL->make_physical(ind_pdl));
break; break;
} }
/*** /***
* Figure sizes of the COMP arrrays and allocate them. * Figure sizes of the COMP arrays and allocate them.
*/ */
{ {
PDL_Indx i,nitems; PDL_Indx i,nitems;
$COMP(rdim) = ind_pdl->ndims ? ind_pdl->dims[0] : 1; $COMP(rdim) = ind_pdl->ndims ? ind_pdl->dims[0] : 1;
for(i=nitems=1; i < ind_pdl->ndims; i++) /* Accumulate item list size */ for(i=nitems=1; i < ind_pdl->ndims; i++) /* Accumulate item list size */
nitems *= ind_pdl->dims[i]; nitems *= ind_pdl->dims[i];
$COMP(nitems) = nitems; $COMP(nitems) = nitems;
$COMP(itdim) = ind_pdl->ndims ? ind_pdl->ndims - 1 : 0; $COMP(itdim) = ind_pdl->ndims ? ind_pdl->ndims - 1 : 0;
$DOCOMPALLOC(); $DOCOMPALLOC();
skipping to change at line 681 skipping to change at line 677
*/ */
{ {
char *bstr; char *bstr;
STRLEN blen; STRLEN blen;
bstr = SvPV(boundary_sv,blen); bstr = SvPV(boundary_sv,blen);
if(blen == 0) { if(blen == 0) {
/* If no boundary is specified then every dim gets forbidden */ /* If no boundary is specified then every dim gets forbidden */
PDL_Indx i; PDL_Indx i;
for (i=0;i<$COMP(rdim);i++) for (i=0;i<$COMP(rdim);i++)
$COMP(boundary[i]) = 0; $COMP(boundary)[i] = 0;
} else { } else {
PDL_Indx i; PDL_Indx i;
for(i=0;i<$COMP(rdim);i++) { for(i=0;i<$COMP(rdim);i++) {
switch(bstr[i < blen ? i : blen-1 ]) { switch(bstr[PDLMIN(i, blen-1)]) {
case '0': case 'f': case 'F': /* forbid */ case '0': case 'f': case 'F': /* forbid */
$COMP(boundary[i]) = 0; $COMP(boundary)[i] = 0;
break; break;
case '1': case 't': case 'T': /* truncate */ case '1': case 't': case 'T': /* truncate */
$COMP(boundary[i]) = 1; $COMP(boundary)[i] = 1;
break; break;
case '2': case 'e': case 'E': /* extend */ case '2': case 'e': case 'E': /* extend */
$COMP(boundary[i]) = 2; $COMP(boundary)[i] = 2;
break; break;
case '3': case 'p': case 'P': /* periodic */ case '3': case 'p': case 'P': /* periodic */
$COMP(boundary[i]) = 3; $COMP(boundary)[i] = 3;
break; break;
case '4': case 'm': case 'M': /* mirror */ case '4': case 'm': case 'M': /* mirror */
$COMP(boundary[i]) = 4; $COMP(boundary)[i] = 4;
break; break;
default: default:
{ {
/* No need to check if i < blen -- this will barf out the /* No need to check if i < blen -- this will barf out the
* first time it gets hit. * first time it gets hit.
*/ */
$CROAK("Unknown boundary condition '%c' in range",bstr[i]); $CROAK("Unknown boundary condition '%c' in range",bstr[i]);
} }
break; break;
} // end of switch } // end of switch
} }
} }
} }
/*** /***
* Store the sizes of the index-broadcast dims * Store the sizes of the index-broadcast dims
*/ */
{ {
PDL_Indx i; PDL_Indx i;
PDL_Indx nd = ind_pdl->ndims - 1; PDL_Indx nd = ind_pdl->ndims - 1;
for(i=0; i < nd ; i++) for(i=0; i < nd ; i++)
$COMP(itdims[i]) = ind_pdl->dims[i+1]; $COMP(itdims)[i] = ind_pdl->dims[i+1];
} }
/*** /***
* Check and condition the size ndarray, and store sizes of the ranges * Check and condition the size ndarray, and store sizes of the ranges
*/ */
{ {
PDL_Indx i,ntsize; PDL_Indx i,ntsize;
if( (size == NULL) || (size == &PL_sv_undef) ) { if( (size_sv == NULL) || !SvOK(size_sv) ) {
// NO size was passed in, not normally executed even if you passed in no siz e to range(), // NO size was passed in, not normally executed even if you passed in no siz e to range(),
// as range() generates a size array... // as range() generates a size array...
for(i=0;i<$COMP(rdim);i++) for(i=0;i<$COMP(rdim);i++)
$COMP(sizes[i]) = 0; $COMP(sizes)[i] = 0;
} else { } else {
/* Normal case with sizes present in a PDL */ /* Normal case with sizes present in a PDL */
if(!(size_pdl = PDL->SvPDLV(size))) /* assignment */ if(!(size_pdl = PDL->SvPDLV(size_sv))) /* assignment */
$CROAK("Unable to convert size to a PDL"); $CROAK("Unable to convert size to a PDL");
if(size_pdl->nvals == 0) { if(size_pdl->nvals == 0) {
// no values in the size_pdl - Empty or Null. Just copy 0s to all the ran ge dims // no values in the size_pdl - Empty or Null. Just copy 0s to all the ran ge dims
for(i=0;i<$COMP(rdim);i++) for(i=0;i<$COMP(rdim);i++)
$COMP(sizes[i]) = 0; $COMP(sizes)[i] = 0;
} else { } else {
// Convert size PDL to PDL_IND to support indices // Convert size PDL to PDL_IND to support indices
switch(size_pdl->datatype) { switch(size_pdl->datatype) {
default: /* Most types: */ default: /* Most types: */
size_pdl = PDL->hard_copy(size_pdl); /* copy and fall through */ size_pdl = PDL->hard_copy(size_pdl); /* copy and fall through */
if (!size_pdl) $CROAK("Error in hard_copy"); if (!size_pdl) $CROAK("Error in hard_copy");
$COMP(size_pdl_destroy) = 1; $COMP(size_pdl_destroy) = 1;
case PDL_SB: case PDL_B: case PDL_S: case PDL_US: case PDL_L: case PDL_UL: case PDL_LL: case PDL_ULL: case PDL_SB: case PDL_B: case PDL_S: case PDL_US: case PDL_L: case PDL_UL: case PDL_LL: case PDL_ULL:
skipping to change at line 767 skipping to change at line 763
case PDL_IND: case PDL_IND:
PDL_RETERROR(PDL_err, PDL->make_physical(size_pdl)); PDL_RETERROR(PDL_err, PDL->make_physical(size_pdl));
break; break;
} }
$COMP(nsizes) = size_pdl->nvals; /* Store for later permissiveness check * / $COMP(nsizes) = size_pdl->nvals; /* Store for later permissiveness check * /
/* Copy the sizes, or die if they're the wrong shape */ /* Copy the sizes, or die if they're the wrong shape */
if(size_pdl->nvals == 1) { if(size_pdl->nvals == 1) {
for(i=0;i<$COMP(rdim);i++) { for(i=0;i<$COMP(rdim);i++) {
$COMP(sizes[i]) = *((PDL_Indx *)(size_pdl->data)); $COMP(sizes)[i] = *((PDL_Indx *)(size_pdl->data));
} }
/* Check for nonnegativity of sizes. The rdim>0 mask ensures that */ /* Check for nonnegativity of sizes. The rdim>0 mask ensures that */
/* we don't barf on the Empty PDL (as an index). */ /* we don't barf on the Empty PDL (as an index). */
if( $COMP(rdim) > 0 && $COMP(sizes[0]) < 0 ) { if( $COMP(rdim) > 0 && $COMP(sizes)[0] < 0 ) {
$CROAK("Negative range size is not allowed\n"); $CROAK("Negative range size is not allowed\n");
} }
} }
else if( size_pdl->nvals <= $COMP(rdim) && size_pdl->ndims == 1) { else if( size_pdl->nvals <= $COMP(rdim) && size_pdl->ndims == 1) {
for(i=0;i<$COMP(rdim);i++) { for(i=0;i<$COMP(rdim);i++) {
$COMP(sizes[i]) = ( (i < size_pdl->nvals) ? $COMP(sizes)[i] = ( (i < size_pdl->nvals) ?
((PDL_Indx *)(size_pdl->data))[i] : ((PDL_Indx *)(size_pdl->data))[i] :
0 0
); );
if($COMP(sizes[i]) < 0) if($COMP(sizes)[i] < 0)
$CROAK("Negative range sizes are not allowed\n"); $CROAK("Negative range sizes are not allowed\n");
} }
} }
else { else {
$CROAK("Size must match index's 0th dim\n"); $CROAK("Size must match index's 0th dim\n");
} }
} /* end of nonempty size-ndarray code */ } /* end of nonempty size-ndarray code */
} /* end of defined-size-ndarray code */ } /* end of defined-size-ndarray code */
/* Insert the number of nontrivial sizes (these get output dimensions) */ /* Insert the number of nontrivial sizes (these get output dimensions) */
for(i=ntsize=0;i<$COMP(rdim);i++) for(i=ntsize=0;i<$COMP(rdim);i++)
if($COMP(sizes[i])) if($COMP(sizes)[i])
ntsize++; ntsize++;
$COMP(ntsize) = ntsize; $COMP(ntsize) = ntsize;
} }
/*** /***
* Stash coordinates of the corners * Stash coordinates of the corners
*/ */
PDL_Indx j,k,ioff, iter[$COMP(itdim)]; PDL_Indx j,k,ioff, iter[$COMP(itdim)];
/* initialize iterator to loop over index broadcasts */ /* initialize iterator to loop over index broadcasts */
skipping to change at line 829 skipping to change at line 825
/* Increment the iterator -- the test increments, the body carries. */ /* Increment the iterator -- the test increments, the body carries. */
for(k=0; k<$COMP(itdim) && (++(iter[k]))>=($COMP(itdims)[k]) ;k++) for(k=0; k<$COMP(itdim) && (++(iter[k]))>=($COMP(itdims)[k]) ;k++)
iter[k] = 0; iter[k] = 0;
} while(k<$COMP(itdim)); } while(k<$COMP(itdim));
if ($COMP(ind_pdl_destroy)) if ($COMP(ind_pdl_destroy))
PDL->destroy(ind_pdl); /* finished with our copy */ PDL->destroy(ind_pdl); /* finished with our copy */
if ($COMP(size_pdl_destroy)) if ($COMP(size_pdl_destroy))
PDL->destroy(size_pdl); /* finished with our copy */ PDL->destroy(size_pdl); /* finished with our copy */
EOD-MakeComp EOD-MakeComp
RedoDims => <<'EOD-RedoDims' , RedoDims => pp_line_numbers(__LINE__, <<'EOD-RedoDims'),
PDL_Indx stdim = $PDL(PARENT)->ndims - $COMP(rdim); PDL_Indx stdim = $PDL(PARENT)->ndims - $COMP(rdim);
PDL_Indx dim,inc; PDL_Indx dim,inc;
PDL_Indx i,rdvalid; PDL_Indx i,rdvalid;
// Speed bump for ludicrous cases // Speed bump for ludicrous cases
if( $COMP(rdim) > $PDL(PARENT)->ndims+5 && $COMP(nsizes) != $COMP(rdim)) { if( $COMP(rdim) > $PDL(PARENT)->ndims+5 && $COMP(nsizes) != $COMP(rdim)) {
$CROAK( $CROAK(
"Ludicrous number of extra dims in range index; leaving child null.\n" "Ludicrous number of extra dims in range index; leaving child null.\n"
" (%"IND_FLAG" implicit dims is > 5; index has %"IND_FLAG" dims; source h as %"IND_FLAG" dim%s.)\n" " (%"IND_FLAG" implicit dims is > 5; index has %"IND_FLAG" dims; source h as %"IND_FLAG" dim%s.)\n"
" This often means that your index PDL is incorrect.\n" " This often means that your index PDL is incorrect.\n"
skipping to change at line 858 skipping to change at line 854
stdim = 0; stdim = 0;
/* Set dimensionality of child */ /* Set dimensionality of child */
$PDL(CHILD)->ndims = $COMP(itdim) + $COMP(ntsize) + stdim; $PDL(CHILD)->ndims = $COMP(itdim) + $COMP(ntsize) + stdim;
$SETNDIMS($COMP(itdim)+$COMP(ntsize)+stdim); $SETNDIMS($COMP(itdim)+$COMP(ntsize)+stdim);
inc = 1; inc = 1;
/* Copy size dimensions to child, crunching as we go. */ /* Copy size dimensions to child, crunching as we go. */
dim = $COMP(itdim); dim = $COMP(itdim);
for(i=rdvalid=0;i<$COMP(rdim);i++) { for(i=rdvalid=0;i<$COMP(rdim);i++) {
if($COMP(sizes[i])) { if($COMP(sizes)[i]) {
rdvalid++; rdvalid++;
$PDL(CHILD)->dimincs[dim] = inc; $PDL(CHILD)->dimincs[dim] = inc;
inc *= ($PDL(CHILD)->dims[dim++] = $COMP(sizes[i])); /* assignment */ inc *= ($PDL(CHILD)->dims[dim++] = $COMP(sizes)[i]); /* assignment */
} }
} }
/* Copy index broadcast dimensions to child */ /* Copy index broadcast dimensions to child */
for(dim=0; dim<$COMP(itdim); dim++) { for(dim=0; dim<$COMP(itdim); dim++) {
$PDL(CHILD)->dimincs[dim] = inc; $PDL(CHILD)->dimincs[dim] = inc;
inc *= ($PDL(CHILD)->dims[dim] = $COMP(itdims[dim])); /* assignment */ inc *= ($PDL(CHILD)->dims[dim] = $COMP(itdims)[dim]); /* assignment */
} }
/* Copy source broadcast dimensions to child */ /* Copy source broadcast dimensions to child */
dim = $COMP(itdim) + rdvalid; dim = $COMP(itdim) + rdvalid;
for(i=0;i<stdim;i++) { for(i=0;i<stdim;i++) {
$PDL(CHILD)->dimincs[dim] = inc; $PDL(CHILD)->dimincs[dim] = inc;
inc *= ($PDL(CHILD)->dims[dim++] = $PDL(PARENT)->dims[i+$COMP(rdim)]); /* as signment */ inc *= ($PDL(CHILD)->dims[dim++] = $PDL(PARENT)->dims[i+$COMP(rdim)]); /* as signment */
} }
/* Cover bizarre case where the source PDL is empty - in that case, change */ /* Cover bizarre case where the source PDL is empty - in that case, change */
/* all non-barfing boundary conditions to truncation, since we have no data */ /* all non-barfing boundary conditions to truncation, since we have no data */
/* to reflect, extend, or mirror. */ /* to reflect, extend, or mirror. */
if($PDL(PARENT)->dims[0] == 0) { if($PDL(PARENT)->dims[0] == 0) {
for(dim=0; dim<$COMP(rdim); dim++) { for(dim=0; dim<$COMP(rdim); dim++) {
if($COMP(boundary[dim])) if($COMP(boundary)[dim])
$COMP(boundary[dim]) = 1; // force truncation $COMP(boundary)[dim] = 1; // force truncation
} }
} }
$PDL(CHILD)->datatype = $PDL(PARENT)->datatype; $PDL(CHILD)->datatype = $PDL(PARENT)->datatype;
$SETDIMS(); $SETDIMS();
EOD-RedoDims EOD-RedoDims
EquivCPOffsCode => <<'EOD-EquivCPOffsCode', EquivCPOffsCode => pp_line_numbers(__LINE__, <<'EOD-EquivCPOffsCode'),
PDL_Indx *ip; /* vector iterator */ PDL_Indx *ip; /* vector iterator */
PDL_Indx *sp; /* size vector including stdims */ PDL_Indx *sp; /* size vector including stdims */
PDL_Indx *coords; /* current coordinates */ PDL_Indx *coords; /* current coordinates */
PDL_Indx k; /* index */ PDL_Indx k; /* index */
PDL_Indx item; /* index broadcast iterator */ PDL_Indx item; /* index broadcast iterator */
PDL_Indx pdim = $PDL(PARENT)->ndims; PDL_Indx pdim = $PDL(PARENT)->ndims;
PDL_Indx rdim = $COMP(rdim); PDL_Indx rdim = $COMP(rdim);
PDL_Indx prdim = PDLMIN(rdim,pdim); PDL_Indx prdim = PDLMIN(rdim,pdim);
PDL_Indx iter2[pdim * 2 + rdim]; PDL_Indx iter2[pdim * 2 + rdim];
skipping to change at line 928 skipping to change at line 924
do { do {
PDL_Indx poff = 0; PDL_Indx poff = 0;
PDL_Indx coff; PDL_Indx coff;
PDL_Indx k2; PDL_Indx k2;
char trunc = 0; /* Flag used to skip truncation case */ char trunc = 0; /* Flag used to skip truncation case */
/* Collect and boundary-check the current N-D coords */ /* Collect and boundary-check the current N-D coords */
for(k=0; k < prdim; k++){ for(k=0; k < prdim; k++){
PDL_Indx ck = iter2[k] + $COMP(corners[ item * rdim + k ]) ; PDL_Indx ck = iter2[k] + $COMP(corners)[ item * rdim + k ] ;
/* normal case */ /* normal case */
if(ck < 0 || ck >= $PDL(PARENT)->dims[k]) { if(ck < 0 || ck >= $PDL(PARENT)->dims[k]) {
switch($COMP(boundary[k])) { switch($COMP(boundary)[k]) {
case 0: /* no boundary breakage allowed */ case 0: /* no boundary breakage allowed */
$CROAK("index out-of-bounds in range (index vector #%ld)",item); $CROAK("index out-of-bounds in range (index vector #%ld)",item);
break; break;
case 1: /* truncation */ case 1: /* truncation */
trunc = 1; trunc = 1;
break; break;
case 2: /* extension -- crop */ case 2: /* extension -- crop */
ck = (ck >= $PDL(PARENT)->dims[k]) ? $PDL(PARENT)->dims[k]-1 : 0; ck = (ck >= $PDL(PARENT)->dims[k]) ? $PDL(PARENT)->dims[k]-1 : 0;
break; break;
case 3: /* periodic -- mod it */ case 3: /* periodic -- mod it */
skipping to change at line 972 skipping to change at line 968
} }
coords[k] = ck; coords[k] = ck;
} }
/* Check extra dimensions -- pick up where k left off... */ /* Check extra dimensions -- pick up where k left off... */
for( ; k < rdim ; k++) { for( ; k < rdim ; k++) {
/* Check for indexing off the end of the dimension list */ /* Check for indexing off the end of the dimension list */
PDL_Indx ck = iter2[k] + $COMP(corners[ item * rdim + k ]) ; PDL_Indx ck = iter2[k] + $COMP(corners)[ item * rdim + k ] ;
switch($COMP(boundary[k])) { switch($COMP(boundary)[k]) {
case 0: /* No boundary breakage allowed -- nonzero corners cause barfa ge */ case 0: /* No boundary breakage allowed -- nonzero corners cause barfa ge */
if(ck != 0) if(ck != 0)
$CROAK("Too many dims in range index (and you've forbidden bounda ry violations)"); $CROAK("Too many dims in range index (and you've forbidden bounda ry violations)");
break; break;
case 1: /* truncation - just truncate if the corner is nonzero */ case 1: /* truncation - just truncate if the corner is nonzero */
trunc |= (ck != 0); trunc |= (ck != 0);
break; break;
case 2: /* extension -- ignore the corner (same as 3) */ case 2: /* extension -- ignore the corner (same as 3) */
case 3: /* periodic -- ignore the corner */ case 3: /* periodic -- ignore the corner */
case 4: /* mirror -- ignore the corner */ case 4: /* mirror -- ignore the corner */
skipping to change at line 1002 skipping to change at line 998
/* Find offsets into the child and parent arrays, from the N-D coords */ /* Find offsets into the child and parent arrays, from the N-D coords */
/* Note we only loop over real source dims (prdim) to accumulate -- */ /* Note we only loop over real source dims (prdim) to accumulate -- */
/* because the offset is trivial and/or we're truncating for virtual */ /* because the offset is trivial and/or we're truncating for virtual */
/* dims caused by permissive ranging. */ /* dims caused by permissive ranging. */
coff = $PDL(CHILD)->dimincs[0] * item; coff = $PDL(CHILD)->dimincs[0] * item;
for(k2 = $COMP(itdim), poff = k = 0; for(k2 = $COMP(itdim), poff = k = 0;
k < prdim; k < prdim;
k++) { k++) {
poff += coords[k]*$PDL(PARENT)->dimincs[k]; poff += coords[k]*$PDL(PARENT)->dimincs[k];
if($COMP(sizes[k])) if($COMP(sizes)[k])
coff += iter2[k] * $PDL(CHILD)->dimincs[k2++]; coff += iter2[k] * $PDL(CHILD)->dimincs[k2++];
} }
/* Loop the copy over all the source broadcast dims (above rdim). */ /* Loop the copy over all the source broadcast dims (above rdim). */
do { do {
PDL_Indx poff1 = poff; PDL_Indx poff1 = poff;
PDL_Indx coff1 = coff; PDL_Indx coff1 = coff;
/* Accumulate the offset due to source broadcasting */ /* Accumulate the offset due to source broadcasting */
for(k2 = $COMP(itdim) + $COMP(ntsize), k = rdim; for(k2 = $COMP(itdim) + $COMP(ntsize), k = rdim;
skipping to change at line 1034 skipping to change at line 1030
/* Increment the source broadcast iterator */ /* Increment the source broadcast iterator */
for( k=$COMP(rdim); for( k=$COMP(rdim);
k < pdim && (++(iter2[k]) >= $PDL(PARENT)->dims[k]); k < pdim && (++(iter2[k]) >= $PDL(PARENT)->dims[k]);
k++) k++)
iter2[k] = 0; iter2[k] = 0;
} while(k < pdim); /* end of source-broadcast iteration */ } while(k < pdim); /* end of source-broadcast iteration */
/* Increment the in-range iterator */ /* Increment the in-range iterator */
for(k = 0; for(k = 0;
k < $COMP(rdim) && (++(iter2[k]) >= $COMP(sizes[k])); k < $COMP(rdim) && (++(iter2[k]) >= $COMP(sizes)[k]);
k++) k++)
iter2[k] = 0; iter2[k] = 0;
} while(k < $COMP(rdim)); /* end of main iteration */ } while(k < $COMP(rdim)); /* end of main iteration */
} /* end of item do loop */ } /* end of item do loop */
EOD-EquivCPOffsCode EOD-EquivCPOffsCode
); );
pp_def( pp_def(
'rld', 'rld',
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
Pars=>'indx a(n); b(n); [o]c(m);', Pars=>'indx a(n); b(n); [o]c(m);',
PMCode =><<'EOD', PMCode =>pp_line_numbers(__LINE__, <<'EOD'),
sub PDL::rld { sub PDL::rld {
my ($x,$y) = @_; my ($x,$y) = @_;
my ($c); my ($c);
if ($#_ == 2) { if ($#_ == 2) {
$c = $_[2]; $c = $_[2];
} else { } else {
# XXX Need to improve emulation of broadcasting in auto-generating c # XXX Need to improve emulation of broadcasting in auto-generating c
my ($size) = $x->sumover->max->sclr; my ($size) = $x->sumover->max->sclr;
my (@dims) = $x->dims; my (@dims) = $x->dims;
shift @dims; shift @dims;
$c = $y->zeroes($size,@dims); $c = $y->zeroes($size,@dims);
} }
&PDL::_rld_int($x,$y,$c); &PDL::_rld_int($x,$y,$c);
$c; $c;
} }
EOD EOD
Code=>' Code => pp_line_numbers(__LINE__, <<'EOF'),
PDL_Indx i,j=0,an; PDL_Indx i,j=0,an;
$GENERIC(b) bv; $GENERIC(b) bv;
loop (n) %{ loop (n) %{
an = $a(); an = $a();
bv = $b(); bv = $b();
for (i=0;i<an;i++) { for (i=0;i<an;i++) {
$c(m=>j) = bv; $c(m=>j) = bv;
j++; j++;
} }
%}', %}
EOF
Doc => <<'EOD' Doc => <<'EOD'
=for ref =for ref
Run-length decode a vector Run-length decode a vector
Given a vector C<$x> of the numbers of instances of values C<$y>, run-length Given a vector C<$x> of the numbers of instances of values C<$y>, run-length
decode to C<$c>. decode to C<$c>.
=for example =for example
skipping to change at line 1096 skipping to change at line 1093
=cut =cut
EOD EOD
); );
pp_def( pp_def(
'rle', 'rle',
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
Pars=>'c(n); indx [o]a(m); [o]b(m);', Pars=>'c(n); indx [o]a(m); [o]b(m);',
RedoDimsCode=>'$SIZE(m)=$SIZE(n);', RedoDimsCode=>'$SIZE(m)=$SIZE(n);',
PMCode=><<'EOC', PMCode =>pp_line_numbers(__LINE__, <<'EOC'),
sub PDL::rle { sub PDL::rle {
my $c = shift; my $c = shift;
my ($x,$y) = @_==2 ? @_ : (null,null); my ($x,$y) = @_==2 ? @_ : (null,null);
PDL::_rle_int($c,$x,$y); PDL::_rle_int($c,$x,$y);
my $max_ind = ($c->ndims<2) ? ($x!=0)->sumover-1 : my $max_ind = ($c->ndims<2) ? ($x!=0)->sumover-1 :
($x!=0)->clump(1..$x->ndims-1)->sumover->max->sc lr-1; ($x!=0)->clump(1..$x->ndims-1)->sumover->max->sc lr-1;
return ($x->slice("0:$max_ind"),$y->slice("0:$max_ind")); return ($x->slice("0:$max_ind"),$y->slice("0:$max_ind"));
} }
EOC EOC
Code=>' Code =>pp_line_numbers(__LINE__, <<'EOF'),
PDL_Indx j=0; PDL_Indx j=0;
$GENERIC(c) cv, clv; $GENERIC(c) cv, clv;
clv = $c(n=>0); clv = $c(n=>0);
$b(m=>0) = clv; $b(m=>0) = clv;
$a(m=>0) = 0; $a(m=>0) = 0;
loop (n) %{ loop (n) %{
cv = $c(); cv = $c();
if (cv == clv) { if (cv == clv) {
$a(m=>j)++; $a(m=>j)++;
} else { } else {
j++; j++;
$b(m=>j) = clv = cv; $b(m=>j) = clv = cv;
$a(m=>j) = 1; $a(m=>j) = 1;
} }
%} %}
for (j++;j<$SIZE(m);j++) { for (j++;j<$SIZE(m);j++) {
$a(m=>j) = 0; $a(m=>j) = 0;
$b(m=>j) = 0; $b(m=>j) = 0;
} }
', EOF
Doc => <<'EOD' Doc => <<'EOD'
=for ref =for ref
Run-length encode a vector Run-length encode a vector
Given vector C<$c>, generate a vector C<$x> with the number of each Given vector C<$c>, generate a vector C<$x> with the number of each
element, and a vector C<$y> of the unique values. New in PDL 2.017, element, and a vector C<$y> of the unique values. New in PDL 2.017,
only the elements up to the first instance of C<0> in C<$x> are only the elements up to the first instance of C<0> in C<$x> are
returned, which makes the common use case of a 1-dimensional C<$c> simpler. returned, which makes the common use case of a 1-dimensional C<$c> simpler.
For broadcast operation, C<$x> and C<$y> will be large enough For broadcast operation, C<$x> and C<$y> will be large enough
skipping to change at line 1172 skipping to change at line 1169
$yprob1x = $y->slice('-1:0')->double / $y->slice('(-1)'); $yprob1x = $y->slice('-1:0')->double / $y->slice('(-1)');
$z = cat($x1, 1 / $yprob1x**$nrolls)->transpose; $z = cat($x1, 1 / $yprob1x**$nrolls)->transpose;
=cut =cut
EOD EOD
); );
pp_def('rlevec', pp_def('rlevec',
Pars => "c(M,N); indx [o]a(N); [o]b(M,N)", Pars => "c(M,N); indx [o]a(N); [o]b(M,N)",
Code =><<'EOC', Code =>pp_line_numbers(__LINE__, <<'EOC'),
PDL_Indx cn,bn=0, sn=$SIZE(N), matches; PDL_Indx cn,bn=0, sn=$SIZE(N), matches;
loop (M) %{ $b(N=>0)=$c(N=>0); %} loop (M) %{ $b(N=>0)=$c(N=>0); %}
$a(N=>0) = 1; $a(N=>0) = 1;
for (cn=1; cn<sn; cn++) { for (cn=1; cn<sn; cn++) {
matches=1; matches=1;
loop (M) %{ loop (M) %{
if ($c(N=>cn) != $b(N=>bn)) { if ($c(N=>cn) != $b(N=>bn)) {
matches=0; matches=0;
break; break;
} }
skipping to change at line 1221 skipping to change at line 1218
over a 1d PDL. over a 1d PDL.
See also: L</rle>, L<PDL::Ufunc/qsortvec>, L<PDL::Primitive/uniqvec> See also: L</rle>, L<PDL::Ufunc/qsortvec>, L<PDL::Primitive/uniqvec>
Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>. Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
EOD EOD
); );
pp_def('rldvec', pp_def('rldvec',
Pars => 'indx a(N); b(M,N); [o]c(M,N)', Pars => 'indx a(N); b(M,N); [o]c(M,N)',
PMCode=><<'EOC', PMCode =>pp_line_numbers(__LINE__, <<'EOC'),
sub PDL::rldvec { sub PDL::rldvec {
my ($a,$b,$c) = @_; my ($a,$b,$c) = @_;
if (!defined($c)) { if (!defined($c)) {
# XXX Need to improve emulation of threading in auto-generating c # XXX Need to improve emulation of threading in auto-generating c
my ($rowlen) = $b->dim(0); my ($rowlen) = $b->dim(0);
my ($size) = $a->sumover->max; my ($size) = $a->sumover->max;
my (@dims) = $a->dims; my (@dims) = $a->dims;
shift(@dims); shift(@dims);
$c = $b->zeroes($b->type,$rowlen,$size,@dims); $c = $b->zeroes($b->type,$rowlen,$size,@dims);
} }
&PDL::_rldvec_int($a,$b,$c); &PDL::_rldvec_int($a,$b,$c);
return $c; return $c;
} }
EOC EOC
Code =><<'EOC', Code =>pp_line_numbers(__LINE__, <<'EOC'),
PDL_Indx i,nrows,bn,cn=0, sn=$SIZE(N); PDL_Indx i,nrows,bn,cn=0, sn=$SIZE(N);
for (bn=0; bn<sn; bn++) { for (bn=0; bn<sn; bn++) {
nrows = $a(N=>bn); nrows = $a(N=>bn);
for (i=0; i<nrows; i++) { for (i=0; i<nrows; i++) {
loop (M) %{ $c(N=>cn) = $b(N=>bn); %} loop (M) %{ $c(N=>cn) = $b(N=>bn); %}
cn++; cn++;
} }
} }
EOC EOC
Doc =><<'EOD' Doc =><<'EOD'
skipping to change at line 1263 skipping to change at line 1260
Can be used together with clump() to run-length decode "values" of arbitrary dim ensions. Can be used together with clump() to run-length decode "values" of arbitrary dim ensions.
See also: L</rld>. See also: L</rld>.
Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>. Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
EOD EOD
); );
pp_def('rleseq', pp_def('rleseq',
Pars => "c(N); indx [o]a(N); [o]b(N)", Pars => "c(N); indx [o]a(N); [o]b(N)",
Code=><<'EOC', Code =>pp_line_numbers(__LINE__, <<'EOC'),
PDL_Indx j=0, sizeN=$SIZE(N); PDL_Indx j=0, sizeN=$SIZE(N);
$GENERIC(c) coff; $GENERIC(c) coff;
coff = $c(N=>0); coff = $c(N=>0);
$b(N=>0) = coff; $b(N=>0) = coff;
$a(N=>0) = 0; $a(N=>0) = 0;
loop (N) %{ loop (N) %{
if ($c() == coff+$a(N=>j)) { if ($c() == coff+$a(N=>j)) {
$a(N=>j)++; $a(N=>j)++;
} else { } else {
j++; j++;
skipping to change at line 1300 skipping to change at line 1297
and a vector $b containing the subsequence offsets. and a vector $b containing the subsequence offsets.
As for rle(), only the elements up to the first instance of 0 in $a should be co nsidered. As for rle(), only the elements up to the first instance of 0 in $a should be co nsidered.
See also L</rle>. See also L</rle>.
Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>. Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
EOD EOD
); );
pp_def('rldseq', pp_def('rldseq',
Pars => 'indx a(N); b(N); [o]c(M)', Pars => 'indx a(N); b(N); [o]c(M)',
PMCode=><<'EOC', PMCode =>pp_line_numbers(__LINE__, <<'EOC'),
sub PDL::rldseq { sub PDL::rldseq {
my ($a,$b,$c) = @_; my ($a,$b,$c) = @_;
if (!defined($c)) { if (!defined($c)) {
my $size = $a->sumover->max; my $size = $a->sumover->max;
my (@dims) = $a->dims; my (@dims) = $a->dims;
shift(@dims); shift(@dims);
$c = $b->zeroes($b->type,$size,@dims); $c = $b->zeroes($b->type,$size,@dims);
} }
&PDL::_rldseq_int($a,$b,$c); &PDL::_rldseq_int($a,$b,$c);
return $c; return $c;
} }
EOC EOC
Code =><<'EOC', Code =>pp_line_numbers(__LINE__, <<'EOC'),
size_t mi=0; size_t mi=0;
loop (N) %{ loop (N) %{
size_t len = $a(), li; size_t len = $a(), li;
for (li=0; li < len; ++li, ++mi) { for (li=0; li < len; ++li, ++mi) {
$c(M=>mi) = $b() + li; $c(M=>mi) = $b() + li;
} }
%} %}
EOC EOC
Doc =><<'EOD' Doc =><<'EOD'
=for ref =for ref
skipping to change at line 1427 skipping to change at line 1424
return $data; return $data;
} }
EOF EOF
# the perl wrapper clump is now defined in Core.pm # the perl wrapper clump is now defined in Core.pm
# this is just the low level interface # this is just the low level interface
pp_def( pp_def(
'_clump_int', '_clump_int',
OtherPars => 'PDL_Indx n', OtherPars => 'PDL_Indx n',
P2Child => 1, P2Child => 1,
RedoDims => 'PDL_Indx i, d1; RedoDims => pp_line_numbers(__LINE__, <<'EOF'),
/* truncate overly long clumps to just clump existing dimensions */ /* truncate overly long clumps to just clump existing dimensions */
if($COMP(n) > $PDL(PARENT)->ndims) if($COMP(n) > $PDL(PARENT)->ndims)
$COMP(n) = $PDL(PARENT)->ndims; $COMP(n) = $PDL(PARENT)->ndims;
if($COMP(n) < -1) if($COMP(n) < -1)
$COMP(n) = $PDL(PARENT)->ndims + $COMP(n) + 1; $COMP(n) = $PDL(PARENT)->ndims + $COMP(n) + 1;
PDL_Indx nrem = ($COMP(n) == -1 ? $PDL(PARENT)->broadcastids[0] : $COMP(n)); PDL_Indx nrem = ($COMP(n) == -1 ? $PDL(PARENT)->broadcastids[0] : $COMP(n));
$SETNDIMS($PDL(PARENT)->ndims - nrem + 1); $SETNDIMS($PDL(PARENT)->ndims - nrem + 1);
d1=1; PDL_Indx i, d1=1;
for(i=0; i<nrem; i++) { for(i=0; i<nrem; i++) {
d1 *= $PDL(PARENT)->dims[i]; d1 *= $PDL(PARENT)->dims[i];
} }
$PDL(CHILD)->dims[0] = d1; $PDL(CHILD)->dims[0] = d1;
for(; i<$PDL(PARENT)->ndims; i++) { for(; i<$PDL(PARENT)->ndims; i++) {
$PDL(CHILD)->dims[i-nrem+1] = $PDL(PARENT)->dims[i]; $PDL(CHILD)->dims[i-nrem+1] = $PDL(PARENT)->dims[i];
} }
$SETDIMS(); $SETDIMS();
$SETDELTABROADCASTIDS(1-nrem); $SETDELTABROADCASTIDS(1-nrem);
', EOF
EquivCPOffsCode => ' EquivCPOffsCode => pp_line_numbers(__LINE__, <<'EOF'),
PDL_Indx i; PDL_Indx i;
for(i=0; i<$PDL(CHILD)->nvals; i++) { for(i=0; i<$PDL(CHILD)->nvals; i++) {
$EQUIVCPOFFS(i,i); $EQUIVCPOFFS(i,i);
} }
', EOF
TwoWay => 1, TwoWay => 1,
Doc => 'internal', Doc => 'internal',
); );
pp_def( pp_def(
'xchg', 'xchg',
OtherPars => 'PDL_Indx n1; PDL_Indx n2;', OtherPars => 'PDL_Indx n1; PDL_Indx n2;',
TwoWay => 1, TwoWay => 1,
P2Child => 1, P2Child => 1,
AffinePriv => 1, AffinePriv => 1,
EquivDimCheck => 'if ($COMP(n1) <0) EquivDimCheck =>pp_line_numbers(__LINE__, <<'EOF'),
$COMP(n1) += $PDL(PARENT)->broadcastids[0]; if ($COMP(n1) <0)
if ($COMP(n2) <0) $COMP(n1) += $PDL(PARENT)->broadcastids[0];
$COMP(n2) += $PDL(PARENT)->broadcastids[0]; if ($COMP(n2) <0)
$COMP(n2) += $PDL(PARENT)->broadcastids[0];
if (PDLMIN($COMP(n1),$COMP(n2)) <0 || if (PDLMIN($COMP(n1),$COMP(n2)) <0 ||
PDLMAX($COMP(n1),$COMP(n2)) >= $PDL(PARENT)->broadcastids[0]) PDLMAX($COMP(n1),$COMP(n2)) >= $PDL(PARENT)->broadcastids[0])
$CROAK("One of dims %"IND_FLAG", %"IND_FLAG" out of range: shoul d be 0<=dim<%"IND_FLAG"", $CROAK("One of dims %"IND_FLAG", %"IND_FLAG" out of range: shoul d be 0<=dim<%"IND_FLAG"",
$COMP(n1),$COMP(n2),$PDL(PARENT)->broadcastids[0]);', $COMP(n1),$COMP(n2),$PDL(PARENT)->broadcastids[0]);
EquivPDimExpr => '(($CDIM == $COMP(n1)) ? $COMP(n2) : ($CDIM == $COMP(n2 EOF
)) ? $COMP(n1) : $CDIM)', EquivPDimExpr =>pp_line_numbers(__LINE__, <<'EOF'),
(($CDIM == $COMP(n1)) ? $COMP(n2) : ($CDIM == $COMP(n2)) ? $COMP(n1) :
$CDIM)
EOF
Doc => <<'EOD', Doc => <<'EOD',
=for ref =for ref
exchange two dimensions exchange two dimensions
Negative dimension indices count from the end. Negative dimension indices count from the end.
The command The command
=for example =for example
skipping to change at line 1497 skipping to change at line 1497
creates C<$y> to be like C<$x> except that the dimensions 2 and 3 creates C<$y> to be like C<$x> except that the dimensions 2 and 3
are exchanged with each other i.e. are exchanged with each other i.e.
$y->at(5,3,2,8) == $x->at(5,3,8,2) $y->at(5,3,2,8) == $x->at(5,3,8,2)
=cut =cut
EOD EOD
); );
pp_addpm(<< 'EOD'); pp_addpm(<<'EOD');
=head2 reorder =head2 reorder
=for ref =for ref
Re-orders the dimensions of a PDL based on the supplied list. Re-orders the dimensions of a PDL based on the supplied list.
Similar to the L</xchg> method, this method re-orders the dimensions Similar to the L</xchg> method, this method re-orders the dimensions
of a PDL. While the L</xchg> method swaps the position of two dimensions, of a PDL. While the L</xchg> method swaps the position of two dimensions,
the reorder method can change the positions of many dimensions at the reorder method can change the positions of many dimensions at
skipping to change at line 1631 skipping to change at line 1631
#define EQUIVDIM(dima,dimb,cdim,inc) \ #define EQUIVDIM(dima,dimb,cdim,inc) \
((cdim < PDLMIN(dima,dimb) || cdim > PDLMAX(dima,dimb)) ? \ ((cdim < PDLMIN(dima,dimb) || cdim > PDLMAX(dima,dimb)) ? \
cdim : ((cdim == dimb) ? dima : cdim + inc)) cdim : ((cdim == dimb) ? dima : cdim + inc))
EOF EOF
pp_def( pp_def(
'mv', 'mv',
OtherPars => 'PDL_Indx n1; PDL_Indx n2;', OtherPars => 'PDL_Indx n1; PDL_Indx n2;',
TwoWay => 1, TwoWay => 1,
P2Child => 1, P2Child => 1,
AffinePriv => 1, AffinePriv => 1,
EquivDimCheck => 'if ($COMP(n1) <0) EquivDimCheck =>pp_line_numbers(__LINE__, <<'EOF'),
$COMP(n1) += $PDL(PARENT)->broadcastids[0]; if ($COMP(n1) <0)
if ($COMP(n2) <0) $COMP(n1) += $PDL(PARENT)->broadcastids[0];
$COMP(n2) += $PDL(PARENT)->broadcastids[0]; if ($COMP(n2) <0)
$COMP(n2) += $PDL(PARENT)->broadcastids[0];
if (PDLMIN($COMP(n1),$COMP(n2)) <0 || if (PDLMIN($COMP(n1),$COMP(n2)) <0 ||
PDLMAX($COMP(n1),$COMP(n2)) >= $PDL(PARENT)->broadcastids[0]) PDLMAX($COMP(n1),$COMP(n2)) >= $PDL(PARENT)->broadcastids[0])
$CROAK("One of dims %"IND_FLAG", %"IND_FLAG" out of range: shoul d be 0<=dim<%"IND_FLAG"", $CROAK("One of dims %"IND_FLAG", %"IND_FLAG" out of range: shoul d be 0<=dim<%"IND_FLAG"",
$COMP(n1),$COMP(n2),$PDL(PARENT)->broadcastids[0]);', $COMP(n1),$COMP(n2),$PDL(PARENT)->broadcastids[0]);
EquivPDimExpr => '( EOF
EquivPDimExpr =>pp_line_numbers(__LINE__, <<'EOF'),
(
$COMP(n1) == $COMP(n2) ? $CDIM : $COMP(n1) == $COMP(n2) ? $CDIM :
$COMP(n1) < $COMP(n2) ? EQUIVDIM($COMP(n1),$COMP(n2),$CDIM,1) : $COMP(n1) < $COMP(n2) ? EQUIVDIM($COMP(n1),$COMP(n2),$CDIM,1) :
EQUIVDIM($COMP(n1),$COMP(n2),$CDIM,-1) EQUIVDIM($COMP(n1),$COMP(n2),$CDIM,-1)
)', )
Doc => << 'EOD', EOF
Doc => <<'EOD',
=for ref =for ref
move a dimension to another position move a dimension to another position
The command The command
=for example =for example
$y = $x->mv(4,1); $y = $x->mv(4,1);
skipping to change at line 1666 skipping to change at line 1670
place 1, so: place 1, so:
$y->at(1,2,3,4,5,6) == $x->at(1,5,2,3,4,6); $y->at(1,2,3,4,5,6) == $x->at(1,5,2,3,4,6);
The other dimensions are moved accordingly. The other dimensions are moved accordingly.
Negative dimension indices count from the end. Negative dimension indices count from the end.
=cut =cut
EOD EOD
); );
pp_addhdr << 'EOH'; pp_addhdr <<'EOH';
#define sign(x) ( (x) < 0 ? -1 : 1) #define sign(x) ( (x) < 0 ? -1 : 1)
EOH EOH
pp_addpm(<<'EOD' pp_addpm(<<'EOD'
=head2 using =head2 using
=for ref =for ref
Returns array of column numbers requested Returns array of column numbers requested
skipping to change at line 1729 skipping to change at line 1733
TwoWay => 1, TwoWay => 1,
AffinePriv => 1, AffinePriv => 1,
OtherPars => 'PDL_Indx whichdims[]', OtherPars => 'PDL_Indx whichdims[]',
MakeComp => pp_line_numbers(__LINE__-1, ' MakeComp => pp_line_numbers(__LINE__-1, '
if ($COMP(whichdims_count) < 1) if ($COMP(whichdims_count) < 1)
$CROAK("must have at least 1 dimension"); $CROAK("must have at least 1 dimension");
qsort($COMP(whichdims), $COMP(whichdims_count), sizeof(PDL_Indx) , qsort($COMP(whichdims), $COMP(whichdims_count), sizeof(PDL_Indx) ,
cmp_pdll); cmp_pdll);
'), '),
RedoDims => pp_line_numbers(__LINE__-1, ' RedoDims => pp_line_numbers(__LINE__-1, '
PDL_Indx nthp,nthc,nthd, cd = $COMP(whichdims[0]); PDL_Indx nthp,nthc,nthd, cd = $COMP(whichdims)[0];
$SETNDIMS($PDL(PARENT)->ndims-$COMP(whichdims_count)+1); $SETNDIMS($PDL(PARENT)->ndims-$COMP(whichdims_count)+1);
$DOPRIVALLOC(); $DOPRIVALLOC();
$PRIV(offs) = 0; $PRIV(offs) = 0;
if ($COMP(whichdims)[$COMP(whichdims_count)-1] >= $PDL(PARENT)-> ndims || if ($COMP(whichdims)[$COMP(whichdims_count)-1] >= $PDL(PARENT)-> ndims ||
$COMP(whichdims)[0] < 0) $COMP(whichdims)[0] < 0)
$CROAK("dim out of range"); $CROAK("dim out of range");
nthd=0; nthc=0; nthd=0; nthc=0;
for(nthp=0; nthp<$PDL(PARENT)->ndims; nthp++) for(nthp=0; nthp<$PDL(PARENT)->ndims; nthp++)
if (nthd < $COMP(whichdims_count) && if (nthd < $COMP(whichdims_count) &&
nthp == $COMP(whichdims)[nthd]) { nthp == $COMP(whichdims)[nthd]) {
skipping to change at line 1763 skipping to change at line 1767
$PDL(PARENT)->dims[nthp]); $PDL(PARENT)->dims[nthp]);
} }
$PRIV(incs)[cd] += $PDL(PARENT)->dimincs[nthp]; $PRIV(incs)[cd] += $PDL(PARENT)->dimincs[nthp];
} else { } else {
$PRIV(incs)[nthc] = $PDL(PARENT)->dimincs[nthp]; $PRIV(incs)[nthc] = $PDL(PARENT)->dimincs[nthp];
$PDL(CHILD)->dims[nthc] = $PDL(PARENT)->dims[nth p]; $PDL(CHILD)->dims[nthc] = $PDL(PARENT)->dims[nth p];
nthc++; nthc++;
} }
$SETDIMS(); $SETDIMS();
'), '),
PMCode => << 'EOD', PMCode =>pp_line_numbers(__LINE__, <<'EOD'),
sub PDL::diagonal { shift->_diagonal_int(my $o=PDL->null, \@_); $o } sub PDL::diagonal { shift->_diagonal_int(my $o=PDL->null, \@_); $o }
EOD EOD
Doc => << 'EOD', Doc => <<'EOD',
=for ref =for ref
Returns the multidimensional diagonal over the specified dimensions. Returns the multidimensional diagonal over the specified dimensions.
The diagonal is placed at the first (by number) dimension that is The diagonal is placed at the first (by number) dimension that is
diagonalized. diagonalized.
The other diagonalized dimensions are removed. So if C<$x> has dimensions The other diagonalized dimensions are removed. So if C<$x> has dimensions
C<(5,3,5,4,6,5)> then after C<(5,3,5,4,6,5)> then after
=for usage =for usage
skipping to change at line 1854 skipping to change at line 1858
C<$step> and C<$nlags> must be positive. C<$nthdim> can be C<$step> and C<$nlags> must be positive. C<$nthdim> can be
negative and will then be counted from the last dim backwards negative and will then be counted from the last dim backwards
in the usual way (-1 = last dim). in the usual way (-1 = last dim).
=cut =cut
EOD EOD
P2Child => 1, P2Child => 1,
TwoWay => 1, TwoWay => 1,
AffinePriv => 1, AffinePriv => 1,
OtherPars => join('', map "PDL_Indx $_;", qw(nthdim step n)), OtherPars => join('', map "PDL_Indx $_;", qw(nthdim step n)),
RedoDims => ' RedoDims => pp_line_numbers(__LINE__, <<'EOF'),
PDL_Indx i; PDL_Indx i;
if ($COMP(nthdim) < 0) /* the usual conventions */ if ($COMP(nthdim) < 0) /* the usual conventions */
$COMP(nthdim) += $PDL(PARENT)->ndims; $COMP(nthdim) += $PDL(PARENT)->ndims;
if ($COMP(nthdim) < 0 || $COMP(nthdim) >= $PDL(PARENT)->ndims) if ($COMP(nthdim) < 0 || $COMP(nthdim) >= $PDL(PARENT)->ndims)
$CROAK("dim out of range"); $CROAK("dim out of range");
if ($COMP(n) < 1) if ($COMP(n) < 1)
$CROAK("number of lags must be positive"); $CROAK("number of lags must be positive");
if ($COMP(step) < 1) if ($COMP(step) < 1)
$CROAK("step must be positive"); $CROAK("step must be positive");
$PRIV(offs) = 0; $PRIV(offs) = 0;
skipping to change at line 1884 skipping to change at line 1888
$PDL(CHILD)->dims[i+1] = $COMP(n); $PDL(CHILD)->dims[i+1] = $COMP(n);
$PRIV(incs)[i] = ($PDL(PARENT)->dimincs[i]); $PRIV(incs)[i] = ($PDL(PARENT)->dimincs[i]);
$PRIV(incs)[i+1] = - $PDL(PARENT)->dimincs[i] * $COMP(step); $PRIV(incs)[i+1] = - $PDL(PARENT)->dimincs[i] * $COMP(step);
$PRIV(offs) += ($PDL(CHILD)->dims[i+1] - 1) * (-$PRIV(incs)[i+1] ); $PRIV(offs) += ($PDL(CHILD)->dims[i+1] - 1) * (-$PRIV(incs)[i+1] );
i++; i++;
for(; i<$PDL(PARENT)->ndims; i++) { for(; i<$PDL(PARENT)->ndims; i++) {
$PDL(CHILD)->dims[i+1] = $PDL(PARENT)->dims[i]; $PDL(CHILD)->dims[i+1] = $PDL(PARENT)->dims[i];
$PRIV(incs)[i+1] = $PDL(PARENT)->dimincs[i]; $PRIV(incs)[i+1] = $PDL(PARENT)->dimincs[i];
} }
$SETDIMS(); $SETDIMS();
' EOF
); );
pp_def( pp_def(
'splitdim', 'splitdim',
Doc => <<'EOD', Doc => <<'EOD',
=for ref =for ref
Splits a dimension in the parent ndarray (opposite of L<clump|PDL::Core/clump>). Splits a dimension in the parent ndarray (opposite of L<clump|PDL::Core/clump>).
As of 2.076, throws exception if non-divisible C<nsp> given, and can As of 2.076, throws exception if non-divisible C<nsp> given, and can
give negative C<nthdim> which then counts backwards. give negative C<nthdim> which then counts backwards.
skipping to change at line 1913 skipping to change at line 1917
$y->at(6,4,m,n,3,6) == $x->at(6,4,m+3*n) $y->at(6,4,m,n,3,6) == $x->at(6,4,m+3*n)
is always true (C<m> has to be less than 3). is always true (C<m> has to be less than 3).
=cut =cut
EOD EOD
P2Child => 1, P2Child => 1,
TwoWay => 1, TwoWay => 1,
OtherPars => join('', map "PDL_Indx $_;", qw(nthdim nsp)), OtherPars => join('', map "PDL_Indx $_;", qw(nthdim nsp)),
AffinePriv => 1, AffinePriv => 1,
RedoDims => ' RedoDims => pp_line_numbers(__LINE__, <<'EOF'),
PDL_Indx i = $COMP(nthdim); PDL_Indx i = $COMP(nthdim);
PDL_Indx nsp = $COMP(nsp); PDL_Indx nsp = $COMP(nsp);
if(nsp == 0) {$CROAK("Cannot split to 0\n");} if(nsp == 0) {$CROAK("Cannot split to 0\n");}
if (i < 0) if (i < 0)
i = $COMP(nthdim) += $PDL(PARENT)->ndims; i = $COMP(nthdim) += $PDL(PARENT)->ndims;
if (i < 0 || i >= $PDL(PARENT)->ndims) if (i < 0 || i >= $PDL(PARENT)->ndims)
$CROAK("nthdim %"IND_FLAG" after adjusting for negative must not be negative or greater or equal to number of dims %"IND_FLAG"\n", $CROAK("nthdim %"IND_FLAG" after adjusting for negative must not be negative or greater or equal to number of dims %"IND_FLAG"\n",
i, $PDL(PARENT)->ndims); i, $PDL(PARENT)->ndims);
if (nsp > $PDL(PARENT)->dims[i]) if (nsp > $PDL(PARENT)->dims[i])
$CROAK("nsp %"IND_FLAG" cannot be greater than dim %"IND _FLAG"\n", $CROAK("nsp %"IND_FLAG" cannot be greater than dim %"IND _FLAG"\n",
skipping to change at line 1945 skipping to change at line 1949
$PDL(CHILD)->dims[i] = $COMP(nsp); $PDL(CHILD)->dims[i] = $COMP(nsp);
$PDL(CHILD)->dims[i+1] = $PDL(PARENT)->dims[i] / $COMP(nsp); $PDL(CHILD)->dims[i+1] = $PDL(PARENT)->dims[i] / $COMP(nsp);
$PRIV(incs)[i] = $PDL(PARENT)->dimincs[i]; $PRIV(incs)[i] = $PDL(PARENT)->dimincs[i];
$PRIV(incs)[i+1] = $PDL(PARENT)->dimincs[i] * $COMP(nsp); $PRIV(incs)[i+1] = $PDL(PARENT)->dimincs[i] * $COMP(nsp);
i++; i++;
for(; i<$PDL(PARENT)->ndims; i++) { for(; i<$PDL(PARENT)->ndims; i++) {
$PDL(CHILD)->dims[i+1] = $PDL(PARENT)->dims[i]; $PDL(CHILD)->dims[i+1] = $PDL(PARENT)->dims[i];
$PRIV(incs)[i+1] = $PDL(PARENT)->dimincs[i]; $PRIV(incs)[i+1] = $PDL(PARENT)->dimincs[i];
} }
$SETDIMS(); $SETDIMS();
', EOF
); );
my $rotate_code = ' my $rotate_code = pp_line_numbers(__LINE__, <<'EOF');
PDL_Indx i,j; PDL_Indx i,j;
PDL_Indx n_size = $SIZE(n); PDL_Indx n_size = $SIZE(n);
if (n_size == 0) if (n_size == 0)
$CROAK("can not shift zero size ndarray (n_size is zero)"); $CROAK("can not shift zero size ndarray (n_size is zero)");
j = ($shift()) % n_size; j = ($shift()) % n_size;
if (j < 0) if (j < 0)
j += n_size; j += n_size;
for(i=0; i<n_size; i++,j++) { for(i=0; i<n_size; i++,j++) {
if (j == n_size) if (j == n_size)
j = 0; j = 0;
'; EOF
pp_def('rotate', pp_def('rotate',
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
Doc => <<'EOD', Doc => <<'EOD',
=for ref =for ref
Shift vector elements along with wrap. Flows data back&forth. Shift vector elements along with wrap. Flows data back&forth.
=cut =cut
EOD EOD
Pars=>'x(n); indx shift(); [oca]y(n)', Pars=>'x(n); indx shift(); [oca]y(n)',
DefaultFlow => 1, DefaultFlow => 1,
skipping to change at line 2002 skipping to change at line 2006
$y = $x->broadcastI(0,1,5); # broadcast over dims 1,5 in id 1 $y = $x->broadcastI(0,1,5); # broadcast over dims 1,5 in id 1
=cut =cut
EOD EOD
P2Child => 1, P2Child => 1,
TwoWay => 1, TwoWay => 1,
AffinePriv => 1, AffinePriv => 1,
OtherPars => "PDL_Indx id; PDL_Indx whichdims[]", OtherPars => "PDL_Indx id; PDL_Indx whichdims[]",
Comp => 'PDL_Indx nrealwhichdims', Comp => 'PDL_Indx nrealwhichdims',
MakeComp => ' MakeComp => pp_line_numbers(__LINE__, <<'EOF'),
PDL_Indx i,j; PDL_Indx i,j;
$COMP(nrealwhichdims) = 0; $COMP(nrealwhichdims) = 0;
for(i=0; i<$COMP(whichdims_count); i++) { for(i=0; i<$COMP(whichdims_count); i++) {
for(j=i+1; j<$COMP(whichdims_count); j++) for(j=i+1; j<$COMP(whichdims_count); j++)
if($COMP(whichdims[i]) == $COMP(whichdims[j]) && if($COMP(whichdims)[i] == $COMP(whichdims)[j] &&
$COMP(whichdims[i]) != -1) { $COMP(whichdims)[i] != -1) {
$CROAK("duplicate arg %"IND_FLAG" %"IND_FLAG" %" IND_FLAG"", $CROAK("duplicate arg %"IND_FLAG" %"IND_FLAG" %" IND_FLAG"",
i,j,$COMP(whichdims[i])); i,j,$COMP(whichdims)[i]);
} }
if($COMP(whichdims)[i] != -1) { if($COMP(whichdims)[i] != -1) {
$COMP(nrealwhichdims) ++; $COMP(nrealwhichdims) ++;
} }
} }
', EOF
RedoDims => ' RedoDims => pp_line_numbers(__LINE__, <<'EOF'),
PDL_Indx nthc,i,j,flag; PDL_Indx nthc,i,j,flag;
$SETNDIMS($PDL(PARENT)->ndims); $SETNDIMS($PDL(PARENT)->ndims);
$DOPRIVALLOC(); $DOPRIVALLOC();
$PRIV(offs) = 0; $PRIV(offs) = 0;
nthc=0; nthc=0;
for(i=0; i<$PDL(PARENT)->ndims; i++) { for(i=0; i<$PDL(PARENT)->ndims; i++) {
flag=0; flag=0;
if($PDL(PARENT)->nbroadcastids > $COMP(id) && $COMP(id) >= 0 && if($PDL(PARENT)->nbroadcastids > $COMP(id) && $COMP(id) >= 0 &&
i == $PDL(PARENT)->broadcastids[$COMP(id])) { i == $PDL(PARENT)->broadcastids[$COMP(id])) {
nthc += $COMP(whichdims_count); nthc += $COMP(whichdims_count);
} }
for(j=0; j<$COMP(whichdims_count); j++) { for(j=0; j<$COMP(whichdims_count); j++) {
if($COMP(whichdims[j] == i)) {flag=1; break;} if($COMP(whichdims)[j] == i) {flag=1; break;}
} }
if(flag) { if(flag) {
continue; continue;
} }
$PDL(CHILD)->dims[nthc] = $PDL(PARENT)->dims[i]; $PDL(CHILD)->dims[nthc] = $PDL(PARENT)->dims[i];
$PRIV(incs)[nthc] = $PDL(PARENT)->dimincs[i]; $PRIV(incs)[nthc] = $PDL(PARENT)->dimincs[i];
nthc++; nthc++;
} }
for(i=0; i<$COMP(whichdims_count); i++) { for(i=0; i<$COMP(whichdims_count); i++) {
PDL_Indx cdim,pdim; PDL_Indx cdim,pdim;
cdim = i + cdim = i +
($PDL(PARENT)->nbroadcastids > $COMP(id) && $COMP(id) > = 0? ($PDL(PARENT)->nbroadcastids > $COMP(id) && $COMP(id) > = 0?
$PDL(PARENT)->broadcastids[$COMP(id)] : $PDL(PARENT)-> ndims) $PDL(PARENT)->broadcastids[$COMP(id)] : $PDL(PARENT)-> ndims)
- $COMP(nrealwhichdims); - $COMP(nrealwhichdims);
pdim = $COMP(whichdims[i]); pdim = $COMP(whichdims)[i];
if(pdim == -1) { if(pdim == -1) {
$PDL(CHILD)->dims[cdim] = 1; $PDL(CHILD)->dims[cdim] = 1;
$PRIV(incs)[cdim] = 0; $PRIV(incs)[cdim] = 0;
} else { } else {
$PDL(CHILD)->dims[cdim] = $PDL(PARENT)->dims[pdi m]; $PDL(CHILD)->dims[cdim] = $PDL(PARENT)->dims[pdi m];
$PRIV(incs)[cdim] = $PDL(PARENT)->dimincs[pdim]; $PRIV(incs)[cdim] = $PDL(PARENT)->dimincs[pdim];
} }
} }
$SETDIMS(); $SETDIMS();
PDL_RETERROR(PDL_err, PDL->reallocbroadcastids($PDL(CHILD), PDL_RETERROR(PDL_err, PDL->reallocbroadcastids($PDL(CHILD),
PDLMAX($COMP(id)+1, $PDL(PARENT)->nbroadcastids) )); PDLMAX($COMP(id)+1, $PDL(PARENT)->nbroadcastids) ));
for(i=0; i<$PDL(CHILD)->nbroadcastids-1; i++) { for(i=0; i<$PDL(CHILD)->nbroadcastids-1; i++) {
$PDL(CHILD)->broadcastids[i] = $PDL(CHILD)->broadcastids[i] =
($PDL(PARENT)->nbroadcastids > i ? ($PDL(PARENT)->nbroadcastids > i ?
$PDL(PARENT)->broadcastids[i] : $PDL(PARENT)->ndims) + $PDL(PARENT)->broadcastids[i] : $PDL(PARENT)->ndims) +
(i <= $COMP(id) ? - $COMP(nrealwhichdims) : (i <= $COMP(id) ? - $COMP(nrealwhichdims) :
$COMP(whichdims_count) - $COMP(nrealwhichdims)); $COMP(whichdims_count) - $COMP(nrealwhichdims));
} }
$PDL(CHILD)->broadcastids[$PDL(CHILD)->nbroadcastids-1] = $PDL(C HILD)->ndims; $PDL(CHILD)->broadcastids[$PDL(CHILD)->nbroadcastids-1] = $PDL(C HILD)->ndims;
', EOF
); );
pp_def( pp_def(
'unbroadcast', 'unbroadcast',
Doc => <<'EOD', Doc => <<'EOD',
=for ref =for ref
All broadcasted dimensions are made real again. All broadcasted dimensions are made real again.
See [TBD Doc] for details and examples. See [TBD Doc] for details and examples.
=cut =cut
EOD EOD
P2Child => 1, P2Child => 1,
TwoWay => 1, TwoWay => 1,
AffinePriv => 1, AffinePriv => 1,
OtherPars => 'PDL_Indx atind;', OtherPars => 'PDL_Indx atind;',
RedoDims => ' RedoDims => pp_line_numbers(__LINE__, <<'EOF'),
PDL_Indx i; PDL_Indx i;
$SETNDIMS($PDL(PARENT)->ndims); $SETNDIMS($PDL(PARENT)->ndims);
$DOPRIVALLOC(); $DOPRIVALLOC();
$PRIV(offs) = 0; $PRIV(offs) = 0;
for(i=0; i<$PDL(PARENT)->ndims; i++) { for(i=0; i<$PDL(PARENT)->ndims; i++) {
PDL_Indx corc; PDL_Indx corc;
if(i<$COMP(atind)) { if(i<$COMP(atind)) {
corc = i; corc = i;
} else if(i < $PDL(PARENT)->broadcastids[0]) { } else if(i < $PDL(PARENT)->broadcastids[0]) {
corc = i + $PDL(PARENT)->ndims-$PDL(PARENT)->bro adcastids[0]; corc = i + $PDL(PARENT)->ndims-$PDL(PARENT)->bro adcastids[0];
} else { } else {
corc = i - $PDL(PARENT)->broadcastids[0] + $COMP (atind); corc = i - $PDL(PARENT)->broadcastids[0] + $COMP (atind);
} }
$PDL(CHILD)->dims[corc] = $PDL(PARENT)->dims[i]; $PDL(CHILD)->dims[corc] = $PDL(PARENT)->dims[i];
$PRIV(incs)[corc] = $PDL(PARENT)->dimincs[i]; $PRIV(incs)[corc] = $PDL(PARENT)->dimincs[i];
} }
$SETDIMS(); $SETDIMS();
', EOF
); );
pp_add_exported('', 'dice dice_axis'); pp_add_exported('', 'dice dice_axis');
pp_addpm(<<'EOD'); pp_addpm(<<'EOD');
=head2 dice =head2 dice
=for ref =for ref
Dice rows/columns/planes out of a PDL using indexes for Dice rows/columns/planes out of a PDL using indexes for
skipping to change at line 2421 skipping to change at line 2425
=for example =for example
$x->slice('1:3'); # return the second to fourth elements of $x $x->slice('1:3'); # return the second to fourth elements of $x
$x->slice('3:1'); # reverse the above $x->slice('3:1'); # reverse the above
$x->slice('-2:1'); # return last-but-one to second elements of $x $x->slice('-2:1'); # return last-but-one to second elements of $x
$x->slice([1,3]); # Same as above three calls, but using array ref syntax $x->slice([1,3]); # Same as above three calls, but using array ref syntax
$x->slice([3,1]); $x->slice([3,1]);
$x->slice([-2,1]); $x->slice([-2,1]);
EOD-slice EOD-slice
PMCode => <<'EOD-slice', PMCode => pp_line_numbers(__LINE__, <<'EOD-slice'),
sub PDL::slice { sub PDL::slice {
my ($source, @others) = @_; my ($source, @others) = @_;
for my $i(0..$#others) { for my $i(0..$#others) {
my $idx = $others[$i]; my $idx = $others[$i];
if (ref $idx eq 'ARRAY') { if (ref $idx eq 'ARRAY') {
my @arr = map UNIVERSAL::isa($_, 'PDL') ? $_->flat->at(0) : $_, @{$other s[$i]}; my @arr = map UNIVERSAL::isa($_, 'PDL') ? $_->flat->at(0) : $_, @{$other s[$i]};
$others[$i] = \@arr; $others[$i] = \@arr;
next; next;
} }
next if !( blessed($idx) && $idx->isa('PDL') ); next if !( blessed($idx) && $idx->isa('PDL') );
skipping to change at line 2483 skipping to change at line 2487
PDL_Indx odim[$COMP(nargs)]; PDL_Indx odim[$COMP(nargs)];
PDL_Indx idim[$COMP(nargs)]; PDL_Indx idim[$COMP(nargs)];
PDL_Indx idim_top; PDL_Indx idim_top;
PDL_Indx odim_top; PDL_Indx odim_top;
PDL_Indx start[$COMP(nargs)]; PDL_Indx start[$COMP(nargs)];
PDL_Indx inc[$COMP(nargs)]; PDL_Indx inc[$COMP(nargs)];
PDL_Indx end[$COMP(nargs)]; PDL_Indx end[$COMP(nargs)];
', ',
AffinePriv => 1, AffinePriv => 1,
TwoWay => 1, TwoWay => 1,
MakeComp => <<'SLICE-MC' MakeComp => pp_line_numbers(__LINE__, <<'SLICE-MC'),
PDL_Indx nargs = 0; PDL_Indx nargs = 0;
pdl_slice_args *argsptr = arglist; pdl_slice_args *argsptr = arglist;
while (argsptr) nargs++, argsptr = argsptr->next; while (argsptr) nargs++, argsptr = argsptr->next;
$COMP(nargs) = nargs; $COMP(nargs) = nargs;
$DOCOMPALLOC(); $DOCOMPALLOC();
PDL_Indx i, idim, odim, imax; PDL_Indx i, idim, odim, imax;
argsptr = arglist; argsptr = arglist;
for(odim=idim=i=0; i<nargs; i++) { for(odim=idim=i=0; i<nargs; i++) {
/* Copy parsed values into the limits */ /* Copy parsed values into the limits */
$COMP(start)[i] = argsptr->start; $COMP(start)[i] = argsptr->start;
$COMP(end)[i] = argsptr->end; $COMP(end)[i] = argsptr->end;
$COMP(inc)[i] = argsptr->inc; $COMP(inc)[i] = argsptr->inc;
/* Deal with dimensions */ /* Deal with dimensions */
$COMP(odim)[i] = argsptr->squish ? -1 : odim++; $COMP(odim)[i] = argsptr->squish ? -1 : odim++;
$COMP(idim)[i] = argsptr->dummy ? -1 : idim++; $COMP(idim)[i] = argsptr->dummy ? -1 : idim++;
argsptr = argsptr->next; argsptr = argsptr->next;
} /* end of arg-parsing loop */ } /* end of arg-parsing loop */
$COMP(idim_top) = idim; $COMP(idim_top) = idim;
$COMP(odim_top) = odim; $COMP(odim_top) = odim;
SLICE-MC SLICE-MC
, RedoDims => pp_line_numbers(__LINE__, <<'EOF'),
RedoDims => q{
PDL_Indx i; PDL_Indx i;
PDL_Indx PDIMS; PDL_Indx PDIMS;
PDL_Indx o_ndims_extra = PDLMAX(0, $PDL(PARENT)->ndims - $COMP (idim_top)); PDL_Indx o_ndims_extra = PDLMAX(0, $PDL(PARENT)->ndims - $COMP (idim_top));
/* slurped dims from the arg parsing, plus any extra broadcast dims */ /* slurped dims from the arg parsing, plus any extra broadcast dims */
$SETNDIMS( $COMP(odim_top) + o_ndims_extra ); $SETNDIMS( $COMP(odim_top) + o_ndims_extra );
$DOPRIVALLOC(); $DOPRIVALLOC();
$PRIV(offs) = 0; /* Offset vector to start of slice */ $PRIV(offs) = 0; /* Offset vector to start of slice */
for(i=0; i<$COMP(nargs); i++) { for(i=0; i<$COMP(nargs); i++) {
/** Belt-and-suspenders **/ /** Belt-and-suspenders **/
if( ($COMP(idim[i]) < 0) && ($COMP(odim[i]) < 0) ) { if( ($COMP(idim)[i] < 0) && ($COMP(odim)[i] < 0) ) {
PDL->changed($PDL(CHILD), PDL_PARENTDIMSCHANGED, 0); PDL->changed($PDL(CHILD), PDL_PARENTDIMSCHANGED, 0);
$CROAK("Hmmm, both dummy and squished -- this can never happen. I quit."); $CROAK("Hmmm, both dummy and squished -- this can never happen. I quit.");
} }
/** First handle dummy dims since there's no input from th e parent **/ /** First handle dummy dims since there's no input from th e parent **/
if( $COMP(idim[i]) < 0 ) { if( $COMP(idim)[i] < 0 ) {
/* dummy dim - offset or diminc. */ /* dummy dim - offset or diminc. */
$PDL(CHILD)->dims[ $COMP(odim[i]) ] = $COMP(end[i]) - $ $PDL(CHILD)->dims[ $COMP(odim)[i] ] = $COMP(end)[i] - $
COMP(start[i]) + 1; COMP(start)[i] + 1;
$PRIV(incs)[ $COMP(odim[i]) ] = 0; $PRIV(incs)[ $COMP(odim)[i] ] = 0;
} else { } else {
/** This is not a dummy dim -- deal with a regular slice along it. **/ /** This is not a dummy dim -- deal with a regular slice along it. **/
/** Get parent dim size for this idim, and/or allow perm issive slicing **/ /** Get parent dim size for this idim, and/or allow perm issive slicing **/
PDL_Indx pdsize = $COMP(idim[i]) < $PDL(PARENT)->ndims ? PDL_Indx pdsize = $COMP(idim)[i] < $PDL(PARENT)->ndims ?
$PDL(PARENT)->dims[$COMP(idim)[i]] : 1; $PDL(PARENT)->dims[$COMP(idim)[i]] : 1;
PDL_Indx start = $COMP(start)[i]; PDL_Indx start = $COMP(start)[i];
PDL_Indx end = $COMP(end)[i]; PDL_Indx end = $COMP(end)[i];
if( if(
/** Trap special case: full slices of an empty dim are empty **/ /** Trap special case: full slices of an empty dim are empty **/
(pdsize==0 && start==0 && end==-1 && $COMP(inc[i]) == 0) (pdsize==0 && start==0 && end==-1 && $COMP(inc)[i] == 0)
|| ||
/* the values given when PDL::slice gets empty ndarray for index */ /* the values given when PDL::slice gets empty ndarray for index */
(start==1 && end==0 && $COMP(inc[i]) == 1) (start==1 && end==0 && $COMP(inc)[i] == 1)
) { ) {
$PDL(CHILD)->dims[$COMP(odim)[i]] = 0; $PDL(CHILD)->dims[$COMP(odim)[i]] = 0;
$PRIV(incs)[$COMP(odim)[i]] = 0; $PRIV(incs)[$COMP(odim)[i]] = 0;
} else { } else {
/** Regularize and bounds-check the start location **/ /** Regularize and bounds-check the start location **/
if(start < 0) if(start < 0)
start += pdsize; start += pdsize;
if( start < 0 || start >= pdsize ) { if( start < 0 || start >= pdsize ) {
PDL->changed($PDL(CHILD), PDL_PARENTDIMSCHANGED, 0); PDL->changed($PDL(CHILD), PDL_PARENTDIMSCHANGED, 0);
if(i >= $PDL(PARENT)->ndims) { if(i >= $PDL(PARENT)->ndims) {
$CROAK("slice has too many dims (indexes dim %"IND_ FLAG"; highest is %"IND_FLAG")",i,$PDL(PARENT)->ndims-1); $CROAK("slice has too many dims (indexes dim %"IND_ FLAG"; highest is %"IND_FLAG")",i,$PDL(PARENT)->ndims-1);
} else { } else {
$CROAK("slice starts out of bounds in pos %"IND_FLAG" ( start is %"IND_FLAG"; source dim %"IND_FLAG" runs 0 to %"IND_FLAG")",i,start,$CO MP(idim[i]),pdsize-1); $CROAK("slice starts out of bounds in pos %"IND_FLAG" ( start is %"IND_FLAG"; source dim %"IND_FLAG" runs 0 to %"IND_FLAG")",i,start,$CO MP(idim)[i],pdsize-1);
} }
} }
if( $COMP(odim[i]) < 0) { if( $COMP(odim)[i] < 0) {
/* squished dim - just update the offset. */ /* squished dim - just update the offset. */
/* start is always defined and regularized if we are here here, since */ /* start is always defined and regularized if we are here here, since */
/* both idim[i] and odim[i] can't be <0 */ /* both idim[i] and odim[i] can't be <0 */
$PRIV(offs) += start * $PDL(PARENT)->dimincs[ $COMP(i dim[i]) ]; $PRIV(offs) += start * $PDL(PARENT)->dimincs[ $COMP(i dim)[i] ];
} else /* normal operation */ { } else /* normal operation */ {
/** Regularize and bounds-check the end location **/ /** Regularize and bounds-check the end location **/
if(end<0) end += pdsize; if(end<0) end += pdsize;
if( end < 0 || end >= pdsize ) { if( end < 0 || end >= pdsize ) {
PDL->changed($PDL(CHILD), PDL_PARENTDIMSCHANGED, 0) ; PDL->changed($PDL(CHILD), PDL_PARENTDIMSCHANGED, 0) ;
$CROAK("slice ends out of bounds in pos %"IND_FLAG" (end is %"IND_FLAG"; source dim %"IND_FLAG" runs 0 to %"IND_FLAG")",i,end,$COMP (idim[i]),pdsize-1); $CROAK("slice ends out of bounds in pos %"IND_FLAG" (end is %"IND_FLAG"; source dim %"IND_FLAG" runs 0 to %"IND_FLAG")",i,end,$COMP (idim)[i],pdsize-1);
} }
/* regularize inc */ /* regularize inc */
PDL_Indx inc = $COMP(inc)[i]; PDL_Indx inc = $COMP(inc)[i];
if(!inc) if(!inc)
inc = (start <= end) ? 1 : -1; inc = (start <= end) ? 1 : -1;
$PDL(CHILD)->dims[ $COMP(odim)[i] ] = PDLMAX(0, (end - start + inc) / inc); $PDL(CHILD)->dims[ $COMP(odim)[i] ] = PDLMAX(0, (end - start + inc) / inc);
$PRIV( incs )[ $COMP(odim)[i] ] = $PDL(PARENT)->dimi ncs[ $COMP(idim)[i] ] * inc; $PRIV( incs )[ $COMP(odim)[i] ] = $PDL(PARENT)->dimi ncs[ $COMP(idim)[i] ] * inc;
$PRIV(offs) += start * $PDL(PARENT)->dimincs[ $COMP(i dim)[i] ]; $PRIV(offs) += start * $PDL(PARENT)->dimincs[ $COMP(i dim)[i] ];
} /* end of normal slice case */ } /* end of normal slice case */
} /* end of trapped strange slice case */ } /* end of trapped strange slice case */
skipping to change at line 2587 skipping to change at line 2590
} /* end of nargs loop */ } /* end of nargs loop */
/* Now fill in broadcast dimensions as needed. idim and odim persist from the parsing loop */ /* Now fill in broadcast dimensions as needed. idim and odim persist from the parsing loop */
/* up above. */ /* up above. */
for(i=0; i<o_ndims_extra; i++) { for(i=0; i<o_ndims_extra; i++) {
$PDL(CHILD)->dims[ $COMP(odim_top) + i ] = $PDL(PARENT)->dim s[ $COMP(idim_top) + i ]; $PDL(CHILD)->dims[ $COMP(odim_top) + i ] = $PDL(PARENT)->dim s[ $COMP(idim_top) + i ];
$PRIV(incs)[ $COMP(odim_top) + i ] = $PDL(PARENT)->dimincs[ $COMP(idim_top) + i ]; $PRIV(incs)[ $COMP(odim_top) + i ] = $PDL(PARENT)->dimincs[ $COMP(idim_top) + i ];
} }
$SETDIMS(); $SETDIMS();
EOF
} # end of RedoDims for slice
); );
pp_addpm({At => 'Bot'},<< 'EOD'); pp_addpm({At => 'Bot'},<<'EOD');
=head1 BUGS =head1 BUGS
For the moment, you can't slice one of the zero-length dims of an For the moment, you can't slice one of the zero-length dims of an
empty ndarray. It is not clear how to implement this in a way that makes empty ndarray. It is not clear how to implement this in a way that makes
sense. sense.
Many types of index errors are reported far from the indexing Many types of index errors are reported far from the indexing
operation that caused them. This is caused by the underlying architecture: operation that caused them. This is caused by the underlying architecture:
slice() sets up a mapping between variables, but that mapping isn't slice() sets up a mapping between variables, but that mapping isn't
 End of changes. 103 change blocks. 
127 lines changed or deleted 129 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)