# use pp_addbegin() to force =head1 NAME to appear before =head1 FUNCTIONS
pp_addbegin( <<'EOD' );
# ABSTRACT: XS actions for PDL::NDBin
=head1 NAME
PDL::NDBin::Actions_PP - XS actions for PDL::NDBin
=head1 DESCRIPTION
This module contains the implementation of the core loop of the action classes
in the PDL::NDBin::Action namespace. The subroutines are not intended to be
called by the user.
=cut
# silence Test::Pod::Coverage warning
=for Pod::Coverage set_boundscheck set_debugging
=cut
EOD
pp_setversion( '0.008' );
# _flatten_into
#
pp_def( '_flatten_into',
Pars => 'in(m); int b(m); int [o] idx(m)',
OtherPars => 'double step; double min; int n',
HandleBad => 1,
Doc => 'Bin a series of values and flatten into an existing list of indices.',
Code => '
register double min = $COMP( min );
register double step = $COMP( step );
register int j;
register int maxj = $COMP( n ) - 1;
loop(m) %{
j = (( $in() - min )/step);
if( j < 0 ) j = 0;
if( j > maxj ) j = maxj;
$idx() = j + $COMP( n ) * $b();
%}
',
BadCode => '
register double min = $COMP( min );
register double step = $COMP( step );
register int j;
register int maxj = $COMP( n ) - 1;
loop(m) %{
if( ($ISBAD(in())) || ($ISBAD(b())) ) {
$SETBAD( idx() );
continue;
}
j = (( $in() - min )/step);
if( j < 0 ) j = 0;
if( j > maxj ) j = maxj;
$idx() = j + $COMP( n ) * $b();
%}
',
);
# _setnulltobad
#
pp_def( '_setnulltobad',
Pars => 'count(n); [o] out(n)',
HandleBad => 1,
Doc => 'Set empty bins to the bad value.',
Code => '
int flag = 0;
loop(n) %{
if( ! $count() ) { $SETBAD(out()); flag = 1; }
%}
if( flag ) { $PDLSTATESETBAD(out); }
',
BadCode => '
int flag = 0;
loop(n) %{
if( ! $count() ) { $SETBAD(out()); flag = 1; }
%}
if( flag ) { $PDLSTATESETBAD(out); }
',
);
# _icount_loop
#
pp_def( '_icount_loop',
Pars => 'in(n); int idx(n); int [o] out(m)',
OtherPars => 'int msize => m',
HandleBad => 1,
Doc => 'Count the number of elements in each bin.
This function always returns a piddle of type I<int>.',
Code => '
loop(n) %{
register int j = $idx();
++( $out(m => j) );
%}
',
BadCode => '
loop(n) %{
if( ($ISBAD(idx())) || ($ISBAD(in())) ) continue;
register int j = $idx();
++( $out(m => j) );
%}
',
);
# _isum_loop
#
pp_def( '_isum_loop',
Pars => 'in(n); int idx(n); int+ [o] out(m); int [t] count(m)',
OtherPars => 'int msize => m',
HandleBad => 1,
Doc => 'Sum the elements in each bin.
This function returns a piddle of type I<int> or higher, to reduce the risk of
overflow when collecting sums.',
Code => '
loop(n) %{
register int j = $idx();
$out(m => j) += $in();
++( $count(m => j) );
%}
',
BadCode => '
loop(n) %{
if( ($ISBAD(idx())) || ($ISBAD(in())) ) continue;
register int j = $idx();
$out(m => j) += $in();
++( $count(m => j) );
%}
',
);
# _iavg_loop
#
pp_def( '_iavg_loop',
Pars => 'in(n); int idx(n); double [o] out(m); int [t] count(m)',
OtherPars => 'int msize => m',
HandleBad => 1,
Doc => q[Compute the average of the elements in each bin.
Credit for the algorithm goes to
L<I<ashawley> on commandlinefu.com|http://www.commandlinefu.com/commands/view/3437/compute-running-average-for-a-column-of-numbers>:
awk '{ avg += ($1 - avg) / NR } END { print avg }'
This is a wonderful solution solving many of the problems with more naive
implementations:
=over 4
=item 1.
It's numerically well-behaved: out() is always of the order of magnitude of the
values themselves, unlike the sum of the values, which grows very large as the
number of elements grows large.
=item 2.
The subtraction in() - out() guarantees that the computation will be done in
the correct type (i.e., I<double> instead of the type of the input).
=item 3.
Requires only one temporary.
=item 4.
Requires only one pass over the data.
=back
I used to give the output array type I<float+>, but that creates more problems
than it solves. So now, averages are always computed in type I<double>. Besides
being the default type in PDL and the 'natural' floating-point type in C, it
also makes the implementation easier.],
Code => '
loop(n) %{
register int j = $idx();
$out(m => j) += ( $in() - $out(m => j) ) / ++( $count(m => j) );
%}
',
BadCode => '
loop(n) %{
if( ($ISBAD(idx())) || ($ISBAD(in())) ) continue;
register int j = $idx();
$out(m => j) += ( $in() - $out(m => j) ) / ++( $count(m => j) );
%}
',
);
# _istddev_loop
#
pp_def ( '_istddev_loop',
Pars => 'in(n); int idx(n); double [o] out(m); int [t] count(m); double [t] avg(m)',
OtherPars => 'int msize => m',
HandleBad => 1,
Doc => q[Compute the standard deviation of the elements in each bin.
Note that we compute the sample standard deviation, I<not> an estimate of the
population standard deviation (which differs by a factor).
Credit for the algorithm goes to
L<I<ashawley> on commandlinefu.com|http://www.commandlinefu.com/commands/view/3442/display-the-standard-deviation-of-a-column-of-numbers-with-awk>;
awk '{ delta = $1 - avg; avg += delta / NR; mean2 += delta * ($1 - avg) } END { print sqrt(mean2 / NR) }'
This is a wonderful solution solving many of the problems with more naive
implementations:
=over 4
=item 1.
It's numerically well-behaved.
=item 2.
The subtractions guarantee that the computations will be done in the correct
type (i.e., I<double> instead of the type of the input).
=item 3.
Requires only two temporaries (!).
=item 4.
Requires only one pass over the data.
=back
I used to give the output array type I<float+>, but that creates more problems
than it solves. So now, standard deviations are always computed in type
I<double>. Besides being the default type in PDL and the 'natural'
floating-point type in C, it also makes the implementation easier.],
Code => '
loop(n) %{
register int j = $idx();
double delta = $in() - $avg(m => j);
$avg(m => j) += delta / ++( $count(m => j) );
$out(m => j) += delta * ( $in() - $avg(m => j) );
%}
',
BadCode => '
loop(n) %{
if( ($ISBAD(idx())) || ($ISBAD(in())) ) continue;
register int j = $idx();
double delta = $in() - $avg(m => j);
$avg(m => j) += delta / ++( $count(m => j) );
$out(m => j) += delta * ( $in() - $avg(m => j) );
%}
',
);
# _istddev_post
#
pp_def( '_istddev_post',
Pars => 'count(n); [o] out(n)',
HandleBad => 1,
Doc => 'Finalization for _istddev_loop().',
Code => '
int flag = 0;
loop(n) %{
if( ! $count() ) { $SETBAD(out()); flag = 1; }
else { $out() = sqrt( $out() / $count() ); }
%}
if( flag ) { $PDLSTATESETBAD(out); }
',
BadCode => '
int flag = 0;
loop(n) %{
if( ! $count() ) { $SETBAD(out()); flag = 1; }
else { $out() = sqrt( $out() / $count() ); }
%}
if( flag ) { $PDLSTATESETBAD(out); }
',
);
# TODO
# complete with the other functions that xyz2grd provides
# indirect maximum
# indirect minimum
pp_addpm( <<'EOD' );
=head1 AUTHOR
Edward Baudrez <ebaudrez@cpan.org>
=head1 COPYRIGHT AND LICENSE
This software is copyright (c) 2013 by Edward Baudrez.
This is free software; you can redistribute it and/or modify it under
the same terms as the Perl 5 programming language system itself.
=cut
EOD
pp_done();
# vim:set filetype=perl: