#!/usr/bin/perl
use warnings;
use strict;

=head1 NAME

 sim_polls.pl - This program attempts to simulate polling results, as 
 seen in the national news media.

=head1 SYNOPSIS

 sim_polls.pl [delta_popularity] [standard_deviation_of_polls]
 Options:
        --help          brief help message
        --man           full documentation
	--verbose	wordy explanation of the script added.

=head1 VERSION

 author   David Myers
 date     8/24/2008 
 modified 8/24/2008

=head1 COPYRIGHT
 
 This code is copyrighted (C) 2008 by David W. Myers and is licensed
 using the Artistic license, version 2.0.
 A copy of the Artistic license can be found here:
 http://www.opensource.org/licenses/artistic-license-2.0.php

=head1 DESCRIPTION
 
 sim_polls.pl - This program attempts to simulate polling results, as 
 seen in the national news media.

=cut

use Getopt::Long;
use Pod::Usage;
# use Curses;

my $help    = 0;
my $man     = 0;
my $verbose = 0;

GetOptions(
    'help|?' => \$help,
    man      => \$man,
    verbose  => \$verbose,
) or pod2usage(2);

pod2usage( -exitval => 0, -verbose => 1 ) if $help;
pod2usage( -exitval => 0, -verbose => 2 ) if $man;

my $delta_popularity = shift || 2.0;
my $deviation        = shift || 2.0;

#
# number of simulations
#
my $iterations = 10000;
my $max_trials = 70;

#
# printed explanation
#
if ($verbose) {
    print "  This program does a Monte Carlo simulation of the polling of\n";
    print "the current election to determine if the media gabbing about\n";
    print "the polls is in any way representative of the data presented.\n";
    print "Currently - 8-24-2008 - using polls listed on pollster.com,\n";
    print "Barack Obama is leading by a score of 62-3-5 over the last 70\n";
    print "polls. This is a ratio of 90.71%.\n\n";
    print "  Despite this, people continue to insist that the race is dead\n";
    print "even. This program simulates the distribution of ";
    print "$iterations occurrences\n";
    print "of $max_trials polls by assuming that Obama has a ";
    print "$delta_popularity percent lead\n";
    print "with a $deviation percent standard deviation. These popularity\n";
    print "profiles are assumed to be normally distributed.\n\n\n";
}

#
# statistics counts;
#
my $greater_total = 0;
my $lesser_total  = 0;
my $max_greater   = 0;
my $max_lesser    = 0;

#
# monte carlo run
#
for ( 1 .. $iterations ) {
    my $run_greater = 0;
    my $run_lesser  = 0;
    for my $trial ( 1 .. $max_trials ) {
        my @rand    = gaussian_rand();
        my $greater = $rand[0] * $deviation + $delta_popularity;
        my $lesser  = $rand[1] * $deviation;
        if ( $greater > $lesser ) {
            $run_greater++;
        }
        else {
            $run_lesser++;
        }
    }
    $greater_total += $run_greater;
    $lesser_total  += $run_lesser;
    $max_greater = $run_greater if ( $run_greater > $max_greater );
    $max_lesser  = $run_lesser  if ( $run_lesser > $max_lesser );
}

my $greater_pct = 100.0*$greater_total / ( $iterations * $max_trials );
my $lesser_pct  = 100.0*$lesser_total /  ( $iterations * $max_trials );

print "Greater total: $greater_total\n";
print "Lesser  total: $lesser_total\n";
print "Greater pct  : $greater_pct\n";
print "Lesser  pct  : $lesser_pct\n";
print "Greater max  : $max_greater\n";
print "Lesser max   : $max_lesser\n";

#
# Gaussian rand is a routine taken from the Perl Cookbook
# and seems to have been derived from a routine from
# Numerical Recipes by Press et al.
#

sub gaussian_rand {
    my ( $u1, $u2 );    # uniformly distributed random numbers
    my $w;              # variance, then a weight
    my ( $g1, $g2 );    # gaussian-distributed numbers

    do {
        $u1 = 2 * rand() - 1;
        $u2 = 2 * rand() - 1;
        $w  = $u1 * $u1 + $u2 * $u2;
    } while ( $w >= 1 );

    $w  = sqrt( ( -2 * log($w) ) / $w );
    $g2 = $u1 * $w;
    $g1 = $u2 * $w;

    # return both if wanted, else just one
    return wantarray ? ( $g1, $g2 ) : $g1;
}