The Perl Toolchain Summit needs more sponsors. If your company depends on Perl, please support this very important event.
# 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: