#!/usr/bin/perl
#
# Calculate Shannon's entropy to various orders
# for the abacadabb... stream discussed in his paper,
# along with two binary encodings.
#
# Jim Mahoney, April 2005
#

#
# It looks like as the order gets larger 
# we need *lots* of numbers, or H heads to 0.
# I think that for a given H(order), I need
# to keep grabbing more chars until H(order) approaches
# a limit... which is going to be a different approach 
# than the "generate a finite string ahead of time" aproach
# I'm using here.
#
# My guess is that to do the H(o) corrrectly, we need to 
# really get statistics from essentially the whole population
# of substrings of length o.  For the binary strings, that
# means that to do order o I need strings longer that 2**o,
# and for the (abcd) I need strings longer than 4**o.  
#
# This guess looks vaguley consistent with what I see below :
#
#  * with a string length of 2**19, the 4 symbol
#    results are consistent for a while then go south at about H(8) or so,
#
#  * The 0 10 110 111 encoding (true H=1.00) stays right up there
#    until about H(16) or so and then starts to collapse downward
#
#  * the 00 01 10 11 encoding (true H=0.875) is perhaps converging
#    slowly towards this number (the stats on 1st and 2nd bits are
#    different; the program doesn't know this and so has slow convergence),
#    but at around the same point, H(17) or so starts collpasing more
#    quickly downward
# 
# See abcd.png and make_plots; the picture shows that the
# numbers are pretty close to the theoretical values for 
# orders < 15 for the binary encodings.
#

use strict;
use warnings;
use Math::Random::MT qw(srand rand);  # Mersenne Twister pseudo-randoms
#srand(time, time**2, time**3, time**4);  # ? would this work better?
srand(time);

init( a => .5, b => .25, c => .125, d => .125 );
print_theory();

my $stream_original = generate_stream( 2**19 );
analyze($stream_original, 20);

my $coding = {a=>'00', b=>'11', c=>'01', d=>'10'};
print_coding($coding);
my $stream_one = encode($stream_original, $coding );
analyze($stream_one, 20);

$coding = {a=>'0', b=>'10', c=>'110', d=>'111'};
print_coding($coding);
my $stream_two = encode($stream_original, $coding );
analyze($stream_two, 20);

exit;

# === subs ==================================================

# Usage: print_coding( { a=>'00', b=>'01', ... });
sub print_coding {
  my ($coding) = @_;
  my %hash = %$coding;
  my $as_string;
  $as_string .= $_ . "=>'" . $hash{$_} . "', " for sort keys %hash;
  $as_string =~ s/, $//;  # remove trailing comma
  print "\n - - - - using new coding - - - - \n";
  print "   { $as_string } \n";
}

# Usage: $coded = encode($plain, { a=>'00', b=>'11', ...} );
#   $plain = 'aabaccadaccd...'
#   $coded = '000011...';
sub encode {
  my ($input, $coding) = @_;
  my $output;
  for (my $i=0; $i<length($input); $i++){
    $output .= $coding->{ substr($input,$i,1) };
  }
  return $output;
}

# globals for 
#   init( char1=>prob1, char2=>prob2,..);
#   $char = random_symbol();
#   print_theory();
my @symbol;
my @threshold;
my %probability;

# Usage: print_theory();  # only after calling init
sub print_theory {
  print "\n - - - - theory - - - - \n";
  printf(" %6s   %16s   %16s \n",       
	 'symbol', 'probability', 'threshold');
  printf(" %6s   %16.4f   %16.4f \n", 
	 $symbol[$_], $probability{$symbol[$_]}, $threshold[$_]) 
    for 0..$#symbol;
  my $H = 0.0;
  foreach my $p (values %probability){
    $H += $p * log($p);
  }
  $H = - $H/log(2);
  printf("  H exact is $H \n");
}

# Usage: $char = random_symbol();
# Requires %probability and @threshold globals;
# returns a key from the %probability hash consistent with its prob. value.
sub random_symbol {
  my $r = rand();
  for my $i (0..$#symbol){
    return $symbol[$i] if $r <= $threshold[$i];
  }
}

# Usage: init( a => .5, b => .25, c => .125, d => .125 );
# Initializes %probability, @symbol, @threshold globals.
sub init {
  %probability = @_;
  @symbol    = sort keys %probability;
  my $level  = 0;
  @threshold = map {$level+= $probability{$_}} @symbol;
}

# Usage: $str = generate_stream($n_chars);
sub generate_stream {
  my ($n) = @_;
  my $stream = "";
  my $counter = $n;
  $stream .= random_symbol while $counter--;
  return $stream;
}


