"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Basic/Primitive/primitive.pd" between
PDL-2.075.tar.gz and PDL-2.076.tar.gz

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

primitive.pd  (PDL-2.075):primitive.pd  (PDL-2.076)
skipping to change at line 140 skipping to change at line 140
PDL overloads the C<x> operator (normally the repeat operator) for PDL overloads the C<x> operator (normally the repeat operator) for
matrix multiplication. The number of columns (size of the 0 matrix multiplication. The number of columns (size of the 0
dimension) in the left-hand argument must normally equal the number of dimension) in the left-hand argument must normally equal the number of
rows (size of the 1 dimension) in the right-hand argument. rows (size of the 1 dimension) in the right-hand argument.
Row vectors are represented as (N x 1) two-dimensional PDLs, or you Row vectors are represented as (N x 1) two-dimensional PDLs, or you
may be sloppy and use a one-dimensional PDL. Column vectors are may be sloppy and use a one-dimensional PDL. Column vectors are
represented as (1 x N) two-dimensional PDLs. represented as (1 x N) two-dimensional PDLs.
Threading occurs in the usual way, but as both the 0 and 1 dimension Broadcasting occurs in the usual way, but as both the 0 and 1 dimension
(if present) are included in the operation, you must be sure that (if present) are included in the operation, you must be sure that
you don't try to broadcast over either of those dims. you don't try to broadcast over either of those dims.
Of note, due to how Perl v5.14.0 and above implement operator overloading of Of note, due to how Perl v5.14.0 and above implement operator overloading of
the C<x> operator, the use of parentheses for the left operand creates a list the C<x> operator, the use of parentheses for the left operand creates a list
context, that is context, that is
pdl> ( $x * $y ) x $z pdl> ( $x * $y ) x $z
ERROR: Argument "..." isn't numeric in repeat (x) ... ERROR: Argument "..." isn't numeric in repeat (x) ...
skipping to change at line 213 skipping to change at line 213
INTERNALS INTERNALS
The mechanics of the multiplication are carried out by the The mechanics of the multiplication are carried out by the
L</matmult> method. L</matmult> method.
=cut =cut
EOD EOD
pp_add_exported('', 'matmult');
pp_def('matmult', pp_def('matmult',
HandleBad=>0, HandleBad=>0,
Pars => 'a(t,h); b(w,t); [o]c(w,h);', Pars => 'a(t,h); b(w,t); [o]c(w,h);',
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
PMCode => <<'EOPM', PMCode => <<'EOPM',
sub PDL::matmult { sub PDL::matmult {
my ($x,$y,$c) = @_; my ($x,$y,$c) = @_;
$y = PDL->topdl($y);
$y = pdl($y) unless eval { $y->isa('PDL') }; $c = PDL->null unless do { local $@; eval { $c->isa('PDL') } };
$c = PDL->null unless eval { $c->isa('PDL') };
while($x->getndims < 2) {$x = $x->dummy(-1)} while($x->getndims < 2) {$x = $x->dummy(-1)}
while($y->getndims < 2) {$y = $y->dummy(-1)} while($y->getndims < 2) {$y = $y->dummy(-1)}
return ($c .= $x * $y) if( ($x->dim(0)==1 && $x->dim(1)==1) || return ($c .= $x * $y) if( ($x->dim(0)==1 && $x->dim(1)==1) ||
($y->dim(0)==1 && $y->dim(1)==1) ); ($y->dim(0)==1 && $y->dim(1)==1) );
if($y->dim(1) != $x->dim(0)) { barf sprintf 'Dim mismatch in matmult of [%1$dx%2$d] x [%3$dx%4$d]: %1$d !=
barf(sprintf("Dim mismatch in matmult of [%dx%d] x [%dx%d]: %d != %d",$x %4$d',$x->dim(0),$x->dim(1),$y->dim(0),$y->dim(1)
->dim(0),$x->dim(1),$y->dim(0),$y->dim(1),$x->dim(0),$y->dim(1))); if $y->dim(1) != $x->dim(0);
}
PDL::_matmult_int($x,$y,$c); PDL::_matmult_int($x,$y,$c);
$c; $c;
} }
EOPM EOPM
Code => <<'EOC', Code => <<'EOC',
PDL_Indx ih, iw, it, ow, oh, ot, wlim, hlim, tlim; PDL_Indx ih, iw, it, ow, oh, ot, wlim, hlim, tlim;
$GENERIC() *ad, *bd; PDL_Indx tsiz = 8 * sizeof(double) / sizeof($GENERIC());
PDL_Indx atdi, btdi;
PDL_Indx tsiz = 64;
// Zero the output
loop(w) %{
loop(h) %{
$c() = 0;
%}
%}
// Make sure we're physical // Cache the dimincs to avoid constant lookups
// (Not needed if we don't need dimincs, see below) PDL_Indx atdi = PDL_REPRINCS($PDL(a))[0];
// PDL->make_physdims($PDL(a)); PDL_Indx btdi = PDL_REPRINCS($PDL(b))[1];
// PDL->make_physdims($PDL(b));
// Cache the dimincs to avoid constant lookups
// These two lines are what I wanted, but they break sometimes (dimincs n
ot set right despite calling physdims?)
// I deleted them in favor of explicit offset calculation, which appears
more robust.
// atdi = $PDL(a)->dimincs[0];
// btdi = $PDL(b)->dimincs[1];
atdi = &($a(t=>1, h=>0)) - &($a(t=>0,h=>0));
btdi = &($b(t=>1, w=>0)) - &($b(t=>0,w=>0));
// Loop over tiles // Loop over tiles
for( oh=0; oh < $SIZE(h); oh += tsiz ) { for( oh=0; oh < $SIZE(h); oh += tsiz ) {
hlim = PDLMIN(oh + tsiz, $SIZE(h)); hlim = PDLMIN(oh + tsiz, $SIZE(h));
for( ow=0; ow < $SIZE(w); ow += tsiz ) { for( ow=0; ow < $SIZE(w); ow += tsiz ) {
wlim = PDLMIN(ow + tsiz, $SIZE(w)); wlim = PDLMIN(ow + tsiz, $SIZE(w));
// Zero the output for this tile
for( ih=oh; ih<hlim; ih++ )
for( iw=ow; iw<wlim; iw++ )
$c(w=>iw, h=>ih) = 0;
for( ot=0; ot < $SIZE(t); ot += tsiz ) { for( ot=0; ot < $SIZE(t); ot += tsiz ) {
tlim = PDLMIN(ot + tsiz, $SIZE(t)); tlim = PDLMIN(ot + tsiz, $SIZE(t));
for( ih=oh; ih<hlim; ih++ ) { for( ih=oh; ih<hlim; ih++ ) {
for( iw=ow; iw<wlim; iw++ ) { for( iw=ow; iw<wlim; iw++ ) {
$GENERIC() cc; // Cache the accumulated value for the output
$GENERIC() cc = $c(w=>iw, h=>ih);
// Cache data pointers before 't' run through tile // Cache data pointers before 't' run through tile
ad = &($a(t=>ot, h=>ih)); $GENERIC() *ad = &($a(t=>ot, h=>ih));
bd = &($b(w=>iw, t=>ot)); $GENERIC() *bd = &($b(w=>iw, t=>ot));
// Cache the accumulated value for the output
cc = $c(w=>iw, h=>ih);
// Hotspot - run the 't' summation // Hotspot - run the 't' summation
for( it=ot; it<tlim; it++ ) { for( it=ot; it<tlim; it++ ) {
cc += *ad * *bd; cc += *ad * *bd;
ad += atdi; ad += atdi;
bd += btdi; bd += btdi;
} }
// put the output back to be further accumulated later // put the output back to be further accumulated later
$c(w=>iw, h=>ih) = cc; $c(w=>iw, h=>ih) = cc;
} }
} }
} }
} }
skipping to change at line 1021 skipping to change at line 1000
} elsif (\$#_ > 1) {\$c=\$_[2]} else {\$c=PDL->nullcreate(\$x)} } elsif (\$#_ > 1) {\$c=\$_[2]} else {\$c=PDL->nullcreate(\$x)}
PDL::_${name}_int(\$x,\$y,\$c); PDL::_${name}_int(\$x,\$y,\$c);
return \$c; return \$c;
} }
EOD EOD
); # pp_def $name ); # pp_def $name
} # for: my $opt } # for: my $opt
pp_add_exported('', 'clip'); pp_add_exported('', 'clip');
pp_addpm(<<'EOD'); pp_addpm(<<'EOD');
=head2 clip =head2 clip
=for ref =for ref
Clip (threshold) an ndarray by (optional) upper or lower bounds. Clip (threshold) an ndarray by (optional) upper or lower bounds.
=for usage =for usage
$y = $x->clip(0,3); $y = $x->clip(0,3);
skipping to change at line 1910 skipping to change at line 1888
sub randsym { ref($_[0]) && ref($_[0]) ne 'PDL::Type' ? $_[0]->randsym : PDL->ra ndsym(@_) } sub randsym { ref($_[0]) && ref($_[0]) ne 'PDL::Type' ? $_[0]->randsym : PDL->ra ndsym(@_) }
sub PDL::randsym { sub PDL::randsym {
my $class = shift; my $class = shift;
my $x = scalar(@_)? $class->new_from_specification(@_) : $class->new_or_inpla ce; my $x = scalar(@_)? $class->new_from_specification(@_) : $class->new_or_inpla ce;
PDL::_randsym_int($x); PDL::_randsym_int($x);
return $x; return $x;
} }
EOD EOD
); );
pp_add_exported('','grandom');
pp_addpm(<<'EOD'); pp_addpm(<<'EOD');
=head2 grandom =head2 grandom
=for ref =for ref
Constructor which returns ndarray of Gaussian random numbers Constructor which returns ndarray of Gaussian random numbers
=for usage =for usage
$x = grandom([type], $nx, $ny, $nz,...); $x = grandom([type], $nx, $ny, $nz,...);
skipping to change at line 1941 skipping to change at line 1920
=cut =cut
sub grandom { ref($_[0]) && ref($_[0]) ne 'PDL::Type' ? $_[0]->grandom : PDL->gr andom(@_) } sub grandom { ref($_[0]) && ref($_[0]) ne 'PDL::Type' ? $_[0]->grandom : PDL->gr andom(@_) }
sub PDL::grandom { sub PDL::grandom {
my $class = shift; my $class = shift;
my $x = scalar(@_)? $class->new_from_specification(@_) : $class->new_or_inpla ce; my $x = scalar(@_)? $class->new_from_specification(@_) : $class->new_or_inpla ce;
use PDL::Math 'ndtri'; use PDL::Math 'ndtri';
$x .= ndtri(randsym($x)); $x .= ndtri(randsym($x));
return $x; return $x;
} }
EOD EOD
pp_add_exported('','grandom');
############################################################### ###############################################################
# binary searches in an ndarray; various forms # binary searches in an ndarray; various forms
############################################################### ###############################################################
pp_add_exported('','vsearch');
# generic front end; defaults to vsearch_sample for backwards compatibility # generic front end; defaults to vsearch_sample for backwards compatibility
pp_add_exported('','vsearch');
pp_addpm(<<'EOD'); pp_addpm(<<'EOD');
=head2 vsearch =head2 vsearch
=for sig =for sig
Signature: ( vals(); xs(n); [o] indx(); [\%options] ) Signature: ( vals(); xs(n); [o] indx(); [\%options] )
=for ref =for ref
Efficiently search for values in a sorted ndarray, returning indices. Efficiently search for values in a sorted ndarray, returning indices.
skipping to change at line 2873 skipping to change at line 2848
); );
# "Collapse" the subpixel offset... # "Collapse" the subpixel offset...
$y = $y->slice(":,($i)"); $y = $y->slice(":,($i)");
} }
return $samp; return $samp;
} elsif($method =~ m/^f(ft|ourier)?/i) { } elsif($method =~ m/^f(ft|ourier)?/i) {
local $@;
eval "use PDL::FFT;"; eval "use PDL::FFT;";
my $fftref = $opt->{fft}; my $fftref = $opt->{fft};
$fftref = [] unless(ref $fftref eq 'ARRAY'); $fftref = [] unless(ref $fftref eq 'ARRAY');
if(@$fftref != 2) { if(@$fftref != 2) {
my $x = $source->copy; my $x = $source->copy;
my $y = zeroes($source); my $y = zeroes($source);
fftnd($x,$y); fftnd($x,$y);
$fftref->[0] = sqrt($x*$x+$y*$y) / $x->nelem; $fftref->[0] = sqrt($x*$x+$y*$y) / $x->nelem;
$fftref->[1] = - atan2($y,$x); $fftref->[1] = - atan2($y,$x);
} }
skipping to change at line 3120 skipping to change at line 3096
else else
for (i=0; i<dpdl->dims[0]; i++) { for (i=0; i<dpdl->dims[0]; i++) {
$GENERIC() foo = *(m_datap+inc*i+offs); $GENERIC() foo = *(m_datap+inc*i+offs);
if(foo) sum++; if(foo) sum++;
} }
'. $_->{Autosize} '. $_->{Autosize}
); );
} }
pp_add_exported("", 'where');
pp_addpm(<<'EOD'); pp_addpm(<<'EOD');
=head2 where =head2 where
=for ref =for ref
Use a mask to select values from one or more data PDLs Use a mask to select values from one or more data PDLs
C<where> accepts one or more data ndarrays and a mask ndarray. It C<where> accepts one or more data ndarrays and a mask ndarray. It
returns a list of output ndarrays, corresponding to the input data returns a list of output ndarrays, corresponding to the input data
ndarrays. Each output ndarray is a 1-dimensional list of values in its ndarrays. Each output ndarray is a 1-dimensional list of values in its
skipping to change at line 3188 skipping to change at line 3165
} else { } else {
my $mask = $_[-1]->which; my $mask = $_[-1]->which;
return map {$_->index($mask)} @_[0..$#_-1]; return map {$_->index($mask)} @_[0..$#_-1];
} }
} }
} }
*where = \&PDL::where; *where = \&PDL::where;
EOD EOD
pp_add_exported("", 'where'); pp_add_exported("", 'whereND');
pp_addpm(<<'EOD'); pp_addpm(<<'EOD');
=head2 whereND =head2 whereND
=for ref =for ref
C<where> with support for ND masks and broadcasting C<where> with support for ND masks and broadcasting
C<whereND> accepts one or more data ndarrays and a C<whereND> accepts one or more data ndarrays and a
mask ndarray. It returns a list of output ndarrays, mask ndarray. It returns a list of output ndarrays,
corresponding to the input data ndarrays. The values corresponding to the input data ndarrays. The values
skipping to change at line 3279 skipping to change at line 3255
push @to_return, $where_sub_i; push @to_return, $where_sub_i;
} }
return (@to_return == 1) ? $to_return[0] : @to_return; return (@to_return == 1) ? $to_return[0] : @to_return;
} }
*whereND = \&PDL::whereND; *whereND = \&PDL::whereND;
EOD EOD
pp_add_exported("", 'whereND'); pp_add_exported("", 'whichND');
pp_addpm(<<'EOD'); pp_addpm(<<'EOD');
=head2 whichND =head2 whichND
=for ref =for ref
Return the coordinates of non-zero values in a mask. Return the coordinates of non-zero values in a mask.
=for usage =for usage
WhichND returns the N-dimensional coordinates of each nonzero value in WhichND returns the N-dimensional coordinates of each nonzero value in
skipping to change at line 3383 skipping to change at line 3358
} }
for my $i (0..$#mdims) { for my $i (0..$#mdims) {
my($s) = $ind->index($i); my($s) = $ind->index($i);
$s /= $mult->index($i); $s /= $mult->index($i);
$s %= $mdims[$i]; $s %= $mdims[$i];
} }
return $ind; return $ind;
} }
EOD EOD
pp_add_exported("", 'whichND');
# #
# Set operations suited for manipulation of the operations above. # Set operations suited for manipulation of the operations above.
# #
pp_add_exported("", 'setops');
pp_addpm(<<'EOD'); pp_addpm(<<'EOD');
=head2 setops =head2 setops
=for ref =for ref
Implements simple set operations like union and intersection Implements simple set operations like union and intersection
=for usage =for usage
Usage: $set = setops($x, <OPERATOR>, $y); Usage: $set = setops($x, <OPERATOR>, $y);
skipping to change at line 3573 skipping to change at line 3546
return $union->where($union == rotate($union, -1)); return $union->where($union == rotate($union, -1));
} else { } else {
print "The operation $op is not known!"; print "The operation $op is not known!";
return -1; return -1;
} }
} }
EOD EOD
pp_add_exported("", 'setops'); pp_add_exported("", 'intersect');
pp_addpm(<<'EOD'); pp_addpm(<<'EOD');
=head2 intersect =head2 intersect
=for ref =for ref
Calculate the intersection of two ndarrays Calculate the intersection of two ndarrays
=for usage =for usage
Usage: $set = intersect($x, $y); Usage: $set = intersect($x, $y);
skipping to change at line 3609 skipping to change at line 3581
=cut =cut
*intersect = \&PDL::intersect; *intersect = \&PDL::intersect;
sub PDL::intersect { sub PDL::intersect {
return setops($_[0], 'AND', $_[1]); return setops($_[0], 'AND', $_[1]);
} }
EOD EOD
pp_add_exported("", 'intersect');
pp_addpm({At=>'Bot'},<<'EOD'); pp_addpm({At=>'Bot'},<<'EOD');
=head1 AUTHOR =head1 AUTHOR
Copyright (C) Tuomas J. Lukka 1997 (lukka@husc.harvard.edu). Contributions Copyright (C) Tuomas J. Lukka 1997 (lukka@husc.harvard.edu). Contributions
by Christian Soeller (c.soeller@auckland.ac.nz), Karl Glazebrook by Christian Soeller (c.soeller@auckland.ac.nz), Karl Glazebrook
(kgb@aaoepp.aao.gov.au), Craig DeForest (deforest@boulder.swri.edu) (kgb@aaoepp.aao.gov.au), Craig DeForest (deforest@boulder.swri.edu)
and Jarle Brinchmann (jarle@astro.up.pt) and Jarle Brinchmann (jarle@astro.up.pt)
All rights reserved. There is no warranty. You are allowed All rights reserved. There is no warranty. You are allowed
to redistribute this software / documentation under certain to redistribute this software / documentation under certain
 End of changes. 27 change blocks. 
61 lines changed or deleted 28 lines changed or added

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