package Statistics::Test::RandomWalk;

use 5.006;
use strict;
use warnings;

our $VERSION = '0.02';

use Carp qw/croak/;
use Params::Util qw/_POSINT _ARRAY _CODE/;
use Memoize;
use Math::BigFloat;
use Statistics::Test::Sequence;
use Class::XSAccessor {
    constructor => 'new',
    getters => {
        rescale_factor => 'rescale',
    },
    setters => {
        set_rescale_factor => 'rescale',
    },
};

=head1 NAME

Statistics::Test::RandomWalk - Random Walk test for random numbers

=head1 SYNOPSIS

  use Statistics::Test::RandomWalk;
  my $tester = Statistics::Test::RandomWalk->new();
  $tester->set_data( [map {rand()} 1..1000000] );
  
  my $no_bins = 10;
  my ($quant, $got, $expected) = $tester->test($no_bins);
  print $tester->data_to_report($quant, $got, $expected);

=head1 DESCRIPTION

This module implements a Random Walk test of a random number generator
as outlined in Blobel et al (Refer to the SEE ALSO section).

Basically, it tests that the numbers C<[0,1]> generated by a
random number generator are distributed evenly. It divides C<[0,1]>
into C<n> evenly sized bins and calculates the number of expected and
actual random numbers in the bin. (In fact, this counts the
cumulated numbers, but that works the same.)

=head1 METHODS

=head2 new

Creates a new random number tester.

=head2 set_rescale_factor

The default range of the random numbers [0, 1) can be rescaled
by a constant factor. This method is the setter for that factor.

=head2 rescale_factor

Returns the current rescaling factor.

=head2 set_data

Sets the random numbers to operate on. First argument
must be either an array reference to an array of random
numbers or a code reference.

If the first argument is a code reference, the second
argument must be an integer C<n>. The code reference is
called C<n>-times and its return values are used as
random numbers.

The code reference semantics are particularly useful if
you do not want to store all random numbers in memory at
the same time. You can write a subroutine that, for example,
generates and returns batches of 100 random numbers so no
more than 101 of these numbers will be in memory at the
same time. Note that if you return 100 numbers at once and
pass in C<n=50>, you will have a sequence of 5000 random
numbers.

=cut

sub set_data {
    my $self = shift;
    my $data = shift;
    if (_ARRAY($data)) {
        $self->{data} = $data;
        return 1;
    }
    elsif (_CODE($data)) {
        $self->{data} = $data;
        my $n = shift;
        if (not _POSINT($n)) {
            croak("'set_data' needs an integer as second argument if the first argument is a code reference.");
        }
        $self->{n} = $n;
        return 1;
    }
    else {
        croak("Invalid arguments to 'set_data'.");
    }
}

=head2 test

Runs the Random Walk test on the data that was previously set using
C<set_data>.

First argument must be the number of bins.

Returns three array references. First is an array of quantiles.
If the number of bins was ten, this (and all other returned arrays)
will hold ten items.

Second are the determined numbers of random numbers below the
quantiles. Third are the expected counts.

=cut

sub test {
    my $self = shift;
    my $bins = shift;
    if (not _POSINT($bins)) {
        croak("Expecting number of bins as argument to 'test'");
    }
    my $rescale_factor = $self->rescale_factor||1;


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

    if (not defined $data) {
        croak("Set data using 'set_data' first.");
    }

    my $step = 1 / $bins * $rescale_factor;
    my @alpha;
    push @alpha, $_*$step for 1..$bins;

    my @bins = (0) x $bins;
    my $numbers;

    if (_ARRAY($data)) {
        $numbers = @$data;

        foreach my $i (@$data) {
            foreach my $ai (0..$#alpha) {
                if ($i < $alpha[$ai]) {
                    $bins[$_]++ for $ai..$#alpha;
                    last;
                }
            }
        }
    }
    else { # CODE
        my @cache;
        my $calls = $self->{n};
        foreach (1..$calls) {
            # get new data
            push @cache, $data->();
            while (@cache) {
                $numbers++;
                my $this = shift @cache;
                foreach my $ai (0..$#alpha) {
                    if ($this < $alpha[$ai]) {
                        $bins[$_]++ for $ai..$#alpha;
                        last;
                    }
                }
            }
        }
    }

    my @expected_smaller = map Math::BigFloat->new($numbers)*$_/$rescale_factor, @alpha;

    return(
        [map $_/$rescale_factor, @alpha],
        \@bins,
        \@expected_smaller,
    );
}

=head2 data_to_report

From the data returned by the C<test()> method, this
method creates a textual report and returns it as a string.

Do not forget to pass in the data that was returned by C<test()>
or use the C<test_report()> method directly if you do not use
the data otherwise.

=cut

sub data_to_report {
    my $self = shift;
    my $alpha = shift;
    my $got = shift;
    my $expected = shift;
    if (grep {not _ARRAY($_)} ($alpha, $got, $expected)) {
        croak("Please pass the data returned from 'test' to the 'data_to_report' method.");
    }

    my $max_a = _max_length($alpha);
    $max_a = length('Quantile') if length('Quantile') > $max_a;
    my $max_g = _max_length($got);
    $max_g = length('Got') if length('Got') > $max_g;
    my $max_e = _max_length($expected);
    $max_e = length('Expected') if length('Expected') > $max_e;

    my $str = '';

    $str .= sprintf(
        "\%${max_a}s | \%${max_g}s | \%${max_e}s\n",
        qw/Quantile Got Expected/
    );
    $str .= ('-' x (length($str)-1))."\n";
    foreach my $i (0..$#$alpha) {
        $str .= sprintf(
            "\%${max_a}f | \%${max_g}u | \%${max_e}u\n",
            $alpha->[$i], $got->[$i], $expected->[$i]
        );
    }
    return $str;
}

sub _max_length {
    my $max = 0;
    foreach (@{$_[0]}) {
        $max = length $_ if length($_) > $max;
    }
    return $max;

}

=head1 SUBROUTINES

=head2 n_over_k

Computes C<n> over C<k>. Uses Perl's big number support and
returns a L<Math::BigFloat> object.

This sub is memoized.

=cut

memoize('n_over_k');
sub n_over_k {
    my $n = shift;
    my $k = shift;
    my @bits = ((0) x $k, (1) x ($n-$k));
    foreach my $x (1..($n-$k)) {
        $bits[$x-1]--;
    }

    my $o = Math::BigFloat->bone();
    foreach my $i (0..$#bits) {
        $o *= Math::BigFloat->new($i+1)**$bits[$i] if $bits[$i] != 0;
    }

    return $o->ffround(0);
}

1;

__END__

=head1 SEE ALSO

L<Math::BigFloat>, L<Memoize>, L<Params::Util>

Random number generators:
L<Math::Random::MT>, L<Math::Random>, L<Math::Random::OO>,
L<Math::TrulyRandom>, C</dev/random> where available

L<Statistics::Test::Sequence>

The algorithm was taken from: (German)

Blobel, V., and Lohrmann, E. I<Statistische und numerische Methoden
der Datenanalyse>. Stuttgart, Leipzig: Teubner, 1998

=head1 AUTHOR

Steffen Mueller, E<lt>smueller@cpan.orgE<gt>

=head1 COPYRIGHT AND LICENSE

Copyright (C) 2007-2010 by Steffen Mueller

This library is free software; you can redistribute it and/or modify
it under the same terms as Perl itself, either Perl version 5.6 or,
at your option, any later version of Perl 5 you may have available.

=cut