#!/usr/bin/perl -w

# Another look at letter frequencies and 
# shannon's entropy in darwin's origin of the species.
#  - Jim Mahoney, Dec 2001

my $start=time();
my $log2 = log(2);
my (%firstorder, %secondorder);

my @allchars = qw(a b c d e f g h i j k l m n o p q r s t u v w x y z);
push @allchars, ' ';

foreach my $char (@allchars) {
  $firstorder{$char} = 0;
  foreach my $nextchar (@allchars) {
    $secondorder{$char}{$nextchar} = 0;
  }
}

open(DARWIN, 'darwin.txt') or die(" ** OOPS **  couldn't open input");
print "Analyzing 'darwin.txt' for 2nd order character pair counts ... \n";

while (my $line = <DARWIN>) {
  next if $line =~ /^\s*$/;
  $line = lc($line);
  for ($line) {
    s/[0-9.\-',;&\/()?!:]//g;                                #']; #
    s/\s+/ /g;
    s/\s*$/ /;   # make each line end with a space ( i.e. word separator)
  }

  my @chars = split //, $line;
  $secondorder{' '}{$chars[0]}++;   # as if first char came after space 
  while (@chars) {
    my $char = shift @chars;
    my $nextchar = $chars[0];
    $firstorder{$char}++;
    $secondorder{$char}{$nextchar}++ if $nextchar;
  }
}

my ($total, $totalpairs);
foreach my $char (keys %firstorder) {
  $total += $firstorder{$char};
  foreach my $nextchar (keys %{$secondorder{$char}}) {
    $totalpairs += $secondorder{$char}{$nextchar};
  }
}

#    definition of P(b|a) =  probability of b given a
#      P(ab) = P(a)*P(b|a)  ;   prob(a and b) = prob(a) * prob(b given a)


print " char |  count  | probability \n";
print " ---- +  -----  + ----------- \n";
my $entropy;
foreach my $char (keys %firstorder) {
  my $p = $firstorder{$char}/$total;
  $entropy -=  $p * log($p)/$log2 unless $p==0;
  printf( " '$char'  | %7.0f | %8.6f \n", $firstorder{$char},$p);
}

my $maxentropy = log(27.0)/$log2;
my $redundancy = 1.0 - $entropy/$maxentropy;
print "Shannon's 1st order entropy = -sum{p_i log2(p_i)} = $entropy .\n";

print "    (With 27 symbols (a-z, space), \n";
print "      max entropy = log base 2 (27) = $maxentropy .)\n";
printf("So character redundancy (first order) = 1 - entropy/max = %5.2f %% \n", 100*$redundancy);

print "  (But since we're using 8 bits = 256 distinct symbols to encode \n";
print "   these 27 symbols, the file itself has of course a much higher \n";
print "   redundancy factor.) \n\n";

print " ---- Now on to second order ------- \n\n";
print "      Total number of pairs = $totalpairs . \n";
print "      P(a|b) for pairs ab = P(ab)/P(a),  where a=row, b=column : \n";

print "     |";
foreach my $char (@allchars) {
  print " '$char'  ";
}
print "\n";

print " --- +";
foreach my $char (@allchars) {
  print " ---- ";
}
print "\n";

my $entropysecond;
my $tdivp = $total/$totalpairs;
foreach my $char (@allchars) {
  print " '$char' |";
  foreach my $nextchar (@allchars){
    #      P(ab) = N(ab)/totalpairs;  P(a)=N(a)/total
    #      P(b|a) = P(ab)/P(a)
    my $Pab = $secondorder{$char}{$nextchar}/$totalpairs;
    my $Pa = $firstorder{$char}/$total;
    my $Pbgivena = $Pab/$Pa;
    $entropysecond -= $Pab * log($Pbgivena)/$log2 unless $Pbgivena == 0;
    #printf(" %5.0f", $secondorder{$char}{$nextchar});
    printf(" %5.3f", $Pbgivena);
  }
  print "\n";
}
print "\n";

print "Second order entropy is " . $entropysecond . " . \n";
my $redund2 = 1.0 - $entropysecond/$maxentropy;
printf("Symbol redundancy (2nd order) =  %5.2f %% \n\n", 100*$redund2);

print "Elapsed time was " . (time() - $start) . " seconds. \n";


