#!/usr/local/bin/perl

## A simple script to illustrate the Asymptotic Equipartition Property (AEP).
## J. Lafferty, 9/22/95

$input = $length = $ARGV[0];

&loadProbs;

if (-e $input) {
   &test($input);
}
elsif ($length > 0) {
   &sample($length,0);
}
else {
   &help;
}

exit;

sub help {
   printf "Usage: $0 <N | filename>\n\n";
   printf "     - This program can be used to get some intuition about the\n";
   printf "       Asymptotic Equipartition Property (AEP):\n\n";
   printf "            -1/N log P(x_1,...x_N) -> H(X)\n\n";
   printf "     - Specifying a numeric length N causes the program to \n";
   printf "       generate a string of N characters, iid ~ X, whose \n";
   printf "       probability is then compared to the entropy H(X).\n";
   printf "       (X is the empirical distribution of characters in a \n";
   printf "       sample of about 20,000 characters of English text.)\n\n";
   printf "     - If a filename is specified, the probability of its \n";
   printf "       contents is calculated and compared to the entropy.\n";
}

sub test {
    local($file) = $_[0];
    local($N) = 0;
    local($logprob) = 0;
    local($outchar) = 0;
    local($i);

    open (INPUT, "$file") || die "Cannot open $file\n";
    printf "Reading file $file\n";
    while(<INPUT>) {
	s/\n/ /;
        tr/[a-z]/[A-Z]/;
	s/[^A-Z ]//g;
	s/[ ]+/ /g;
	$rec = $_;
	for ($i=0; $i<length($rec); ++$i) {
	    $c = substr($rec,$i,1);
            $N++;
	    if ($N <= 5000) {
               printf "$c";
	       if (++$outchar % 70 == 0) {
		   printf "\n";
	       }
	       if ($N==5000) {
		   printf "...\n";
	       }
	    }
	    $logprob -= log($prob{$c});
	}
    }
    close(INPUT);
    $logprob /= ($N*log(2));
    $epsilon = $logprob - $entropy;
    if ($epsilon < 0) {
       $epsilon = -$epsilon;
    }
    printf "\n\nEntropy: %f  -1/%d*log(c_1,...,c_%d): %f  Epsilon: %f\n\n", 
           $entropy, $N, $N, $logprob, $epsilon;
}

sub sample {
    local($N) = $_[0];
    local($s) = $_[1];
    local($logprob) = 0;
    local($i);
    local($outchar) = 0;

    srand(time|($$+$s));
    for ($i=0; $i<$N; ++$i) {
       $c = &randomchar;
       printf "$c";
       if (++$outchar % 70 == 0) {
	   printf "\n";
       }
       $logprob -= log($prob{$c});
    }
    $logprob /= ($N*log(2));
    $epsilon = $logprob - $entropy;
    if ($epsilon < 0) {
       $epsilon = -$epsilon;
    }
    printf "\n\nEntropy: %f  -1/%d*log(c_1,...,c_%d): %f  Epsilon: %f\n\n", 
           $entropy, $N, $N, $logprob, $epsilon;
}

sub randomchar {
    $p = rand;
    $cum = 0;
    foreach $char (keys(%prob)) {
        $cum += $prob{$char};
	if ($p <= $cum) {
           return $char;
       }
    }
    printf "Error: did not find character\n";
    exit -1;
}

sub bycount { $count{$b} <=> $count{$a} }

sub loadProbs {    

### The following probabilities were extracted from a
### a sample of 20,000 characters

    $prob{" "} = 0.176887;
    $prob{"E"} = 0.107721;
    $prob{"T"} = 0.076553;
    $prob{"A"} = 0.065624;
    $prob{"O"} = 0.063095;
    $prob{"N"} = 0.054898;
    $prob{"H"} = 0.051305;
    $prob{"I"} = 0.051204;
    $prob{"S"} = 0.049535;
    $prob{"R"} = 0.046549;
    $prob{"L"} = 0.037543;
    $prob{"D"} = 0.036177;
    $prob{"U"} = 0.024843;
    $prob{"C"} = 0.021504;
    $prob{"W"} = 0.021200;
    $prob{"M"} = 0.020138;
    $prob{"Y"} = 0.017709;	
    $prob{"G"} = 0.015938;
    $prob{"F"} = 0.015331;
    $prob{"B"} = 0.014370;
    $prob{"P"} = 0.013509;
    $prob{"K"} = 0.008500;
    $prob{"V"} = 0.007235;
    $prob{"X"} = 0.000961;
    $prob{"Q"} = 0.000759;
    $prob{"J"} = 0.000708;
    $prob{"Z"} = 0.000202;

    foreach $char (keys(%prob)) {
	$entropy -= $prob{$char} * log($prob{$char});
    }
    $entropy /= log(2);
}
