"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

---