(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 |
|---|
"Fossies" - the Fresh Open Source Software Archive 