# Usage: analyze($stream, $max_order);
sub analyze {
  my ($stream, $max_order) = @_;

  my $n = length($stream);                   # use all chars
  $stream .= substr($stream,0,$max_order);   # ... and wrap around

  # ... or, just use a bit less for $n and don't wrap.
  # Doesn't seem to matter much which way
  # my $n = length($stream) - $max_order;

  print "\n - - - - analysis - - - - \n";

  my $log2 = log($n)/log(2);
  print "  $n = 2**$log2 symbol stream :\n  ";
  if ( $n < 70){
    print $stream;
  }
  else {
    print substr($stream,0,30), " ... ", substr($stream,$n-30,30);
  }
  print "\n";

  my $count = [ { "" => $n } ];  # should give 1 when divided by $n

  for (my $order=1; $order<=$max_order; $order++){
    for (my $index=0; $index<$n; $index++){
      my $str = substr($stream,$index,$order);
      # print "order=$order, str=$str\n";
      $count->[$order]{substr($stream,$index,$order)}++;
    }
  }

  #use Data::Dumper;
  #print " -- count -- \n";
  #print Dumper($count);
  #exit;

  #printf(" %7s = %9.6f \n", 'H(0)' , log(scalar @symbol)/log(2) );
  my $n_symbols = scalar(keys %{$count->[1]});
  print_Hn( 0, log($n_symbols)/log(2) );

  for (my $order=1; $order<=$max_order; $order++){
    my $probability_sum = 0;
    my $Hn = 0;  # order'th approx to Shannon's entropy
    #print " order $order :\n";
    #printf(" %16s   %16s \n", 'symbol', 'probability');
    foreach my $substr (sort keys %{$count->[$order]}){
      # substr       = string of length "order" 
      # substr_base  = all but last symbol of substr
      # Hn = sum ( prob(substr) * log_2( prob(substr)/prob(substr_begin)) )
      my $prob_substr = $count->[$order]{$substr}/$n;
      my $substr_base = substr($substr,0,$order-1);
      my $prob_base = $count->[$order-1]{$substr_base}/$n;
      $Hn += $prob_substr * log( $prob_substr / $prob_base );
      $probability_sum += $prob_substr;
      #printf(" %16s  %16.4f \n", $substr, $prob_substr),
    }
    $Hn = - $Hn/log(2);  # convert log base 10 to log base 2
    #printf(" %16s   %16s \n",      '', '------');
    #printf(" %16s   %16.4f \n",    '', $probability_sum);
    print_Hn($order,$Hn);

  }
  
}

sub print_Hn {
  my ($order, $Hn) = @_;
  printf(" %5s = %9.6f \n", 'H('.$order.')' , $Hn );
}

__END__


mahoney@cs entropy$ ./stream.pl 

 - - - - theory - - - - 
 symbol        probability          threshold 
      a             0.5000             0.5000 
      b             0.2500             0.7500 
      c             0.1250             0.8750 
      d             0.1250             1.0000 
  H exact is 1.75 

 - - - - analysis - - - - 
  524288 = 2**19 symbol stream :
  dabababaaaacbbabdbbbbabbdaadab ... cdcccbaaadaaaacabaaaaabaaabaaa
  H(0) =  2.000000 
  H(1) =  1.750659 
  H(2) =  1.750647 
  H(3) =  1.750605 
  H(4) =  1.750424 
  H(5) =  1.749648 
  H(6) =  1.746395 
  H(7) =  1.732622 
  H(8) =  1.682836 
  H(9) =  1.549388 
 H(10) =  1.299425 
 H(11) =  0.964730 
 H(12) =  0.621541 
 H(13) =  0.350193 
 H(14) =  0.174415 
 H(15) =  0.077378 
 H(16) =  0.030893 
 H(17) =  0.011813 
 H(18) =  0.004208 
 H(19) =  0.001433 
 H(20) =  0.000496 

 - - - - using new coding - - - - 
   { a=>'00', b=>'11', c=>'01', d=>'10' } 

 - - - - analysis - - - - 
  1048576 = 2**20 symbol stream :
  100011001100110000000001111100 ... 001100000000001100000011000000
  H(0) =  1.000000 
  H(1) =  0.954646 
  H(2) =  0.915948 
  H(3) =  0.913313 
  H(4) =  0.911194 
  H(5) =  0.909370 
  H(6) =  0.907787 
  H(7) =  0.906234 
  H(8) =  0.904753 
  H(9) =  0.903347 
 H(10) =  0.901965 
 H(11) =  0.900517 
 H(12) =  0.898649 
 H(13) =  0.896111 
 H(14) =  0.892297 
 H(15) =  0.885092 
 H(16) =  0.871824 
 H(17) =  0.847548 
 H(18) =  0.807236 
 H(19) =  0.747887 
 H(20) =  0.670719 

 - - - - using new coding - - - - 
   { a=>'0', b=>'10', c=>'110', d=>'111' } 

 - - - - analysis - - - - 
  917850 = 2**19.8078988743306 symbol stream :
  111010010010000011010100101111 ... 001110000110010000001000010000
  H(0) =  1.000000 
  H(1) =  1.000000 
  H(2) =  1.000000 
  H(3) =  1.000000 
  H(4) =  0.999999 
  H(5) =  0.999987 
  H(6) =  0.999980 
  H(7) =  0.999957 
  H(8) =  0.999906 
  H(9) =  0.999829 
 H(10) =  0.999612 
 H(11) =  0.999242 
 H(12) =  0.998440 
 H(13) =  0.996753 
 H(14) =  0.993489 
 H(15) =  0.986807 
 H(16) =  0.973374 
 H(17) =  0.945961 
 H(18) =  0.883725 
 H(19) =  0.748011 
 H(20) =  0.542006 

