#! /usr/local/bin/perl
#
if(grep(/^\-(h|\-help)$/i, @ARGV))
{
print STDOUT << "ENDOFHELP";
Read N-gram tables and caluclate the probability
of a (partial) sentence.

use:
ngramprobability [--verbose] [--count <k>] [--n-gram <n>] [--help] \
                 <N-gramtable> [<N-gramtable...>] "<sentence>"

--verbose
-v
Print all contributions

--perplexity
-p
Print Log Perplexity (Perplexity) over the input with respect to the
language model (i.e., the cross-entropy). Surpresses normal output.

--count <k>
-c <k>
Maximum count for Good-Turing smoothing.  Actually, this is
a variant of Good-Turing discounting: Katz smoothing.
Default -1, i.e., NO SMOOTHING. Suggested value: --count = 5

--n-gram <n>
-n <n>
The N-gram used. If n is larger than the depth of the
table, the highest possible value is used.

--license
-l
Print license information

--help
-h
This message

N-gramtable is a list of files (or - for STDIN) with whitespace 
delimited lines of N-gram counts and N-grams with the lowercase 
words in reverse order. Lines starting with '#' are removed.
From an N-gram table (N<= 4) all smaller Ngrams are reconstructed.
You can enter several tables, even mixing different N-grams.
For example:
     435	not	do	i

sentence is a list of words enclosed in quotes. It MUST be the
following the tables. All non-alphabetic characters are removed.
For the start of the sentence, the largest possible N-gram is used,
i.e., unigrams for the first word, bigrams for the second etc.

ENDOFHELP
exit;
}
#
#
###############################################################################
if(grep(/^\-(l|\-license)$/i, @ARGV))
{
print STDOUT << "ENDOFLICENSE" ;

Copyright R.J.J.H. van Son © 20002

Author Rob van Son
Institute of Phonetic Sciences & ACLC
University of Amsterdam
Herengracht 338
NL-1016CG Amsterdam, The Netherlands
Email: Rob.van.Son\@hum.uva.nl
	   rob.van.son\@workmail.com
WWW  : http://www.fon.hum.uva.nl/rob/
mail:  Institute of Phonetic Sciences
	   University of Amsterdam
	   Herengracht 338
	   NL-1016CG Amsterdam
	   The Netherlands
	   tel +31 205252183
	   fax +31 205252197

License for use and disclaimers

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.


ENDOFLICENSE
exit;
};
#######################################################
#
# Get Smoothing factors (Good-Turing)
# Default is NO smoothing
my $k=-1;
my $i;
for($i=0; $i<scalar(@ARGV); ++$i)
{
	if($ARGV[$i] =~ /^\-(c|\-count)$/i)
	{
		$k = $ARGV[$i+1];
		# Remove option
		splice(@ARGV, $i, 2);
	};
};

# Get N-gram power
# Default maximum
my $n=0;
for($i=0; $i<scalar(@ARGV); ++$i)
{
	if($ARGV[$i] =~ /^\-(n|\-n\-gram)$/i)
	{
		$n = $ARGV[$i+1];
		# Remove option
		splice(@ARGV, $i, 2);
	};
};

# Print all? Verbose
my $verbose = scalar grep(/^\-(v|\-verbose)$/i, @ARGV);
@ARGV = grep(!/^\-(v|\-verbose)$/i, @ARGV);

# Check all? Verbose
my $CheckAll = scalar grep(/^\-\-CheckAll$/i, @ARGV);
@ARGV = grep(!/^\-\-CheckAll$/i, @ARGV);

# Perplexity only? 
my $Perplexity = scalar grep(/^\-(p|\-perplexity)$/i, @ARGV);
@ARGV = grep(!/^\-(p|\-perplexity)$/i, @ARGV);

# Back-off?
my $backoff = scalar grep(/^\-(b|\-backoff)$/i, @ARGV);
@ARGV = grep(!/^\-(b|\-backoff)$/i, @ARGV);
# No backoff without smoothing
$k = 5 if $k <= 0 && $backoff;

#########################################################
#
# Start handling input
#
#########################################################
# Get Sentence(s)
my $Sentence = pop;
# Use a file if Sentence exists
my $SentenceInput = "";
if(-s $Sentence || $Sentence eq '-' || $Sentence eq 'STDIN'
    || $Sentence eq 'stdin')
{
	$SentenceInput = $Sentence;
	$Sentence = "";
};

