# 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
# note that version numbers with trailing zeroes (e.g, 0.010) created problems
# in some of the tests
pp_setversion( '0.015' );
# check for indx/PDL_Indx support
use PDL::Lite; # load PDL to check definedness of &PDL::indx
my $indx_type = defined( &PDL::indx ) ? 'indx' : 'int';
my $PDL_Indx_type = defined( &PDL::indx ) ? 'PDL_Indx' : 'int';
# _flatten_into
#
pp_def( '_flatten_into',
Pars => "in(m); $indx_type b(m); $indx_type [o] idx(m)",
OtherPars => "double step; double min; $PDL_Indx_type 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 '.$PDL_Indx_type.' j;
register '.$PDL_Indx_type.' 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 '.$PDL_Indx_type.' j;
register '.$PDL_Indx_type.' 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
# modelled after setvaltobad()
#
pp_def( '_setnulltobad',
Pars => "in(); $indx_type count(); [o] out()",
HandleBad => 1,
Inplace => [ 'in' ],
Doc => 'Set empty bins to the bad value.',
CopyBadStatusCode => '
if( in == out && $ISPDLSTATEGOOD(in) ) {
PDL->propogate_badflag( out, 1 );
}
$SETPDLSTATEBAD(out);
',
Code => '
if( ! $count() ) { $SETBAD(out()); }
else { $out() = $in(); }
',
);
# _icount_loop
#
pp_def( '_icount_loop',
Pars => "in(n); $indx_type idx(n); $indx_type [o] out(m)",
OtherPars => "$PDL_Indx_type msize => m",
HandleBad => 1,
Doc => 'Count the number of elements in each bin.
This function returns a piddle of type I<indx> if you have a 64-bit-capable
PDL, otherwise it returns a piddle of type I<long>.',
Code => '
loop(n) %{
register '.$PDL_Indx_type.' j = $idx();
++( $out(m => j) );
%}
',
BadCode => '
loop(n) %{
if( ($ISBAD(idx())) || ($ISBAD(in())) ) continue;
register '.$PDL_Indx_type.' j = $idx();
++( $out(m => j) );
%}
',
);
# _imax_loop
#
pp_def( '_imax_loop',
Pars => "in(n); $indx_type idx(n); [o] out(m)",
OtherPars => "$PDL_Indx_type msize => m",
HandleBad => 1,
Doc => 'Find the maximum in each bin.',
Code => '
loop(n) %{
register '.$PDL_Indx_type.' j = $idx();
if( ($ISBAD(out(m => j))) || ($in() > $out(m => j)) ) {
$out(m => j) = $in();
}
%}
',
BadCode => '
loop(n) %{
if( ($ISBAD(idx())) || ($ISBAD(in())) ) continue;
register '.$PDL_Indx_type.' j = $idx();
if( ($ISBAD(out(m => j))) || ($in() > $out(m => j)) ) {
$out(m => j) = $in();
}
%}
',
);
# _imin_loop
#
pp_def( '_imin_loop',
Pars => "in(n); $indx_type idx(n); [o] out(m)",
OtherPars => "$PDL_Indx_type msize => m",
HandleBad => 1,
Doc => 'Find the minimum in each bin.',
Code => '
loop(n) %{
register '.$PDL_Indx_type.' j = $idx();
if( ($ISBAD(out(m => j))) || ($in() < $out(m => j)) ) {
$out(m => j) = $in();
}
%}
',
BadCode => '
loop(n) %{
if( ($ISBAD(idx())) || ($ISBAD(in())) ) continue;
register '.$PDL_Indx_type.' j = $idx();
if( ($ISBAD(out(m => j))) || ($in() < $out(m => j)) ) {
$out(m => j) = $in();
}
%}
',
);
# _isum_loop
#
pp_def( '_isum_loop',
Pars => "in(n); $indx_type idx(n); int+ [o] out(m); $indx_type [t] count(m)",
OtherPars => "$PDL_Indx_type 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 '.$PDL_Indx_type.' j = $idx();
$out(m => j) += $in();
++( $count(m => j) );
%}
',
BadCode => '
loop(n) %{
if( ($ISBAD(idx())) || ($ISBAD(in())) ) continue;
register '.$PDL_Indx_type.' j = $idx();
$out(m => j) += $in();
++( $count(m => j) );
%}
',
);
# _iavg_loop
#
pp_def( '_iavg_loop',
Pars => "in(n); $indx_type idx(n); double [o] out(m); $indx_type [t] count(m)",
OtherPars => "$PDL_Indx_type 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 '.$PDL_Indx_type.' j = $idx();
$out(m => j) += ( $in() - $out(m => j) ) / ++( $count(m => j) );
%}
',
BadCode => '
loop(n) %{
if( ($ISBAD(idx())) || ($ISBAD(in())) ) continue;
register '.$PDL_Indx_type.' j = $idx();
$out(m => j) += ( $in() - $out(m => j) ) / ++( $count(m => j) );
%}
',
);
# _istddev_loop
#
pp_def ( '_istddev_loop',
Pars => "in(n); $indx_type idx(n); double [o] out(m); $indx_type [t] count(m); double [t] avg(m)",
OtherPars => "$PDL_Indx_type 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 '.$PDL_Indx_type.' 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 '.$PDL_Indx_type.' 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 => "in(); $indx_type count(); [o] out()",
HandleBad => 1,
Inplace => [ 'in' ],
Doc => 'Finalization for _istddev_loop().',
CopyBadStatusCode => '
if( in == out && $ISPDLSTATEGOOD(in) ) {
PDL->propogate_badflag( out, 1 );
}
$SETPDLSTATEBAD(out);
',
Code => '
if( ! $count() ) { $SETBAD(out()); }
else { $out() = sqrt( $in() / $count() ); }
',
);
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: