File Coverage

blib/lib/PDL/NDBin.pm

Criterion	Covered	Total	%
statement	299	316	94.6
branch	119	154	77.2
condition	37	43	86.0
subroutine	43	43	100.0
pod	13	13	100.0
total	511	569	89.8

package PDL::NDBin;
# ABSTRACT: Multidimensional binning & histogramming
$PDL::NDBin::VERSION = '0.019';
use strict;
use warnings;
use Exporter;
use List::Util qw( reduce );
use List::MoreUtils qw( pairwise );
use Math::Round qw( nlowmult nhimult );
use PDL::Lite;		# do not import any functions into this namespace
use PDL::NDBin::Iterator;
use PDL::NDBin::Actions_PP;
use PDL::NDBin::Utils_PP;
use Log::Any qw( $log );
use Data::Dumper;
use UUID::Tiny qw( :std );
use POSIX qw( ceil );
use Params::Validate qw( validate validate_pos validate_with ARRAYREF CODEREF HASHREF SCALAR );
use Carp;
use Class::Load qw( load_class );


our @ISA = qw( Exporter );
our @EXPORT = qw( );
our @EXPORT_OK = qw( ndbinning ndbin );
our %EXPORT_TAGS = ( all => [ qw( ndbinning ndbin ) ] );

# the list of valid keys
my %valid_key = map { $_ => 1 } qw( axes vars );


my @axis_params = qw( max min n round step grid );
my ( %axis_params, %axis_flags );
@axis_params{@axis_params} = (1) x @axis_params;
@axis_flags{@axis_params} = map { 1<<$_ } 0..@axis_params-1;

my %axis_allowed =
    map { reduce( sub { $a | $b }, 0, @axis_flags{@$_} ) => 1 }
        [ ],
	[ qw( n ) ],
	[ qw( min ) ],
	[ qw( max ) ],
        [ qw( step ) ],
	[ qw( min step ) ],
	[ qw( max step ) ],
	[ qw( min n ) ],
	[ qw( max n ) ],
	[ qw( round n ) ],
	[ qw( round step ) ],
	[ qw( n step ) ],
	[ qw( n step round ) ],
	[ qw( min n step ) ],
	[ qw( max n step ) ],
	[ qw( min max n ) ],
	[ qw( min max step ) ],
        [ qw( grid ) ];

sub add_axis
{
	my $self = shift;
	my %params = validate( @_, {
			max   => 0,
			min   => 0,
			n     => 0,
			name  => 1,
			pdl   => 0,
			round => 0,
			step  => 0,
			grid  => 0,
		} );
	$log->tracef( 'adding axis with specs %s', \%params );

	my $pmask = reduce { $a | ($b||0) } 0, @axis_flags{ keys %params };
	croak( "inconsistent or incomplete parameters: ", keys %params )
	    unless $axis_allowed{ $pmask };

	push @{ $self->{axes} }, \%params;
}


sub add_var
{
	my $self = shift;
	my %params = validate( @_, {
			action => { type => CODEREF | HASHREF | SCALAR },
			name   => 1,
			pdl    => 0,
		} );
	$log->tracef( 'adding variable with specs %s', \%params );
	push @{ $self->{vars} }, \%params;
}


sub new
{
	my $class = shift;
	my %params = validate( @_, {
			axes => { optional => 1, type => ARRAYREF },
			vars => { optional => 1, type => ARRAYREF },
		} );
	$log->debug( 'new: arguments = ' . Dumper \%params ) if $log->is_debug;
	my $self = bless { axes => [], vars => [] }, $class;
	# axes
	$params{axes} ||= [];		# be sure we can dereference
	my @axes = @{ $params{axes} };
	for my $axis ( @axes ) {
		my @pat = ( 1 );					# one mandatory argument
		if( @$axis > 1 ) { push @pat, (0) x (@$axis - 1) }	# followed by n-1 optional arguments
		my( $name ) = validate_pos( @$axis, @pat );
		shift @$axis; # remove name
		$self->add_axis( name => $name, @$axis );
	}
	# vars
	$params{vars} ||= [];		# be sure we can dereference
	my @vars = @{ $params{vars} };
	for my $var ( @vars ) {
		my( $name, $action ) = validate_pos( @$var, 1, 1 );
		$self->add_var( name => $name, action => $action );
	}
	return $self;
}


