#!/usr/bin/perl

use FileHandle;
STDOUT->format_top_name("TOP");
STDOUT->format_name("LINE");
srand; # seed the random number generator

open(IN, ($datafile = "prob.data"))
               or die('Can\'t open datafile "' . $datafile . "\":\n$!");

while (<IN>) {
    chomp;
    next if /^#/; # skips "commented" lines
    next if !/\S/g; # skips lines consisting only of whitespace
        
    # Keeps track of which experiment we are currently working on.
    $exp_num++;
        
    # Clears out variables from the previous trial.
    $total_inf_count = $total_tpos_count = $total_tneg_count =
        $total_fpos_count = $total_fneg_count = 0;
                
    foreach $item ( ($trials,$pop,$inf_prob,$tpos_prob,$fpos_prob) = split ) {
         $item = eval $item;
    }
        
    for ($j = 0; $j < $trials; $j++) {

        # Clear out variables from values of their previous population trial.
        $inf_count = $tpos_count = $tneg_count = $fpos_count = $fneg_count = 0;

         for ($i = 0; $i < $pop; $i++) {
             if ( $inf_prob > rand ) {
                  $inf_count++;
                  (($tpos_prob > rand) ? $tpos_count : $fneg_count) += 1;
             } else {
                  (($fpos_prob > rand) ? $fpos_count : $tneg_count) += 1;
              }
          }
                  
        # Keep running total of trial-derived information across all trials
        $total_inf_count += $inf_count; $total_tpos_count += $tpos_count;
        $total_tneg_count += $tneg_count; $total_fpos_count += $fpos_count;
        $total_fneg_count += $fneg_count;
                
    }

    # Normalize the total counts based on the number of trials that were run
    $total_inf_count /= $trials;
    $total_tpos_count /= $trials;
    $total_tneg_count /= $trials;
    $total_fpos_count /= $trials;
    $total_fneg_count /= $trials;
        
    # Compute probabilities from the data in the numerical experiment.
    $prob_a = $total_inf_count / $pop;
    $prob_not_a = 1 - $prob_a;
    $prob_b_given_a = $total_tpos_count / $total_inf_count;
    $prob_b_given_not_a = $total_fpos_count / ($pop - $total_inf_count);
    $prob_a_given_b = ($prob_a * $prob_b_given_a) /
                          ($prob_a * $prob_b_given_a +
                            $prob_not_a * $prob_b_given_not_a) ;
                                                        
    # Output the results of this experiment to our formatted report.
        
    write_form_line($exp_num, $trials, $pop,
        $inf_prob, 1 - $inf_prob, $tpos_prob, $fpos_prob,
        ($inf_prob * $tpos_prob) /
        ($inf_prob * $tpos_prob + (1-$inf_prob) * $fpos_prob),
        $prob_a, $prob_not_a, $prob_b_given_a,
        $prob_b_given_not_a, $prob_a_given_b);
}
close(IN);

sub write_form_line {
    my ($experiment, $trial_runs, $population,
        $apA, $apnA, $apBA, $apBnA, $apAB,
        $epA, $epnA, $epBA, $epBnA, $epAB) = @_;
    write;
}

format TOP =
                   Analytic and Experimental Results
                  from our Population Testing Scenario
======================================================================
                            P(A)  P(not-A)  P(B|A)  P(B|not-A)  P(A|B)
======================================================================
.

format LINE =
Experiment #:@>>>>>>>>
$experiment
Trials Run :@>>>>>>>>
$trial_runs
Pop. Size :@>>>>>>>>
$population
Analytic Results :         @.#####  @.#####  @.#####  @.#####  @.#####
$apA, $apnA, $apBA, $apBnA, $apAB
Experimental Results :     @.#####  @.#####  @.#####  @.#####  @.#####
$epA, $epnA, $epBA, $epBnA, $epAB
=======================================================================
.