use strict;
use warnings;
use Math::Histogram qw(make_histogram);
use Test::More;
use File::Spec;
use Data::Dumper;
BEGIN { push @INC, -d "t" ? File::Spec->catdir(qw(t lib)) : "lib"; }
use Math::Histogram::Test;
my @axis_specs3 = ([3, 0., 1.], [3, 0., 1.], [3, 0., 1.]);
my @axis_specs2 = ([3, 0., 1.], [3, 0., 1.]);
my @axis_specs1 = ([3, 0., 1.]);
my $d3 = [
[ [ 1 .. 3 ], [ 4 .. 6 ], [ 7 .. 9 ] ],
[ [ 2 .. 4 ], [ 5 .. 7 ], [ 8 .. 10 ] ],
[ [ 3 .. 5 ], [ 6 .. 8 ], [ 9 .. 11 ] ],
];
my $d2 = [
[ 1, 2, 3 ],
[ 4, 5, 6 ],
[ 7, 8, 9 ],
];
my $d2_out = [
[ 1, 3, 6 ],
[ 5, 12, 21 ],
[ 12, 27, 45 ],
];
my $d1 = [ 1 .. 3 ];
pass();
SCOPE: {
my $h = make_histogram(@axis_specs1);
fill_hist_from_arrays($h, $d1);
my $clone = $h->clone;
$h->cumulate($_) for 0..$h->ndim-1;
# Nuke overflow...
$h->set_bin_content([$h->get_axis(0)->nbins+1], 0);
my $expected_out = dim(3);
cumulate( $d1, $expected_out, [2] );
my $h1 = make_histogram(@axis_specs1);
fill_hist_from_arrays($h1, $expected_out);
ok($h->data_equal_to($h1), "1d cumulated data equal");
my $ary = arrays_from_hist($h);
} # end SCOPE
SCOPE: {
my $h = make_histogram(@axis_specs2);
fill_hist_from_arrays($h, $d2);
my $clone = $h->clone;
$h->cumulate($_) for 0..$h->ndim-1;
# Nuke overflow...
$h->set_bin_content([$_, $h->get_axis(1)->nbins+1], 0) for 0..$h->get_axis(0)->nbins+1;
$h->set_bin_content([$h->get_axis(0)->nbins+1, $_], 0) for 0..$h->get_axis(1)->nbins+1;
#$h->_debug_dump_data;
#print "\n";
#my $expected_out = dim(3,3);
#cumulate( $d2, $expected_out, [2,2] );
my $h2 = make_histogram(@axis_specs2);
fill_hist_from_arrays($h2, $d2_out);
ok($h->data_equal_to($h2));
#my $ary = arrays_from_hist($h);
#use Data::Dumper;
#warn Dumper $ary;
#warn Dumper $expected_out;
#table($ary);
} # end SCOPE
done_testing();
sub fill_hist_from_arrays {
my $h = shift;
my $data = shift;
my $nbins = [map $h->get_axis($_)->nbins, 0..$h->ndim-1];
my $bin_vec = [(0) x $h->ndim];
my $c = [(1) x $h->ndim];
# inefficient as hell, but fine for now.
OUTER: while (1) {
#print "@$c\n";
my $elem = $data;
$elem = $elem->[$_-1] foreach @$c;
$h->fill_bin_w($c, $elem);
foreach my $dim (0..$h->ndim-1) {
$c->[$dim]++;
if ($c->[$dim] > $nbins->[$dim]) {
last OUTER if $dim == $h->ndim-1;
$c->[$dim] = 1;
}
else { last }
}
}
}
sub arrays_from_hist {
my $h = shift;
my $nbins = [map $h->get_axis(0)->nbins, 0..$h->ndim-1];
my $c = [(1) x $h->ndim];
my $data = [];
# inefficient as hell, but fine for now.
OUTER: while (1) {
my $content = $h->get_bin_content($c);
#print $content, "\n";
my $elem = $data;
$elem = ($elem->[$_-1]||=[]) foreach @{$c}[0..$#$c-1];
$elem->[$c->[-1]-1] = $content;
foreach my $dim (0..$h->ndim-1) {
$c->[$dim]++;
if ($c->[$dim] > $nbins->[$dim]) {
last OUTER if $dim == $h->ndim-1;
$c->[$dim] = 1;
}
else { last }
}
}
return $data;
}
############################################
# Alternative cumulation implementation by David Golden
# Note: Appears to yield different (wrong?) results for ndim=2 already?
# given AoA-style multi-dimensional structure, index an
# arbitrary element by successive dereferencing
sub lookup {
my ($ref, $coord) = @_;
$ref = $ref->[$_] for @$coord;
return $ref;
}
# assign to AoA-style multi-dimensional structure
sub assign {
my ($value, $ref, $coord) = @_;
my $orig = $ref;
my $last = pop @$coord;
for my $i ( @$coord ) {
$ref->[$i] ||= [];
$ref = $ref->[$i];
}
$ref->[$last] = $value;
return;
}
# crude matrix dumper
sub table {
my ($data) = @_;
if ( ref($data->[0]) ) {
table($_) for @$data;
print "\n";
}
else {
print "@$data\n";
}
}
# initialize AoA matrix
sub dim {
my ($size, @rest) = @_;
return [ (undef) x $size ] unless @rest;
my $data;
for my $i (0 .. $size-1) {
$data->[$i] = dim(@rest);
}
return $data;
}
# recursively calculate cumulations from maximum element
# (destructive on provided output data matrix)
sub cumulate {
my ($in, $out, $coord) = @_;
# start with bucket count
my $sum = lookup($in, $coord);
# add recursive cumulation from each cardinal direction
for my $i ( 0 .. $#$coord ) {
my $adjacent = [ @$coord ];
next if --$adjacent->[$i] < 0;
my $l = lookup($out, $adjacent);
$sum += defined($l) ? $l : cumulate($in, $out, $adjacent);
}
# update the output cumulation matrix
assign($sum, $out, $coord);
# return the sum for use in recursion
return $sum;
}