sub axes { wantarray ? @{ $_[0]->{axes} } : $_[0]->{axes} }
sub vars { wantarray ? @{ $_[0]->{vars} } : $_[0]->{vars} }

sub _make_instance_hashref
{
	my %params = validate_with(
		params => \@_,
		spec   => {
			N       => 1,
			class   => 1,
			coderef => 0,
		},
		allow_extra => 1,
	);
	my $short_class = delete $params{class};
	my $full_class = substr( $short_class, 0, 1 ) eq '+'
		? substr( $short_class, 1 )
		: "PDL::NDBin::Action::$short_class";
	load_class( $full_class );
	return $full_class->new( %params );
}

sub _make_instance
{
	my %params = validate( @_, {
			action => 1,
			N      => 1,
		} );
	if( ref $params{action} eq 'CODE' ) {
		return _make_instance_hashref(
			class   => '+PDL::NDBin::Action::CodeRef',
			N       => $params{N},
			coderef => $params{action},
		);
	}
	elsif( ref $params{action} eq 'HASH' ) {
		return _make_instance_hashref(
			%{ $params{action} },
			N => $params{N},
		);
	}
	else {
		return _make_instance_hashref(
			class => $params{action},
			N     => $params{N},
		);
	}
}


sub feed
{
	my $self = shift;
	my %pdls = @_;
	while( my( $name, $pdl ) = each %pdls ) {
		for my $v ( $self->axes, $self->vars ) {
			$v->{pdl} = $pdl if $v->{name} eq $name;
		}
	}
}

sub _check_all_pdls_present
{
	my $self = shift;
	my %warned_for;
	for my $v ( $self->axes, $self->vars ) {
		next if defined $v->{pdl};
		next if $v->{action} eq 'Count'; # those variables don't need data
		my $name = $v->{name};
		next if $warned_for{ $name };
		$log->error( "no data for $name" );
		$warned_for{ $name }++;
	}
}

sub _check_pdl_length
{
	my $self = shift;
	# checking whether the lengths of all axes and variables are equal can
	# only be done here (in a loop), and not in autoscale_axis()
	my $length;
	for my $v ( $self->axes, $self->vars ) {
		$length = $v->{pdl}->nelem unless defined $length;
		# variables don't always need a pdl, or may be happy with a
		# null pdl; let the action figure it out.
		# note that the test isempty() is not a good test for null
		# pdls, but until I have a better one, this will have to do
		next if $v->{action} && ( ! defined $v->{pdl} || $v->{pdl}->isempty );
		if( $v->{pdl}->nelem != $length ) {
			croak( join '', 'number of elements (',
				$v->{pdl}->nelem, ") of '$v->{name}'",
				" is different from previous ($length)" );
		}
	}
}


