#364 type promotion in whistogram is based upon the index, not the weight

[Diab Jerius]

There's an annoying behavior of whistogram with regards to
auto-promotion of argument types which can lead to significant loss of
precision for the unwary.

Type promotion is based upon the first parameter in the function
signature. For whistogram, this is the index, not the weight.
So, if you pass in a double weight, but your index is, say, a long,
you get the default float calculation, not the double calculation that
would be expected based upon the weight's type.

Here's sample code, with a seed chosen to maximize the effect:

use PDL;

srand(372);
my $wt = random( PDL::double, 1000 );
my $lidx = zeroes( PDL::long, 1000 );
my $didx = zeroes( PDL::double, 1000 );

my $lh = whistogram( $lidx, $wt, 1, 0, 1 );
my $dh = whistogram( $didx, $wt, 1, 0, 1 );

my $exp = $wt->dsum;

print "delta, long index: ", $lh - $exp, "\n";
print "delta, double index: ", $dh - $exp, "\n";

And the result:

% perl histo_bug.pl
delta, long index: [0.00069109649]
delta, double index: [0]

Which is ginormous.

Since whistogram's public API can't be changed, I think a proper fix
would be to turn it into a Perl wrapper for an internal routine with
weight the first argument, so that type promotion works as expected.

[Chris Marshall]

Why can't the bug be fixed?

[Chris Marshall]

Printing out the value of $exp you can see that the relative error is very small:

Weight sum: 505.118205387888
delta, long index: [0.00069109649]
delta, double index: [0]

A patch for the modified behavior that is desired would be welcome.

[mohawk]

This matter will now be tracked at https://github.com/PDLPorters/pdl/issues/45

[mohawk]

The input was never an "index", it was just input data. The operation's "type-promotion" (really, determination) was never based on only the first parameter. The output's type is written "float+", which is what happened here (the output got promoted to float), which can be seen easily by using this updated code:

use PDL;

srand(372);
my $wt = random( PDL::double, 1000 );
my $lidx = zeroes( PDL::long, 1000 );
my $fidx = zeroes( PDL::float, 1000 );
my $didx = zeroes( PDL::double, 1000 );

my $lh = whistogram( $lidx, $wt, 1, 0, 1 );
my $fh = whistogram( $fidx, $wt, 1, 0, 1 );
my $dh = whistogram( $didx, $wt, 1, 0, 1 );

my $exp = $wt->dsum;
print "$_->[0]=", $_->[1]->info, "\n" for ['lh',$lh], ['fh',$fh], ['dh',$dh], ['exp',$exp];

print "delta, long index: ", $lh - $exp, "\n";
print "delta, float index: ", $fh - $exp, "\n";
print "delta, double index: ", $dh - $exp, "\n";

You get an identical difference between a "long" input type, and a "float" input type. The difference between those and double is simply the result of the difference between single and double precision. There is no bug here.

[Diab Jerius]

Sorry for the late reply, but I believe that there's still an issue here.

Here's the output of your script
(With a slight fix on PDL 2.038 setting $exp = PDL($wt->dsum), as dsum
doesn't return a piddle):

lh=PDL: Float D [1]
fh=PDL: Float D [1]
dh=PDL: Double D [1]
exp=PDL: Double D []
delta, long index: [0.00069109649]
delta, float index: [0.00069109649]
delta, double index: [0]

As you point out, the difference is due to the difference in precision
between float
and double, but that's the actual problem. In each case the weight
parameter is double
precision, but the determined type of the output histogram does not take
that into account,
leading to a catastrophic loss of precision.

You state "the operation's "type-promotion" (really, determination) was
never based on only the first parameter", but that seems to be exactly
what is happening. There are only two parameters of import,
the input data and the weight, and it doesn't seem to pay attention to
the latter.

[mohawk]

I have posted your comment on https://github.com/PDLPorters/pdl/issues/45 and will deal with it there. This bug-tracker is not being used anymore.