"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Basic/Primitive/primitive.pd" between
PDL-2.077.tar.gz and PDL-2.078.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.077):primitive.pd  (PDL-2.078)
skipping to change at line 58 skipping to change at line 58
# available via operator overloading) # available via operator overloading)
################################################################ ################################################################
pp_def( pp_def(
'inner', 'inner',
HandleBad => 1, HandleBad => 1,
Pars => 'a(n); b(n); [o]c();', Pars => 'a(n); b(n); [o]c();',
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
Code => Code =>
'double tmp = 0; 'double tmp = 0;
loop(n) %{ tmp += $a() * $b(); %}
$c() = tmp;',
BadCode =>
'double tmp = 0;
int badflag = 0; int badflag = 0;
loop(n) %{ loop(n) %{
if ( $ISGOOD(a()) && $ISGOOD(b()) ) { tmp += $a() * $b(); } else { ba PDL_IF_BAD(if (!($ISGOOD(a()) && $ISGOOD(b()))) { badflag = 1; break;
dflag = 1; } }
else,) { tmp += $a() * $b();
}
%} %}
if ( badflag ) { $SETBAD(c()); $PDLSTATESETBAD(c); } PDL_IF_BAD(if (badflag) { $SETBAD(c()); $PDLSTATESETBAD(c); }
else { $c() = tmp; }', else,) { $c() = tmp; }',
CopyBadStatusCode => '', CopyBadStatusCode => '',
Doc => ' Doc => '
=for ref =for ref
Inner product over one dimension Inner product over one dimension
c = sum_i a_i * b_i c = sum_i a_i * b_i
=cut =cut
skipping to change at line 100 skipping to change at line 97
', ',
); # pp_def( inner ) ); # pp_def( inner )
pp_def( pp_def(
'outer', 'outer',
HandleBad => 1, HandleBad => 1,
Pars => 'a(n); b(m); [o]c(n,m);', Pars => 'a(n); b(m); [o]c(n,m);',
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
Code => Code =>
'loop(n,m) %{ 'loop(n,m) %{
$c() = $a() * $b(); PDL_IF_BAD(if ( $ISBAD(a()) || $ISBAD(b()) ) {
%}',
BadCode =>
'loop(n,m) %{
if ( $ISBAD(a()) || $ISBAD(b()) ) {
$SETBAD(c()); $SETBAD(c());
} else { } else,) {
$c() = $a() * $b(); $c() = $a() * $b();
} }
%}', %}',
Doc => ' Doc => '
=for ref =for ref
outer product over one dimension outer product over one dimension
Naturally, it is possible to achieve the effects of outer Naturally, it is possible to achieve the effects of outer
product simply by broadcasting over the "C<*>" product simply by broadcasting over the "C<*>"
skipping to change at line 304 skipping to change at line 297
EOD EOD
); );
pp_def( pp_def(
'innerwt', 'innerwt',
HandleBad => 1, HandleBad => 1,
Pars => 'a(n); b(n); c(n); [o]d();', Pars => 'a(n); b(n); c(n); [o]d();',
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
Code => Code =>
'double tmp = 0; 'long double tmp = 0;
loop(n) %{
tmp += $a() * $b() * $c();
%}
$d() = tmp;',
BadCode =>
'double tmp = 0;
int flag = 0; int flag = 0;
loop(n) %{ loop(n) %{
if ( $ISGOOD(a()) && $ISGOOD(b()) && $ISGOOD(c()) ) { PDL_IF_BAD(if ( $ISGOOD(a()) && $ISGOOD(b()) && $ISGOOD(c()) ),) {
tmp += $a() * $b() * $c(); tmp += $a() * $b() * $c();
flag = 1; flag = 1;
} }
%} %}
if ( flag ) { $d() = tmp; } PDL_IF_BAD(if (!flag) { $SETBAD(d()); }
else { $SETBAD(d()); }', else,) { $d() = tmp; }',
Doc => ' Doc => '
=for ref =for ref
Weighted (i.e. triple) inner product Weighted (i.e. triple) inner product
d = sum_i a(i) b(i) c(i) d = sum_i a(i) b(i) c(i)
=cut =cut
' '
); );
pp_def( pp_def(
'inner2', 'inner2',
HandleBad => 1, HandleBad => 1,
Pars => 'a(n); b(n,m); c(m); [o]d();', Pars => 'a(n); b(n,m); c(m); [o]d();',
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
Code => Code =>
'double tmp=0; 'long double tmp = 0;
loop(n,m) %{
tmp += $a() * $b() * $c();
%}
$d() = tmp;',
BadCode =>
'double tmp = 0;
int flag = 0; int flag = 0;
loop(n,m) %{ loop(n,m) %{
if ( $ISGOOD(a()) && $ISGOOD(b()) && $ISGOOD(c()) ) { PDL_IF_BAD(if ( $ISGOOD(a()) && $ISGOOD(b()) && $ISGOOD(c()) ),) {
tmp += $a() * $b() * $c(); tmp += $a() * $b() * $c();
flag = 1; flag = 1;
} }
%} %}
if ( flag ) { $d() = tmp; } PDL_IF_BAD(if (!flag) { $SETBAD(d()); }
else { $SETBAD(d()); }', else,) { $d() = tmp; }',
Doc => ' Doc => '
=for ref =for ref
Inner product of two vectors and a matrix Inner product of two vectors and a matrix
d = sum_ij a(i) b(i,j) c(j) d = sum_ij a(i) b(i,j) c(j)
Note that you should probably not broadcast over C<a> and C<c> since that would be Note that you should probably not broadcast over C<a> and C<c> since that would be
very wasteful. Instead, you should use a temporary for C<b*c>. very wasteful. Instead, you should use a temporary for C<b*c>.
' '
); );
pp_def( pp_def(
'inner2d', 'inner2d',
HandleBad => 1, HandleBad => 1,
Pars => 'a(n,m); b(n,m); [o]c();', Pars => 'a(n,m); b(n,m); [o]c();',
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
Code => Code =>
'double tmp=0; 'long double tmp = 0;
loop(n,m) %{
tmp += $a() * $b();
%}
$c() = tmp;',
BadCode =>
'double tmp = 0;
int flag = 0; int flag = 0;
loop(n,m) %{ loop(n,m) %{
if ( $ISGOOD(a()) && $ISGOOD(b()) ) { PDL_IF_BAD(if ( $ISGOOD(a()) && $ISGOOD(b()) ),) {
tmp += $a() * $b(); tmp += $a() * $b();
flag = 1; flag = 1;
} }
%} %}
if ( flag ) { $c() = tmp; } PDL_IF_BAD(if (!flag) { $SETBAD(c()); }
else { $SETBAD(c()); }', else,) { $c() = tmp; }',
Doc => ' Doc => '
=for ref =for ref
Inner product over 2 dimensions. Inner product over 2 dimensions.
Equivalent to Equivalent to
$c = inner($x->clump(2), $y->clump(2)) $c = inner($x->clump(2), $y->clump(2))
skipping to change at line 412 skipping to change at line 386
' '
); );
pp_def( pp_def(
'inner2t', 'inner2t',
HandleBad => 1, HandleBad => 1,
Pars => 'a(j,n); b(n,m); c(m,k); [t]tmp(n,k); [o]d(j,k));', Pars => 'a(j,n); b(n,m); c(m,k); [t]tmp(n,k); [o]d(j,k));',
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
Code => Code =>
'loop(n,k) %{ 'loop(n,k) %{
double tmp0 = 0; long double tmp0 = 0;
loop(m) %{
tmp0 += $b() * $c();
%}
$tmp() = tmp0;
%}
loop(j,k) %{
double tmp1 = 0;
loop(n) %{
tmp1 += $a() * $tmp();
%}
$d() = tmp1;
%}',
BadCode =>
'loop(n,k) %{
double tmp0 = 0;
int flag = 0; int flag = 0;
loop(m) %{ loop(m) %{
if ( $ISGOOD(b()) && $ISGOOD(c()) ) { PDL_IF_BAD(if ( $ISGOOD(b()) && $ISGOOD(c()) ),) {
tmp0 += $b() * $c(); tmp0 += $b() * $c();
flag = 1; flag = 1;
} }
%} %}
if ( flag ) { $tmp() = tmp0; } PDL_IF_BAD(if (!flag) { $SETBAD(tmp()); }
else { $SETBAD(tmp()); } else,) { $tmp() = tmp0; }
%} %}
loop(j,k) %{ loop(j,k) %{
double tmp1 = 0; long double tmp1 = 0;
int flag = 0; int flag = 0;
loop(n) %{ loop(n) %{
if ( $ISGOOD(a()) && $ISGOOD(tmp()) ) { if ( $ISGOOD(a()) && $ISGOOD(tmp()) ) {
tmp1 += $a() * $tmp(); tmp1 += $a() * $tmp();
flag = 1; flag = 1;
} }
%} %}
if ( flag ) { $d() = tmp1; } PDL_IF_BAD(if (!flag) { $SETBAD(d()); }
else { $SETBAD(d()); } else,) { $d() = tmp1; }
%}', %}',
Doc => ' Doc => '
=for ref =for ref
Efficient Triple matrix product C<a*b*c> Efficient Triple matrix product C<a*b*c>
Efficiency comes from by using the temporary C<tmp>. This operation only Efficiency comes from by using the temporary C<tmp>. This operation only
scales as C<N**3> whereas broadcasting using L</inner2> would scale scales as C<N**3> whereas broadcasting using L</inner2> would scale
as C<N**4>. as C<N**4>.
The reason for having this routine is that you do not need to The reason for having this routine is that you do not need to
have the same broadcast-dimensions for C<tmp> as for the other arguments, have the same broadcast-dimensions for C<tmp> as for the other arguments,
which in case of large numbers of matrices makes this much more which in case of large numbers of matrices makes this much more
memory-efficient. memory-efficient.
It is hoped that things like this could be taken care of as a kind of It is hoped that things like this could be taken care of as a kind of
closures at some point. closures at some point.
=cut =cut
' '
); # pp_def inner2t() ); # pp_def inner2t()
# a helper function for the cross product definition # a helper function for the cross product definition
sub crassgn { sub crassgn {
"\$c(tri => $_[0]) = \$a(tri => $_[1])*\$b(tri => $_[2]) - "\$c(tri => $_[0]) = \$a(tri => $_[1])*\$b(tri => $_[2]) -
\$a(tri => $_[2])*\$b(tri => $_[1]);" \$a(tri => $_[2])*\$b(tri => $_[1]);"
} }
pp_def('crossp', pp_def('crossp',
skipping to change at line 556 skipping to change at line 513
# this one was motivated by the need to compute # this one was motivated by the need to compute
# the circular mean efficiently # the circular mean efficiently
# without it could not be done efficiently or without # without it could not be done efficiently or without
# creating large intermediates (check pdl-porters for # creating large intermediates (check pdl-porters for
# discussion) # discussion)
# see PDL::ImageND for info about the circ_mean function # see PDL::ImageND for info about the circ_mean function
pp_def( pp_def(
'indadd', 'indadd',
HandleBad => 1, HandleBad => 1,
Pars => 'a(n); indx ind(n); [o] sum(m)', Pars => 'input(n); indx ind(n); [io] sum(m)',
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
Code => Code =>
'loop(n) %{ 'loop(n) %{
register PDL_Indx foo = $ind(); register PDL_Indx this_ind = $ind();
if ( foo<0 || foo>=$SIZE(m) ) { if ( PDL_IF_BAD($ISBADVAR(this_ind,ind) ||,) this_ind<0 || this_ind>=$SI
$CROAK("PDL::indadd: invalid index"); ZE(m) )
} $CROAK("invalid index");
$sum(m => foo) += $a(); PDL_IF_BAD(
if ( $ISBAD(input()) ) { $SETBAD(sum(m => this_ind)); }
else,) { $sum(m => this_ind) += $input(); }
%}', %}',
BadCode =>
'loop(n) %{
register PDL_Indx foo = $ind();
if( $ISBADVAR(foo,ind) || foo<0 || foo>=$SIZE(m) ) {
$CROAK("PDL::indadd: invalid index");
}
if ( $ISBAD(a()) ) { $SETBAD(sum(m => foo)); }
else { $sum(m => foo) += $a(); }
%}',
BadDoc => ' BadDoc => '
=for bad =for bad
The routine barfs if any of the indices are bad. The routine barfs on bad indices, and bad inputs set target outputs bad.
=cut =cut
', ',
Doc=>' Doc=>'
=for ref =for ref
Threaded Index Add: Add C<a> to the C<ind> element of C<sum>, i.e: Broadcasting index add: add C<input> to the C<ind> element of C<sum>, i.e:
sum(ind) += a sum(ind) += input
=for example =for example
Simple Example: Simple example:
$x = 2; $x = 2;
$ind = 3; $ind = 3;
$sum = zeroes(10); $sum = zeroes(10);
indadd($x,$ind, $sum); indadd($x,$ind, $sum);
print $sum print $sum
#Result: ( 2 added to element 3 of $sum) #Result: ( 2 added to element 3 of $sum)
# [0 0 0 2 0 0 0 0 0 0] # [0 0 0 2 0 0 0 0 0 0]
Threaded Example: Broadcasting example:
$x = pdl( 1,2,3); $x = pdl( 1,2,3);
$ind = pdl( 1,4,6); $ind = pdl( 1,4,6);
$sum = zeroes(10); $sum = zeroes(10);
indadd($x,$ind, $sum); indadd($x,$ind, $sum);
print $sum."\n"; print $sum."\n";
#Result: ( 1, 2, and 3 added to elements 1,4,6 $sum) #Result: ( 1, 2, and 3 added to elements 1,4,6 $sum)
# [0 1 0 0 2 0 3 0 0 0] # [0 1 0 0 2 0 3 0 0 0]
=cut =cut
'); ');
# 1D convolution # 1D convolution
# useful for broadcasted 1D filters # useful for broadcasted 1D filters
pp_addhdr(' pp_addhdr('
/* Fast Modulus with proper negative behaviour */ /* Fast Modulus with proper negative behaviour */
#define REALMOD(a,b) while ((a)>=(b)) (a) -= (b); while ((a)<0) (a) += (b); #define REALMOD(a,b) while ((a)>=(b)) (a) -= (b); while ((a)<0) (a) += (b);
'); ');
pp_def('conv1d', pp_def('conv1d',
Doc => << 'EOD', Doc => << 'EOD',
skipping to change at line 972 skipping to change at line 917
['hclip','PDLMIN'], ['hclip','PDLMIN'],
['lclip','PDLMAX'] ['lclip','PDLMAX']
) { ) {
my $name = $opt->[0]; my $name = $opt->[0];
my $op = $opt->[1]; my $op = $opt->[1];
my $code = '$c() = '.$op.'($b(), $a());'; my $code = '$c() = '.$op.'($b(), $a());';
pp_def( pp_def(
$name, $name,
HandleBad => 1, HandleBad => 1,
Pars => 'a(); b(); [o] c()', Pars => 'a(); b(); [o] c()',
Code => $code, Code =>
BadCode => 'PDL_IF_BAD(if ( $ISBAD(a()) || $ISBAD(b()) ) {
'if ( $ISBAD(a()) || $ISBAD(b()) ) {
$SETBAD(c()); $SETBAD(c());
} else { '.$code.' }', } else,) { '.$code.' }',
Doc => 'clip (threshold) C<$a> by C<$b> (C<$b> is '. Doc => 'clip (threshold) C<$a> by C<$b> (C<$b> is '.
($name eq 'hclip' ? 'upper' : 'lower').' bound)', ($name eq 'hclip' ? 'upper' : 'lower').' bound)',
PMCode=><<"EOD", PMCode=><<"EOD",
sub PDL::$name { sub PDL::$name {
my (\$x,\$y) = \@_; my (\$x,\$y) = \@_;
my \$c; my \$c;
if (\$x->is_inplace) { if (\$x->is_inplace) {
\$x->set_inplace(0); \$c = \$x; \$x->set_inplace(0); \$c = \$x;
} elsif (\@_ > 2) {\$c=\$_[2]} else {\$c=PDL->nullcreate(\$x)} } elsif (\@_ > 2) {\$c=\$_[2]} else {\$c=PDL->nullcreate(\$x)}
PDL::_${name}_int(\$x,\$y,\$c); PDL::_${name}_int(\$x,\$y,\$c);
skipping to change at line 1020 skipping to change at line 964
wrapper around L</hclip> and wrapper around L</hclip> and
L</lclip>. L</lclip>.
=cut =cut
EOD EOD
pp_def( pp_def(
'clip', 'clip',
HandleBad => 1, HandleBad => 1,
Pars => 'a(); l(); h(); [o] c()', Pars => 'a(); l(); h(); [o] c()',
Code => Code => <<'EOBC',
'$c() = PDLMIN($h(), PDLMAX($l(), $a()));', PDL_IF_BAD(
BadCode => <<'EOBC',
if( $ISBAD(a()) || $ISBAD(l()) || $ISBAD(h()) ) { if( $ISBAD(a()) || $ISBAD(l()) || $ISBAD(h()) ) {
$SETBAD(c()); $SETBAD(c());
} else { } else,) {
$c() = PDLMIN($h(), PDLMAX($l(), $a())); $c() = PDLMIN($h(), PDLMAX($l(), $a()));
} }
EOBC EOBC
PMCode => <<'EOPM', PMCode => <<'EOPM',
*clip = \&PDL::clip; *clip = \&PDL::clip;
sub PDL::clip { sub PDL::clip {
my($x, $l, $h) = @_; my($x, $l, $h) = @_;
my $d; my $d;
unless(defined($l) || defined($h)) { unless(defined($l) || defined($h)) {
# Deal with pathological case # Deal with pathological case
skipping to change at line 1079 skipping to change at line 1022
pp_def( pp_def(
'wtstat', 'wtstat',
HandleBad => 1, HandleBad => 1,
Pars => 'a(n); wt(n); avg(); [o]b();', Pars => 'a(n); wt(n); avg(); [o]b();',
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
OtherPars => 'int deg', OtherPars => 'int deg',
Code => Code =>
'double wtsum = 0; 'double wtsum = 0;
double statsum = 0; double statsum = 0;
loop(n) %{
register double tmp;
register int i;
wtsum += $wt();
tmp=1;
for(i=0; i<$COMP(deg); i++)
tmp *= $a();
statsum += $wt() * (tmp - $avg());
%}
$b() = statsum / wtsum;',
BadCode =>
'double wtsum = 0;
double statsum = 0;
int flag = 0; int flag = 0;
loop(n) %{ loop(n) %{
if ( $ISGOOD(wt()) && $ISGOOD(a()) && $ISGOOD(avg()) ) { PDL_IF_BAD(if ($ISGOOD(wt()) && $ISGOOD(a()) && $ISGOOD(avg())),) {
register double tmp; long double tmp;
register int i; PDL_Indx i;
wtsum += $wt(); wtsum += $wt();
tmp=1; tmp=1;
for(i=0; i<$COMP(deg); i++) for(i=0; i<$COMP(deg); i++)
tmp *= $a(); tmp *= $a();
statsum += $wt() * (tmp - $avg()); statsum += $wt() * (tmp - $avg());
flag = 1; flag = 1;
} }
%} %}
if ( flag ) { $b() = statsum / wtsum; } PDL_IF_BAD(if (!flag) { $SETBAD(b()); $PDLSTATESETBAD(b); }
else { $SETBAD(b()); $PDLSTATESETBAD(b); }', else,) { $b() = statsum / wtsum; }',
CopyBadStatusCode => '', CopyBadStatusCode => '',
Doc => ' Doc => '
=for ref =for ref
Weighted statistical moment of given degree Weighted statistical moment of given degree
This calculates a weighted statistic over the vector C<a>. This calculates a weighted statistic over the vector C<a>.
The formula is The formula is
b() = (sum_i wt_i * (a_i ** degree - avg)) / (sum_i wt_i) b() = (sum_i wt_i * (a_i ** degree - avg)) / (sum_i wt_i)
skipping to change at line 1130 skipping to change at line 1060
Bad values are ignored in any calculation; C<$b> will only Bad values are ignored in any calculation; C<$b> will only
have its bad flag set if the output contains any bad data. have its bad flag set if the output contains any bad data.
', ',
); );
pp_def('statsover', pp_def('statsover',
HandleBad => 1, HandleBad => 1,
Pars => 'a(n); w(n); float+ [o]avg(); float+ [o]prms(); int+ [o]median(); int+ [o]min(); int+ [o]max(); float+ [o]adev(); float+ [o]rms()', Pars => 'a(n); w(n); float+ [o]avg(); float+ [o]prms(); int+ [o]median(); int+ [o]min(); int+ [o]max(); float+ [o]adev(); float+ [o]rms()',
Code => Code =>
'$GENERIC(avg) tmp = 0; '$GENERIC(avg) tmp = 0;
$GENERIC(avg) tmp1 = 0;
$GENERIC(avg) diff = 0;
$GENERIC(min) curmin = 0, curmax = 0;
$GENERIC(avg) norm = 0;
loop(n) %{ /* Accumulate sum and summed weight. */
tmp += $a()*$w();
norm += ($GENERIC(avg)) $w();
if (!n) { curmin = $a(); curmax = $a();}
if ($a() < curmin) {
curmin = $a();
} else if ($a() > curmax) {
curmax = $a();
}
%}
$avg() = tmp / norm; /* Find mean */
$min() = curmin;
$max() = curmax;
/* Calculate the RMS and standard deviation. */
tmp = 0;
loop(n) %{
diff = ($a() - $avg());
tmp += diff * diff * $w();
tmp1 += fabs(diff) * $w();
%}
$rms() = sqrt ( tmp/norm );
$prms() = (norm>1) ? sqrt( tmp/(norm-1) ) : 0;
$adev() = tmp1/norm ;
',
BadCode =>
'$GENERIC(avg) tmp = 0;
$GENERIC(avg) tmp1 = 0; $GENERIC(avg) tmp1 = 0;
$GENERIC(avg) diff = 0; $GENERIC(avg) diff = 0;
$GENERIC(min) curmin = 0, curmax = 0; $GENERIC(min) curmin = 0, curmax = 0;
$GENERIC(w) norm = 0; $GENERIC(w) norm = 0;
int flag = 0; int flag = 0;
loop(n) %{ loop(n) %{ /* Accumulate sum and summed weight. */
/* perhaps should check w() for bad values too ? */ /* perhaps should check w() for bad values too ? */
if ( $ISGOOD(a()) ) { PDL_IF_BAD(if ( $ISGOOD(a()) ),) {
tmp += $a()*$w(); tmp += $a()*$w();
norm += $w(); norm += ($GENERIC(avg)) $w();
if (!flag) { curmin = $a(); curmax = $a(); flag=1; } if (!flag) { curmin = $a(); curmax = $a(); flag=1; }
if ($a() < curmin) { if ($a() < curmin) {
curmin = $a(); curmin = $a();
} else if ($a() > curmax) { } else if ($a() > curmax) {
curmax = $a(); curmax = $a();
} }
} }
%} %}
/* have at least one valid point if flag == 1 */ /* have at least one valid point if flag true */
if ( flag ) { PDL_IF_BAD(if ( !flag ) {
$SETBAD(avg()); $PDLSTATESETBAD(avg);
$SETBAD(rms()); $PDLSTATESETBAD(rms);
$SETBAD(adev()); $PDLSTATESETBAD(adev);
$SETBAD(min()); $PDLSTATESETBAD(min);
$SETBAD(max()); $PDLSTATESETBAD(max);
$SETBAD(prms()); $PDLSTATESETBAD(prms);
} else,) {
$avg() = tmp / norm; /* Find mean */ $avg() = tmp / norm; /* Find mean */
$min() = curmin; $min() = curmin;
$max() = curmax; $max() = curmax;
/* Calculate the RMS and standard deviation. */
/* Calculate the RMS and standard deviation. */
tmp = 0; tmp = 0;
loop(n) %{ loop(n) %{
if ($ISGOOD(a())) { if ($ISGOOD(a())) {
diff = $a()-$avg(); diff = $a()-$avg();
tmp += diff * diff * $w(); tmp += diff * diff * $w();
tmp1 += fabs(diff) * $w(); tmp1 += fabs(diff) * $w();
} }
%} %}
$rms() = sqrt( tmp/norm ); $rms() = sqrt( tmp/norm );
if(norm>1) if(norm>1)
$prms() = sqrt( tmp/(norm-1) ); $prms() = sqrt( tmp/(norm-1) );
else else
$SETBAD(prms()); PDL_IF_BAD($SETBAD(prms()),$prms() = 0);
$adev() = tmp1 / norm ; $adev() = tmp1 / norm ;
} else {
$SETBAD(avg()); $PDLSTATESETBAD(avg);
$SETBAD(rms()); $PDLSTATESETBAD(rms);
$SETBAD(adev()); $PDLSTATESETBAD(adev);
$SETBAD(min()); $PDLSTATESETBAD(min);
$SETBAD(max()); $PDLSTATESETBAD(max);
$SETBAD(prms()); $PDLSTATESETBAD(prms);
}', }',
CopyBadStatusCode => '', CopyBadStatusCode => '',
PMCode => ' PMCode => '
sub PDL::statsover { sub PDL::statsover {
barf(\'Usage: ($mean,[$prms, $median, $min, $max, $adev, $rms]) = statsover($ data,[$weights])\') if @_>2; barf(\'Usage: ($mean,[$prms, $median, $min, $max, $adev, $rms]) = statsover($ data,[$weights])\') if @_>2;
my ($data, $weights) = @_; my ($data, $weights) = @_;
$weights = $data->ones() if !defined($weights); $weights //= $data->ones();
my $median = $data->medover;
my $median = $data->medover();
my $mean = PDL->nullcreate($data); my $mean = PDL->nullcreate($data);
my $rms = PDL->nullcreate($data); my $rms = PDL->nullcreate($data);
my $min = PDL->nullcreate($data); my $min = PDL->nullcreate($data);
my $max = PDL->nullcreate($data); my $max = PDL->nullcreate($data);
my $adev = PDL->nullcreate($data); my $adev = PDL->nullcreate($data);
my $prms = PDL->nullcreate($data); my $prms = PDL->nullcreate($data);
PDL::_statsover_int($data, $weights, $mean, $prms, $median, $min, $max, $adev , $rms); PDL::_statsover_int($data, $weights, $mean, $prms, $median, $min, $max, $adev , $rms);
wantarray ? ($mean, $prms, $median, $min, $max, $adev, $rms) : $mean;
return $mean unless wantarray;
return ($mean, $prms, $median, $min, $max, $adev, $rms);
} }
', ',
Doc => ' Doc => '
=for ref =for ref
Calculate useful statistics over a dimension of an ndarray Calculate useful statistics over a dimension of an ndarray
=for usage =for usage
($mean,$prms,$median,$min,$max,$adev,$rms) = statsover($ndarray, $weights); ($mean,$prms,$median,$min,$max,$adev,$rms) = statsover($ndarray, $weights);
This utility function calculates various useful This utility function calculates various useful
quantities of an ndarray. These are: quantities of an ndarray. These are:
skipping to change at line 1288 skipping to change at line 1180
=back =back
This operator is a projection operator so the calculation This operator is a projection operator so the calculation
will take place over the final dimension. Thus if the input will take place over the final dimension. Thus if the input
is N-dimensional each returned value will be N-1 dimensional, is N-dimensional each returned value will be N-1 dimensional,
to calculate the statistics for the entire ndarray either to calculate the statistics for the entire ndarray either
use C<clump(-1)> directly on the ndarray or call C<stats>. use C<clump(-1)> directly on the ndarray or call C<stats>.
=cut =cut
', ',
BadDoc =>' BadDoc =>'
=for bad =for bad
Bad values are simply ignored in the calculation, effectively reducing Bad values are simply ignored in the calculation, effectively reducing
the sample size. If all data are bad then the output data are marked bad. the sample size. If all data are bad then the output data are marked bad.
=cut =cut
', ',
); );
pp_add_exported('','stats'); pp_add_exported('','stats');
pp_addpm(<<'EOD'); pp_addpm(<<'EOD');
=head2 stats =head2 stats
=for ref =for ref
Calculates useful statistics on an ndarray Calculates useful statistics on an ndarray
skipping to change at line 1426 skipping to change at line 1315
Doc => $histogram_doc, Doc => $histogram_doc,
}, },
{Name => 'whistogram', {Name => 'whistogram',
WeightPar => 'float+ wt(n);', WeightPar => 'float+ wt(n);',
HistType => 'float+', HistType => 'float+',
HistOp => '+= $wt()', HistOp => '+= $wt()',
Doc => $whistogram_doc, Doc => $whistogram_doc,
} }
) )
{ {
my $p1 = 'register int j; my $code = 'register int j;
register int maxj = $SIZE(m)-1; register int maxj = $SIZE(m)-1;
register double min = $COMP(min); register double min = $COMP(min);
register double step = $COMP(step); register double step = $COMP(step);
broadcastloop %{ broadcastloop %{
loop(m) %{ $hist() = 0; %} loop(m) %{ $hist() = 0; %}
%} %}
broadcastloop %{ broadcastloop %{
loop(n) %{'; loop(n) %{
my $p2 = ' PDL_IF_BAD(if ( $ISGOOD(in()) ),) {
j = (int) (($in()-min)/step); j = (int) (($in()-min)/step);
if (j<0) j=0; if (j<0) j=0;
if (j > maxj) j = maxj; if (j > maxj) j = maxj;
($hist(m => j))'.$_->{HistOp}.';'; ($hist(m => j))'.$_->{HistOp}.';
my $p3 = ' }
%} %}
%}'; %}';
pp_def($_->{Name}, pp_def($_->{Name},
Pars => 'in(n); '.$_->{WeightPar}.$_->{HistType}. '[o] hist(m)', Pars => 'in(n); '.$_->{WeightPar}.$_->{HistType}. '[o] hist(m)',
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
# set outdim by Par! # set outdim by Par!
OtherPars => 'double step; double min; int msize => m', OtherPars => 'double step; double min; int msize => m',
HandleBad => 1, HandleBad => 1,
Code => $p1.$p2.$p3, Code => $code,
BadCode => $p1.'if ( $ISGOOD(in()) ) {'.$p2.'}'.$p3, Doc=>$_->{Doc});
Doc=>$_->{Doc});
} }
my $histogram2d_doc = <<'EOD'; my $histogram2d_doc = <<'EOD';
=for ref =for ref
Calculates a 2d histogram. Calculates a 2d histogram.
=for usage =for usage
$h = histogram2d($datax, $datay, $stepx, $minx, $h = histogram2d($datax, $datay, $stepx, $minx,
skipping to change at line 1537 skipping to change at line 1425
Doc => $histogram2d_doc, Doc => $histogram2d_doc,
}, },
{Name => 'whistogram2d', {Name => 'whistogram2d',
WeightPar => 'float+ wt(n);', WeightPar => 'float+ wt(n);',
HistType => 'float+', HistType => 'float+',
HistOp => '+= $wt()', HistOp => '+= $wt()',
Doc => $whistogram2d_doc, Doc => $whistogram2d_doc,
} }
) )
{ {
my $p1 = 'register int ja,jb; my $code = 'register int ja,jb;
register int maxja = $SIZE(ma)-1; register int maxja = $SIZE(ma)-1;
register int maxjb = $SIZE(mb)-1; register int maxjb = $SIZE(mb)-1;
register double mina = $COMP(mina); register double mina = $COMP(mina);
register double minb = $COMP(minb); register double minb = $COMP(minb);
register double stepa = $COMP(stepa); register double stepa = $COMP(stepa);
register double stepb = $COMP(stepb); register double stepb = $COMP(stepb);
broadcastloop %{ broadcastloop %{
loop(ma,mb) %{ $hist() = 0; %} loop(ma,mb) %{ $hist() = 0; %}
%} %}
broadcastloop %{ broadcastloop %{
loop(n) %{'; loop(n) %{
my $p2 = ' PDL_IF_BAD(if ( $ISGOOD(ina()) && $ISGOOD(inb()) ),) {
ja = (int) (($ina()-mina)/stepa); ja = (int) (($ina()-mina)/stepa);
jb = (int) (($inb()-minb)/stepb); jb = (int) (($inb()-minb)/stepb);
if (ja<0) ja=0; if (ja<0) ja=0;
if (ja > maxja) ja = maxja; if (ja > maxja) ja = maxja;
if (jb<0) jb=0; if (jb<0) jb=0;
if (jb > maxjb) jb = maxjb; if (jb > maxjb) jb = maxjb;
($hist(ma => ja,mb => jb))'.$_->{HistOp}.';'; ($hist(ma => ja,mb => jb))'.$_->{HistOp}.';
my $p3 = ' }
%} %}
%}'; %}';
pp_def($_->{Name}, pp_def($_->{Name},
Pars => 'ina(n); inb(n); '.$_->{WeightPar}.$_->{HistType}. '[o] hist(ma, mb)', Pars => 'ina(n); inb(n); '.$_->{WeightPar}.$_->{HistType}. '[o] hist(ma, mb)',
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
# set outdim by Par! # set outdim by Par!
OtherPars => 'double stepa; double mina; int masize => ma; OtherPars => 'double stepa; double mina; int masize => ma;
double stepb; double minb; int mbsize => mb;', double stepb; double minb; int mbsize => mb;',
HandleBad => 1, HandleBad => 1,
Code => $p1.$p2.$p3, Code => $code,
BadCode => $p1.'if ( $ISGOOD(ina()) && $ISGOOD(inb()) ) {'.$p2.'}'.$p3, Doc=> $_->{Doc});
Doc=> $_->{Doc});
} }
########################################################### ###########################################################
# a number of constructors: fibonacci, append, axisvalues & # a number of constructors: fibonacci, append, axisvalues &
# random numbers # random numbers
########################################################### ###########################################################
pp_def('fibonacci', pp_def('fibonacci',
Pars => 'i(n); indx [o]x(n)', Pars => 'i(n); indx [o]x(n)',
Inplace => 1, Inplace => 1,
skipping to change at line 3041 skipping to change at line 2928
$outi //= $this->nullcreate; $outi //= $this->nullcreate;
$outni //= $this->nullcreate; $outni //= $this->nullcreate;
PDL::_which_both_int($this,$outi,$outni); PDL::_which_both_int($this,$outi,$outni);
return wantarray ? ($outi,$outni) : $outi; return wantarray ? ($outi,$outni) : $outi;
} }
*PDL::which_both = \&which_both; *PDL::which_both = \&which_both;
EOD EOD
} }
) )
{ {
my $p1 = $_->{Variables} .' my $code = $_->{Variables} .'
loop(n) %{ loop(n) %{
if ( $mask() '; if ( $mask() PDL_IF_BAD(&& $ISGOOD($mask()),) ) {
my $p2 = '
) {
$inds(m => dm) = n; $inds(m => dm) = n;
dm++; dm++;
}'.$_->{Elseclause} . "\n". }'.$_->{Elseclause}.'
' %}'; %}';
pp_def($_->{Name}, pp_def($_->{Name},
HandleBad => 1, HandleBad => 1,
Doc => $_->{Doc}, Doc => $_->{Doc},
Pars => $_->{Pars}, Pars => $_->{Pars},
GenericTypes => [ppdefs_all], GenericTypes => [ppdefs_all],
PMCode => $_->{PMCode}, PMCode => $_->{PMCode},
Code => $p1.$p2, Code => $code,
BadCode => $p1.' && $ISGOOD($mask())'.$p2,
# the next one is currently a dirty hack # the next one is currently a dirty hack
# this will probably break once dataflow is enabled again # this will probably break once dataflow is enabled again
# *unless* we have made sure that mask is physical by now!!! # *unless* we have made sure that mask is physical by now!!!
RedoDimsCode => ' RedoDimsCode => '
PDL_Indx sum = 0; PDL_Indx sum = 0;
/* not sure if this is necessary */ /* not sure if this is necessary */
pdl * dpdl = $PDL(mask); pdl * dpdl = $PDL(mask);
$GENERIC() *m_datap = (($GENERIC() *)(PDL_REPRP(dpdl))); $GENERIC() *m_datap = (($GENERIC() *)(PDL_REPRP(dpdl)));
PDL_Indx inc = PDL_REPRINC(dpdl,0); PDL_Indx inc = PDL_REPRINC(dpdl,0);
PDL_Indx offs = PDL_REPROFFS(dpdl); PDL_Indx offs = PDL_REPROFFS(dpdl);
 End of changes. 69 change blocks. 
214 lines changed or deleted 100 lines changed or added

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