sub autoscale_axis
{
	my $axis = shift;
	# return early if step, min, and n have already been calculated
	if( defined $axis->{step} && defined $axis->{min} && defined $axis->{n} ) {
		$log->tracef( 'step, min, n already calculated for %s; not recalculating', $axis );
		return;
	}
	# first get & sanify the arguments
	croak( 'need coordinates' ) unless defined $axis->{pdl};
	# return if axis is empty
	if( $axis->{pdl}->isempty ) {
		$axis->{n} = 0;
		return;
	}

	# return early if a grid has been supplied
	if( defined $axis->{grid} ) {

	        $axis->{grid} = PDL::Core::topdl( $axis->{grid} );
	        croak( "grid supplied for %s must be one-dimensional with at least two elements", $axis )
		  if $axis->{grid}->nelem < 2 || $axis->{grid}->ndims > 1;
		_validate_grid( $axis->{grid} );
		# number of bins is one less than number of bin edges
		$axis->{n} = $axis->{grid}->nelem - 1;
		$log->tracef( 'grid supplied for %s; no need to autoscale', $axis );
		return;
	}


	$axis->{min} = $axis->{pdl}->min unless defined $axis->{min};


	$axis->{max} = $axis->{pdl}->max unless defined $axis->{max};
	if( defined $axis->{round} and $axis->{round} > 0 ) {
		$axis->{min} = nlowmult( $axis->{round}, $axis->{min} );
		$axis->{max} = nhimult(  $axis->{round}, $axis->{max} );
	}
	croak( 'max < min is invalid' ) if $axis->{max} < $axis->{min};
	if( $axis->{pdl}->type >= PDL::float ) {
		croak( 'cannot bin with min = max' ) if $axis->{min} == $axis->{max};
	}
	# calculate the range
	# for floating-point data, we need to augment the range by 1 unit - see
	# the discussion under IMPLEMENTATION NOTES for more details
	my $range = $axis->{max} - $axis->{min};
	if( $axis->{pdl}->type < PDL::float ) {
		$range += 1;
	}
	# if step size has been supplied by user, check it
	if( defined $axis->{step} ) {
		croak( 'step size must be > 0' ) unless $axis->{step} > 0;
		if( $axis->{pdl}->type < PDL::float && $axis->{step} < 1 ) {
			croak( "step size = $axis->{step} < 1 is not allowed when binning integral data" );
		}
	}
	# number of bins I
	if( defined $axis->{n} ) {
		croak( 'number of bins must be > 0' ) unless $axis->{n} > 0;
		croak( 'number of bins must be integral' ) if ceil( $axis->{n} ) - $axis->{n} > 0;
	}
	else {
		if( defined $axis->{step} ) {
			# data range and step size were verified above,
			# so the result of this calculation is
			# guaranteed to be > 0
			$axis->{n} = ceil( $range / $axis->{step} );
		}
		else {
			# if neither number of bins nor step size is defined,
			# use some reasonable default (which used to be the
			# behaviour of hist() in versions of PDL inferior to
			# 2.4.12) (see F)
			$axis->{n} = $axis->{pdl}->nelem > 100 ? 100 : $axis->{pdl}->nelem;
		}
	}
	# step size I
	# if we get here, the data range is certain to be larger than
	# zero, and I is sure to be defined and valid (either
	# because it was supplied explicitly and verified to be valid,
	# or because it was calculated automatically)
	if( ! defined $axis->{step} ) {
		# result of this calculation is guaranteed to be > 0
		$axis->{step} = $range / $axis->{n};
		if( $axis->{pdl}->type < PDL::float ) {
			croak( 'there are more bins than distinct values' ) if $axis->{step} < 1;
		}
	}
}


sub autoscale
{
	my $self = shift;
	$self->feed( @_ );
	$self->_check_all_pdls_present;
	$self->_check_pdl_length;
	autoscale_axis( $_ ) for $self->axes;
}


sub labels
{
	my $self = shift;
	$self->autoscale( @_ );
	my @list = map {
		my $axis = $_;

		if ( defined $axis->{grid} ) {

		     [ map {
			      { range => [ $axis->{grid}->at($_), $axis->{grid}->at($_+1) ] }
		           } 0..$axis->{grid}->nelem -2
		     ]

		}

		else {

			my ( $pdl, $min, $step ) = @{ $axis }{ qw( pdl min step ) };
			[ map {
				{ # anonymous hash
					range => $pdl->type() >= PDL::float()
						? [ $min + $step*$_, $min + $step*($_+1) ]
						: $step > 1
							? [ nhimult( 1, $min + $step*$_ ), nhimult( 1, $min + $step*($_+1) - 1 ) ]
							: $min + $step*$_
				}
			} 0 .. $axis->{n}-1 ];

	    }

	} $self->axes;

	return wantarray ? @list : \@list;
}


sub process
{
	my $self = shift;

	# sanity check
	croak( 'no axes supplied' ) unless @{ $self->axes };
	# default action, when no variables are given, is to produce a histogram
	$self->add_var( name => 'histogram', action => 'Count' ) unless @{ $self->vars };

	#
	$self->autoscale( @_ );

	# process axes
	my $idx = 0;		# flattened bin number
	my @n;			# number of bins in each direction
	# find the last axis and flatten all axes into one dimension, working
	# our way backwards from the last to the first axis
	for my $axis ( reverse $self->axes ) {
	        if ( $log->is_debug ) {
		    $log->debug( 'input (' . $axis->{pdl}->info . ') = ' . $axis->{pdl} );
		    if ( ! defined $axis->{grid} ) {
			$log->debug( "bin with parameters step=$axis->{step}, min=$axis->{min}, n=$axis->{n}" );
		    }
		    else {
			$log->debug( "bin with parameters grid=$axis->{grid}" );
		    }
		}
		croak( 'I cannot bin unless n > 0' ) unless $axis->{n} > 0;
		unshift @n, $axis->{n};			# remember that we are working backwards!
		if ( defined $axis->{grid} ) {
		    $idx = $axis->{pdl}->_flatten_into_grid( $idx, $axis->{grid} );
		}
		else {
		    $idx = $axis->{pdl}->_flatten_into( $idx, $axis->{step}, $axis->{min}, $axis->{n} );
		}
	}
	$log->debug( 'idx (' . $idx->info . ') = ' . $idx ) if $log->is_debug;
	$self->{n} = \@n;

	my $N = reduce { $a * $b } @n; # total number of bins
	croak( 'I need at least one bin' ) unless $N;
	my @vars = map $_->{pdl}, $self->vars;
	$self->{instances} ||= [ map { _make_instance( N => $N, action => $_->{action} ) } $self->vars ];

	#
	{
		local $Data::Dumper::Terse = 1;
		$log->trace( 'process: $self = ' . Dumper $self );
	}

	# now visit all the bins
	my $iter = PDL::NDBin::Iterator->new( bins => \@n, array => \@vars, idx => $idx );
	$log->debug( 'iterator object created: ' . Dumper $iter );
	while( $iter->advance ) {
		my $i = $iter->var;
		$self->{instances}->[ $i ]->process( $iter );
	}

	return $self;
}


sub output
{
	my $self = shift;
	return unless defined wantarray;
	unless( defined $self->{result} ) {
		# reshape output
		my $n = $self->{n};
		my @output = map { $_->result } @{ $self->{instances} };
		for my $pdl ( @output ) { $pdl->reshape( @$n ) }
		if( $log->is_debug ) { $log->debug( 'output: output (' . $_->info . ') = ' . $_ ) for @output }
		$self->{result} = { pairwise { $a->{name} => $b } @{ $self->vars }, @output };
		if( $log->is_debug ) { $log->debug( 'output: result = ' . Dumper $self->{result} ) }
	}
	return wantarray ? %{ $self->{result} } : $self->{result};
}


sub _consume (&\@)
{
	my ( $f, $list ) = @_;
	for my $i ( 0 .. $#$list ) {
		local *_ = \$list->[$i];
		if( not $f->() ) { return splice @$list, 0, $i }
	}
	# If we get here, either the list is empty, or all values in the list
	# meet the condition. In either case, splicing the entire list does
	# what we want.
	return splice @$list;
}

sub _expand_axes
{
	my ( @out, $hash, @num );
	while( @_ ) {
		if( eval { $_[0]->isa('PDL') } ) {
			# a new axis; push the existing one on the output list
			push @out, $hash if $hash;
			$hash = { pdl => shift };
		}
		elsif( ref $_[0] eq 'HASH' ) {
			# the user has supplied a hash directly, which may or
			# may not yet contain a key-value pair pdl => $pdl
			$hash = { } unless $hash;
			push @out, { %$hash, %{ +shift } };
			undef $hash; # do not collapse consecutive hashes into one, too confusing
		}
		elsif( @num = _consume { /^[-+]?(\d+(\.\d*)?|\.\d+)([Ee][-+]?\d+)?$/ } @_ ) {
			croak( 'no axis given' ) unless $hash;
			croak( "too many arguments to axis in `@num'" ) if @num > 3;
			# a series of floating-point numbers
			$hash->{min}  = $num[0] if @num > 0;
			$hash->{max}  = $num[1] if @num > 1;
			$hash->{step} = $num[2] if @num > 2;
		}
		#elsif( @num = ( $_[0] =~ m{^((?:\d+(?:\.\d*)?|\.\d+)(?:[Ee][-+]?\d+)?/)+$}g ) and shift ) {
		#	DOES NOT WORK YET - TODO
		#	print "GMT-style axis spec found! (@num)\n";
		#	croak( 'no axis given' ) unless $hash;
		#	croak( "too many arguments to axis in `@num'" ) if @num > 3;
		#	# a string specification of the form 'min/max/step', a la GMT
		#	$hash->{min}  = $num[0] if @num > 0;
		#	$hash->{max}  = $num[1] if @num > 1;
		#	$hash->{step} = $num[2] if @num > 2;
		#}
		else {
			croak( "while expanding axes: invalid argument at `@_'" );
		}
	}
	push @out, $hash if $hash;
	return @out;
}