# Get table
my @NgramtableList = @ARGV;

# Define Ngramtables
my @NgramTables = ({});
$NgramTables[0]->{0} = 0;

# Read in the N-gram tables
my $NgramInputTable;
foreach $NgramInputTable (@NgramtableList)
{
	# Allow standard input
	$NgramInputTable = "STDIN" if $NgramInputTable eq '-';
	
	# Open table or stop with error message
	open(INPUTTABLE, "<$NgramInputTable") || die "$NgramInputTable: $!\n" unless $NgramInputTable eq 'STDIN';
	open(INPUTTABLE, "<&STDIN") || die "$NgramInputTable: $!\n" if $NgramInputTable eq 'STDIN';

	# Read table
	while(<INPUTTABLE>)
	{
		# Skip empty lines
		next unless /\S/;
		# Skip comments
		next if /^\s*\#/;

		# Split input (input stored in $_, split on whitespace: split(' ', $_))
		my ($Count, @W) = split;

		# Store Count, if all relevant words are present
		$n = scalar(@W) unless $n > 0;
		my $NgramPower = $n < scalar(@W) ? $n : scalar(@W);
		my $i;
		for($i=$NgramPower; $i>=0; --$i)
		{
			$NgramTables[$i] = {} unless ref($NgramTables[$i]) eq 'HASH';
			
			$Key = $i ? join(" ", reverse(@W[0..$i-1])) : 'Total';
			$NgramTables[$i]->{$Key} += $Count;
		};
	};
};

# Calculate the maximum number of N-grams
# and the class size distribution of counts Nc
my @TotalNgrams = ();
my @MaxNgrams = ();
my @NgramNc = ();
for($i=0; $i<=$n; ++$i)
{
	$TotalNgrams[$i] = scalar(keys(%{$NgramTables[$i]}));
	$MaxNgrams[$i] = $i ? $TotalNgrams[1]**$i : 1;
	$NgramNc[$i][0] = ($MaxNgrams[$i] - $TotalNgrams[$i]);
};

# Distribution of counts from N1-Nk
# Iterate over Ngram depth
for($i=0;$i<=$n; ++$i)
{
	my $NgramEntry;
	foreach $NgramEntry (keys(%{$NgramTables[$i]}))
	{
		my $count = $NgramTables[$i]->{$NgramEntry};
		if($count && $count <= $k + 1)
		{
			++$NgramNc[$i][$count];
		};
	};
};

# Good-Turing discounting. 0* = N1 / N0
# Actually, this is Katz smoothing, a variant of Good-Turing
# discounting
sub GoodTuringCount	# ($count, $Nc, $Nc1, $k, $N1, $Nk1) -> $newcount
{
	my $count = shift;
	my $Nc = shift;
	my $Nc1 = shift; 
	my $k = shift;
	my $N1 = shift;
	my $Nk1 = shift;
	
	# Check for existence of counts
	return 0 unless($Nc && $N1 && (1 - ($k+1)*$Nk1/$N1));
	
	return $count ? (($count+1)*$Nc1/$Nc - $count*($k+1)*$Nk1/$N1) 
			/ (1 - ($k+1)*$Nk1/$N1) : $Nc1 / $Nc;
}

# Smoothed N-gram counts
my @SmoothedNgramCounts = ();
for($j=0; $j<=$k; ++$j)
{
	for($i=0;$i<=$n; ++$i)
	{
		# Smooth counts when there is a need for it, i.e., $NgramNc[$i][0] > 0
		$SmoothedNgramCounts[$i][$j] = $NgramNc[$i][0] ?
					GoodTuringCount($j, $NgramNc[$i][$j], $NgramNc[$i][$j+1], 
									$k, $NgramNc[$i][1], $NgramNc[$i][$k+1])
					: $j;
	};
};

# Store smoothed counts
for($i=0;$i<=$n; ++$i)
{
	foreach $NgramEntry (keys(%{$NgramTables[$i]}))
	{
		my $count = $NgramTables[$i]->{$NgramEntry};
		if($count <= $k)
		{
			$NgramTables[$i]->{$NgramEntry} = $SmoothedNgramCounts[$i][$count];
		};
	};
};

