File Coverage

lib/Statistics/Sampler/Multinomial/AliasMethod.pm

Criterion	Covered	Total	%
statement	92	92	100.0
branch	17	22	77.2
condition	6	8	75.0
subroutine	15	15	100.0
pod	2	2	100.0
total	132	139	94.9

package Statistics::Sampler::Multinomial::AliasMethod;

use 5.014;
use warnings;
use strict;

our $VERSION = '0.7';

use Carp;
use Ref::Util qw /is_arrayref/;
use List::Util qw /min sum/;
use List::MoreUtils qw /first_index/;
use Scalar::Util qw /blessed/;

use parent qw /Statistics::Sampler::Multinomial/;

sub _validate_prng_object {
    my ($self, $prng) = @_;

    #  Math::Random::MT::Auto has boolean op overloading
    #  so make sure we don't trigger it or our tests fail
    #  i.e. don't use "if $prng" directly
    #  (and we waste a random number, but that's less of an issue)
    return 1 if !defined $prng;

    croak 'prng arg is not an object'
      if not blessed $prng;
    croak 'prng arg does not have rand() method'
      if not $prng->can('rand');

    return 1;
}

sub _initialise_alias_tables {
    my ($self, %args) = @_;

    my $probs = $self->{data};

    #  caller has not promised they sum to 1
    if (!$self->{data_sum_to_one}) {
        my $sum = sum (@$probs);
        if ($sum != 1) {
            my @scaled_probs = map {$_ / $sum} @$probs;
            $probs = \@scaled_probs;
        }
    }

    #  algorithm and comments stolen/adapted from
    #  https://hips.seas.harvard.edu/blog/2013/03/03/the-alias-method-efficient-sampling-with-many-discrete-outcomes/

    my (@smaller, @larger);
    my @J = (0) x scalar @$probs;
    my @q = (0) x scalar @$probs;
    my $kk = -1;
    my $K = scalar @$probs;

    foreach my $prob (@$probs){
        $kk++;
        $q[$kk] = $K * $prob;
        if ($q[$kk] < 1.0) {
            push @smaller, $kk
        }
        else {
            push @larger, $kk;
        }
    }

    # Loop though and create little binary mixtures that
    # appropriately allocate the larger outcomes over the
    # overall uniform mixture.
    while (scalar @smaller && scalar @larger) {
        my $small = pop @smaller;
        my $large = pop @larger;

        $J[$small] = $large;
        $q[$large] = ($q[$large] + $q[$small]) - 1;

        if ($q[$large] < 1.0) {
            push @smaller, $large;
        }
        else {
            push @larger, $large;
        }
    }
    # handle numeric stability issues
    # courtesy http://www.keithschwarz.com/darts-dice-coins/
    while (scalar @larger) {
        my $g  = shift @larger;
        $q[$g] = 1;
    }
    while (scalar @smaller) {
        my $l  = shift @smaller;
        $q[$l] = 1;
    }

    #  need better names for these,
    $self->{J} = \@J;
    $self->{q} = \@q;

    return if !defined wantarray;

    #  should not expose these to the caller
    my %results = (J => \@J, q => \@q);
    return wantarray ? %results : \%results;
}

sub draw {
    my ($self, $args) = @_;

    my $prng = $self->{prng};

    my $q  = $self->{q}
      // do {$self->_initialise_alias_tables; $self->{q}};

    my $J  = $self->{J};
    my $K  = scalar @$J;
    my $kk = int ($prng->rand * $K);

    # Draw from the binary mixture, either keeping the
    # small one, or choosing the associated larger one.
    return ($prng->rand < $q->[$kk]) ? $kk : $J->[$kk];
}

sub draw_n_samples {
    my ($self, $n) = @_;

    my $prng = $self->{prng};

    my $q  = $self->{q}
      // do {$self->_initialise_alias_tables; $self->{q}};
    my $J  = $self->{J};
    my $K  = scalar @$J;

    my @draws = (0) x $K;
    for (1..$n) {
        my $kk = int ($prng->rand * $K);
        # Draw from the binary mixture, either keeping the
        # small one, or choosing the associated larger one.
        # {SWL: could try to use Data::Alias or refaliasing here
        # as the derefs cause overhead, albeit the big overhead
        # is the method calls}
        #push @draws, ($prng->rand < $q->[$kk]) ? $kk : $J->[$kk];
        $draws[($prng->rand < $q->[$kk]) ? $kk : $J->[$kk]]++;
    }

    return \@draws;
}

#  Cuckoo package to act as a method wrapper
#  to use the perl PRNG stream by default.
package Statistics::Sampler::Multinomial::AliasMethod::DefaultPRNG {
    sub new {
        return bless {}, __PACKAGE__;
    }
    sub rand {
        rand();
    }
    1;
}

1;
__END__