sub _random_name { create_uuid( UUID_RANDOM ) }


sub ndbinning
{
	#
	my $binner = __PACKAGE__->new;

	# leading arguments are axes and axis specifications
	#
	# PDL overloads the `eq' and `ne' operators; by checking for a PDL
	# first, we avoid (invalid) comparisons between piddles and strings in
	# the `grep'
	my @leading = _consume { eval { $_->isa('PDL') } || ! $valid_key{ $_ } } @_;

	# consume and process axes
	# axes require three numerical specifications following it
	while( @leading > 3 && eval { $leading[0]->isa('PDL') } && ! grep ref, @leading[ 1 .. 3 ] ) {
		my( $pdl, $step, $min, $n ) = splice @leading, 0, 4;
		$binner->add_axis( name => _random_name, pdl => $pdl, step => $step, min => $min, n => $n );
	}
	if( @leading ) { croak( "error parsing arguments in `@leading'" ) }

	# remaining arguments are key => value pairs
	my $args = { @_ };
	my @invalid_keys = grep ! $valid_key{ $_ }, keys %$args;
	croak( "invalid key(s) @invalid_keys" ) if @invalid_keys;

	# axes
	$args->{axes} ||= [];
	my @axes = @{ $args->{axes} };
	for my $axis ( @axes ) {
		my $pdl = shift @$axis;
		$binner->add_axis( name => _random_name, pdl => $pdl, @$axis );
	}

	# variables
	$args->{vars} ||= [];
	for my $var ( @{ $args->{vars} } ) {
		if( @$var == 2 ) {
			my( $pdl, $action ) = @$var;
			$binner->add_var( name => _random_name, pdl => $pdl, action => $action );
		}
		else { croak( "wrong number of arguments for var: @$var" ) }
	}

	#
	$binner->process;
	my $output = $binner->output;
	my @result = map $output->{ $_->{name} }, @{ $binner->vars };
	return wantarray ? @result : $result[0];
}


sub ndbin
{
	#
	my $binner = __PACKAGE__->new;

	# leading arguments are axes and axis specifications
	#
	# PDL overloads the `eq' and `ne' operators; by checking for a PDL
	# first, we avoid (invalid) comparisons between piddles and strings in
	# the `grep'
	if( my @leading = _consume { eval { $_->isa('PDL') } || ! $valid_key{ $_ } } @_ ) {
		my @axes = _expand_axes( @leading );
		$binner->add_axis( name => _random_name, %$_ ) for @axes;
	}

	# remaining arguments are key => value pairs
	my $args = { @_ };
	my @invalid_keys = grep ! $valid_key{ $_ }, keys %$args;
	croak( "invalid key(s) @invalid_keys" ) if @invalid_keys;

	# axes
	$args->{axes} ||= [];
	my @axes = @{ $args->{axes} };
	for my $axis ( @axes ) {
		my $pdl = shift @$axis;
		$binner->add_axis( name => _random_name, pdl => $pdl, @$axis );
	}

	# variables
	$args->{vars} ||= [];
	for my $var ( @{ $args->{vars} } ) {
		if( @$var == 2 ) {
			my( $pdl, $action ) = @$var;
			$binner->add_var( name => _random_name, pdl => $pdl, action => $action );
		}
		else { croak( "wrong number of arguments for var: @$var" ) }
	}

	$binner->process;
	my $output = $binner->output;
	my @result = map $output->{ $_->{name} }, @{ $binner->vars };
	return wantarray ? @result : $result[0];
}

1;

__END__