# 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
our $VERSION = '0.020';
pp_setversion( $VERSION );
# _flatten_into
#
pp_def( '_flatten_into',
Pars => "in(m); indx b(m); indx [o] idx(m)",
OtherPars => "double step; double min; PDL_Indx 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 j;
register PDL_Indx maxj = $COMP( n ) - 1;
loop(m) %{
PDL_IF_BAD(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();
%}
',
);
# _flatten_into_grid
#
pp_def( '_flatten_into_grid',
Pars => "in(); indx b(); double edge(n); indx [o] idx()",
HandleBad => 1,
Doc => 'Bin a series of values on a specified grid and flatten into an existing list of indices.',
Code => '
PDL_Indx n = $SIZE(n);
PDL_Indx n1 = n-1;
PDL_Indx mid;
/* determine sort order of edges */
int up = ($edge(n => n1) >= $edge(n => 0));
broadcastloop %{
PDL_Indx low = 0;
PDL_Indx high = n1;
$GENERIC() value = $in();
PDL_IF_BAD(if( ($ISBAD(in())) || ($ISBAD(b())) ) {
$SETBAD( idx() );
continue;
},)
while (low <= high) {
/* ensure we do not overflow (low+high)>>1 for large values of low + high */
mid = low + ((high - low) >> 1);
if (($edge(n => mid) <= value) == up) low = mid + 1;
else high = mid - 1;
}
$idx() = (up
? high < 0 ? 0 : high
: low > n1 ? n1 : low
) + n * $b();
%}
',
);
# _setnulltobad
# modelled after setvaltobad()
#
pp_def( '_setnulltobad',
Pars => "in(); indx count(); [o] out()",
HandleBad => 1,
Inplace => [ 'in' ],
Doc => 'Set empty bins to the bad value.',
Code => '
char anybad = 0;
broadcastloop %{
if( ! $count() ) { anybad = 1; $SETBAD(out()); }
else { $out() = $in(); }
%}
if (anybad) $PDLSTATESETBAD(out);
',
);
# _icount_loop
#
pp_def( '_icount_loop',
Pars => "in(n); indx idx(n); indx [o] out(m)",
OtherPars => "PDL_Indx msize => m",
HandleBad => 1,
Doc => 'Count the number of elements in each bin.
This function returns an ndarray of type I<indx> if you have a 64-bit-capable
PDL, otherwise it returns an ndarray of type I<long>.',
Code => '
loop(n) %{
PDL_IF_BAD(if( ($ISBAD(idx())) || ($ISBAD(in())) ) continue;,)
register PDL_Indx j = $idx();
++( $out(m => j) );
%}
',
);
# _imax_loop
#
pp_def( '_imax_loop',
Pars => "in(n); indx idx(n); [o] out(m)",
OtherPars => "PDL_Indx msize => m",
HandleBad => 1,
Doc => 'Find the maximum in each bin.',
Code => '
loop(n) %{
PDL_IF_BAD(if ($ISBAD(idx()) || $ISBAD(in())) continue;,)
register PDL_Indx 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 idx(n); [o] out(m)",
OtherPars => "PDL_Indx msize => m",
HandleBad => 1,
Doc => 'Find the minimum in each bin.',
Code => '
loop(n) %{
if (PDL_IF_BAD($ISBAD(idx()) ||,) $ISBAD(in())) continue;
register PDL_Indx 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 idx(n); int+ [o] out(m); indx [o] count(m)",
OtherPars => "PDL_Indx msize => m",
HandleBad => 1,
Doc => 'Sum the elements in each bin.
This function returns an ndarray of type I<int> or higher, to reduce the risk of
overflow when collecting sums.',
Code => '
loop(n) %{
PDL_IF_BAD(if ($ISBAD(idx()) || $ISBAD(in())) continue;,)
register PDL_Indx j = $idx();
$out(m => j) += $in();
++( $count(m => j) );
%}
',
);
# _iavg_loop
#
pp_def( '_iavg_loop',
Pars => "in(n); indx idx(n); double [o] out(m); indx [o] count(m)",
OtherPars => "PDL_Indx 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) %{
PDL_IF_BAD(if ($ISBAD(idx()) || $ISBAD(in())) continue;,)
register PDL_Indx j = $idx();
$out(m => j) += ( $in() - $out(m => j) ) / ++( $count(m => j) );
%}
',
);
# _istddev_loop
#
pp_def ( '_istddev_loop',
Pars => "in(n); indx idx(n); double [o] out(m); indx [o] count(m); double [o] avg(m)",
OtherPars => "PDL_Indx 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) %{
PDL_IF_BAD(if ($ISBAD(idx()) || $ISBAD(in())) continue;,)
register PDL_Indx 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 count(); [o] out()",
HandleBad => 1,
Inplace => [ 'in' ],
Doc => 'Finalization for _istddev_loop().',
Code => '
char anybad = 0;
broadcastloop %{
if( ! $count() ) { anybad = 1; $SETBAD(out()); }
else { $out() = sqrt( $in() / $count() ); }
%}
if (anybad) $PDLSTATESETBAD(out);
',
);
pp_addpm( <<'EOD' );
=head1 AUTHOR
Edward Baudrez <ebaudrez@cpan.org>
=head1 COPYRIGHT AND LICENSE
This software is copyright (c) 2014 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: