#!/usr/bin/env perl

use strict;
use warnings;

use  Math::Random::MT;
use Crypt::Random;
use Data::Dumper;

my $MTgen = Math::Random::MT->new();

# number of repeats
my $NUM_REPEATS ||= 2000000;

# range of rand ints to produce is 0..$MAX_RAND INCL!
my $MAX_RAND ||= 30;

# for LAG=1 it calculates the prob between two consecutive random numbers
# for LAG=2, it checks when there is one-in-between the two rands. e.g. A x B
# e.g. LAG=1 rand[0] -> rand[1]
# LAG=3 rand[0] -> rand[3] (rand [1] and [2] are skipped)
my $LAG ||= 1;

sub     perls_rand { int(rand($MAX_RAND+1)) }
sub     mt_rand { int($MTgen->rand($MAX_RAND+1)) }
# size: number of bits in base 10 i.e. a base10 digit
# uniform for uniform rn
# do not use these two because it will block after a few rand requests
sub     dev_rand { Crypt::Random::makerandom(Size=>1, Uniform=>1) }
sub     dev_u_rand { Crypt::Random::makerandom(Size=>1, Uniform=>1, Strength=>0) }

my %results = ();
my $best_results = [ undef, undef ];
foreach( 
        (["transition matrix test : Perl's rand()", \&perls_rand, $NUM_REPEATS, $MAX_RAND, $LAG],
         ["transition matrix test : Math::Random::MT's rand()", \&mt_rand, $NUM_REPEATS, $MAX_RAND, $LAG],
# do not use this because it will block after a few rand requests:
#        ["transition matrix test : /dev/random's rand()", \&dev_rand, $NUM_REPEATS, $MAX_RAND, $LAG],
#        ["transition matrix test : /dev/urandom's rand()", \&dev_u_rand, $NUM_REPEATS, $MAX_RAND, $LAG],
        )
){
        my $aresult = do_a_test(@$_);
        $results{$aresult->[1]} = $aresult;
        print $aresult->[0]; # the report
        print "====================================\n\n\n";
        if( ! defined($best_results->[0]) or $best_results->[0] > $aresult->[2] ){
                $best_results->[0] = $aresult->[2]; # mean deviation from expected count
                $best_results->[1] = $_->[0]; # the label
        }
}
print "\n\n$0 : done, here is a summary:\n";
foreach(values %results){
        print "$0 : ".$_->[1]." with mean %diff from expected count to be ".sprintf("%.3f",$_->[2])." %\n";
}
print "\n\n$0 : best results for ".$best_results->[1]." with mean %diff from expected count to be ".sprintf("%.3f",$best_results->[0])." %\n";
print "$0 : end.\n";
exit(0);

sub     do_a_test {
        my ($test_label, $rand_sub, $nrepeats, $maxrand, $lag) = @_;
        my @mat = ( map { [(0) x ($maxrand+1)] } (0..$maxrand) );

        my $expected = $nrepeats / ($maxrand+1)**2;
        my $expected_probabilities = $expected / $nrepeats;
        my @allrands = (0..$maxrand);
        my ($i, $j);
        my $lag1 = $lag > 1 ? $lag-1 : 0;
        print "do_a_test() : $test_label :";
        for(1..$nrepeats){
                $i = $rand_sub->();
                for(1..$lag1){ $rand_sub->() }
                $j = $rand_sub->();
                $mat[$i][$j]++;
                if( $_ % 100000 == 0 ){ print " ".$_; select()->flush(); }
        }

        my $deviations = 0;
        foreach $i (@allrands){
                foreach $j (@allrands){
                        $mat[$i][$j] /= $nrepeats; # convert to probabilities
                        $deviations += percentage_difference($expected_probabilities, $mat[$i][$j]);
                }
        }
        $deviations /= ($maxrand+1)**2;

        my $report = "$test_label\n           after $nrepeats repeats with a lag of $lag and using random integers from 0 to $maxrand:\n";
        $report .= "              ".join("          ", @allrands)."\n";
        foreach $i (@allrands){
                $report .= sprintf("%-5d", $i);
                foreach $j (@allrands){
                        $report .= sprintf("%10.7f ", $mat[$i][$j]);
                }
                $report .= "\n";
        }
        $report .= "\nmean %diff from expected count ($expected) / expected transition probability ($expected_probabilities) is ".sprintf("%.5f", $deviations)." % (instead of 0).\n";

        return [$report,
                $test_label,
                $deviations
        ]
}
sub     percentage_difference {
        my ($a, $b) = @_;
        return 200 * abs($a-$b) / ($a+$b)
}
sub     percentage_change {
        my ($a, $b) = @_;
        return 100 - 100 * $a / $b
}