# Check smoothed counts
for($i=0;$i<=$n && $CheckAll; ++$i)
{
	my $sumcount = 0;
	foreach $NgramEntry (keys(%{$NgramTables[$i]}))
	{
		$sumcount += $NgramTables[$i]->{$NgramEntry};
	};
	print "Total $i: $NgramTables[0]->{'Total'} = ", $sumcount + $NgramNc[$i][0]*$SmoothedNgramCounts[$i][0], 
	" = $sumcount + $NgramNc[$i][0]*$SmoothedNgramCounts[$i][0] (Counts* + 0 counts)\n";
};

#############################################################
#
# Ready with the preparations. The language model has been constructed,
#  now we are going to use it.
#
#############################################################

# Handle file input
if($SentenceInput)
{
	if($SentenceInput eq '-' || $SentenceInput eq 'STDIN'  || $SentenceInput eq 'stdin')
	{
		open(TEXTINPUT, "<&STDIN") || die "STDIN $SentenceInput: $!\n";
	}
	else
	{
		open(TEXTINPUT, "<$SentenceInput") || die "$SentenceInput: $!\n";
	};
};

# Iterate over input
my @AlphaContext = ();
my @NNonZeroDenominatorCount = ();
my @NDenominatorCount = ();
my $SumLogProb = 0;
my $Nwords = 0;
my $LogProbability = 0;
my $CurrentNgramWords = "";
while($Sentence || ($Sentence = <TEXTINPUT>))
{
	# Skip empty lines and comments
	next unless $Sentence =~ /\S/;
	next if $Sentence =~ /\s*\#/;
	
	# Clean up sentence (lowercase, only alphabetic characters)
	$Sentence = lc($Sentence);
	
	# Protect word internal '
	s/([a-z])\'([a-z])/\1__\2/ig;
	
	# Remove non-sentence ending periods
	s/\b(dr|mr|mrs|drv|blvd)\./\1/ig;
	
	# Restore word internal '
	s/([a-z])\_\_([a-z])/\1\'\2/ig;
	
	# Remove unwanted characters
	$Sentence =~ s/[^a-z\']+/ /ig;
	
	# Add enough sentence start symbols (.) and a sentence END symbol
	my @WordSequence = split(' ', (". "x($n-1)).$Sentence.' .');
	
	# Clean out $Sentence
	$Sentence = "";
	
	# Calculate probabilities of the words
	$LogProbability = 0;
	$CurrentNgramWords = "";
	# Start at first real word
	for($j=$n-1; $j<scalar(@WordSequence); ++$j)
	{
    	# N-gram count
		my $CurrentCount = 0;
		my $CurrentDenominatorCount = 0;
		my $CurrentNonZeroDenominatorCount = 0;
		my $Alpha = 1;

		# Unnecessary precautions
		if($j >= $n-1 && $TotalNgrams[$n])
		{
			my $CurrentNgramWords = join(" ", @WordSequence[($j-($n-1))..$j]);
			my $NonZeroCount = $NgramTables[$n]->{$CurrentNgramWords};
			my $smoothedCount =  ($NonZeroCount || $SmoothedNgramCounts[$n][0]);
			$CurrentCount = $backoff ? $NonZeroCount : $smoothedCount; 

			# Calculate Exact normalization 
			$CurrentDenominatorCount = 0;
			my $CurrentUnigram;
			# Iterate over ALL smoothed Ngrams
			my $N1context = join(" ", @WordSequence[($j-($n-1))..($j-1)]);
			if(!$NDenominatorCount[$n]->{$N1context} > 0 || !$NNonZeroDenominatorCount[$n]->{$N1context} > 0)
			{
				foreach $CurrentUnigram (keys(%{$NgramTables[1]}))
				{
					my $CurrentN1gramWords = "$N1context $CurrentUnigram";
					my $NonZeroCount = $NgramTables[$n]->{$CurrentN1gramWords};

					# Unless Backoff, use smoothed Zero counts
					$CurrentDenominatorCount += $NonZeroCount || $SmoothedNgramCounts[$n][0];
					$CurrentNonZeroDenominatorCount += $NonZeroCount;
				};
			};
			$NDenominatorCount[$n]->{$N1context} = $CurrentDenominatorCount;
			$NNonZeroDenominatorCount[$n]->{$N1context} = $CurrentNonZeroDenominatorCount;

			# Backoff if necessary
			# For backoff, calculate the sum over all words of the probability on a full n-gram starting
			# with this word, but only over those n-grams that have been seen (counts > 0);
			# Calculate the sum over all words of the probability on a full (n-b)-gram starting
			# with this word, but only over those (n-b)-grams for which the n-grams that have 
			# been seen (counts > 0). Then calculate Alpha = (1-SumNgrams)/(1-Sum(N-b)grams)
			# If the (n-b)-gram doesn't exist, switch to a higher b.
			if(!$CurrentCount)
			{
				# Calculate probability mass to redistribute
				my $Pnonzero = 0;
				$Pnonzero = $CurrentNonZeroDenominatorCount / $CurrentDenominatorCount if $CurrentNonZeroDenominatorCount;

				my $BackoffDepth = $n-1;
				my $N1context = join(" ", @WordSequence[($j-($n-1))..($j-1)]);
				while(!$CurrentCount && $BackoffDepth)
				{
					# Sum smoothed counts of all and all non-zero Ngrams ending in N gram
					my $NonZeroN1Count = 0;
					my $TotalN1Count = 0;
				    my $Nbcontext = join(" ", @WordSequence[($j-($BackoffDepth-1))..($j-1)]);
					unless($AlphaContext[$BackoffDepth]->{$N1context} > 0)
					{
						foreach $CurrentUnigram (keys(%{$NgramTables[1]}))
						{
							# Count all N-1grams for which the Ngram exists
							my $TestNgramWords = "$N1context $CurrentUnigram";

							# Count the  N-1grams
							my $CurrentN1gramWords = $BackoffDepth > 1 ?
											"$Nbcontext $CurrentUnigram"
											: $CurrentUnigram;
							my $CurrentN1gramCount = $NgramTables[$BackoffDepth]->{$CurrentN1gramWords};
							# The normalizing count of ALL N-1grams
							$TotalN1Count += $CurrentN1gramCount || $SmoothedNgramCounts[$BackoffDepth][0];

							# Skip if Ngram doesn't exist
							next unless $NgramTables[$n]->{$TestNgramWords};
							# These count!!!
							$NonZeroN1Count += $CurrentN1gramCount;
						};
						$Alpha = $TotalN1Count ? (1 - $Pnonzero) / (1 - $NonZeroN1Count / $TotalN1Count) : 0;
						$AlphaContext[$BackoffDepth]->{$N1context} = $Alpha;
					}
					else
					{
						$Alpha = $AlphaContext[$BackoffDepth]->{$N1context};
					};
					$CurrentNgramWords = join(" ", @WordSequence[($j-($BackoffDepth-1))..$j]);
					$CurrentCount = $NgramTables[$BackoffDepth]->{$CurrentNgramWords};
					$CurrentDenominatorCount = $TotalN1Count;
					# Step down a level
					--$BackoffDepth;
				};
			};
			
			# If backoff did NOT help, use the original smoothed count
			$CurrentCount = $smoothedCount if !$CurrentCount && $backoff;
			
			# Print out detailed calculations
			if($verbose)
			{
				printf("Prob: %5.3g * %7.3g / %7.3g\t", $Alpha, $CurrentCount, $CurrentDenominatorCount) if $Alpha != 1;
				printf("Prob:          %7.3g / %7.3g\t", $CurrentCount, $CurrentDenominatorCount) unless $Alpha != 1;
				print "\"$CurrentNgramWords\" ";
			};
			
		};
		
		# Check for zero (unsmoothened or OOV) counts
		if($CurrentCount <= 0)
		{
			print STDERR "Word ", $j - ($n-1), ": \"$WordSequence[$j]\" ",
						  "resulted in an invalid count of $CurrentCount, ", 
						  "N-gram discarded: '$CurrentNgramWords'\n";
			next;
		};

		# New probability
		my $PartialLogProbability = -(log($Alpha) + log($CurrentCount) - log($CurrentDenominatorCount)) / log(2.0);
		print "\tLogProb: $LogProbability += $PartialLogProbability\n"
				if $verbose;
		$LogProbability += $PartialLogProbability;
		
		# Store perplexity
		$SumLogProb += $PartialLogProbability;
		++$Nwords;
	};
	# Print output current line
	if(!$Perplexity || $verbose)
	{
		print "$LogProbability (", $LogProbability / (scalar(@WordSequence) - ($n-1)),"/word) = -log2(Prob(\"@WordSequence[$n-1..$#WordSequence]\"))\n";
	};
};

# Print total perplexity
print "Log Perplexity over $Nwords words = ", $SumLogProb / $Nwords, " (", 2**($SumLogProb / $Nwords), ")\n" if $Perplexity;
