#! /usr/local/bin/perl
#
if(grep(/^\-(h|\-help)$/i, @ARGV))
{
print STDOUT << "ENDOFHELP";
NCD_Ngramcount.pl Ngramcountfile1 Ngramcountfile2

Use N-gram counts to estimate a 
Normalized Compression Distance
and write to STDOUT.
Ngramcountfiles are two plain ASCII count files.

--license
-l
Print license information

--help
-h
This message

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: R.J.J.H.vanSon\@uva.nl
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;
};

##############################################################################

my %CountFile1 = ();
my %CountFile2 = ();
my %CombinedCountFile = ();

my ($Counts1, $Counts2, $CombinedCounts) = (0, 0, 0);

my $File1 = shift(@ARGV);
my $File2 = shift(@ARGV);

open(FILE1, "<$File1") || die "$File1: $!\n";
while(<FILE1>)
{
    # Skip comments
    next if /^\s*\# /;
    # Skip empty lines
    next unless /\S/;
    
    if(/^\s*(\d+)\t/)
    {
        my $Count = $1;
        my $Value = $';
        $Value =~ s/\s+//ig;
        
        $CountFile1{$Value} += $Count;
        $Counts1 += $Count;
        
        $CombinedCountFile{$Value} += $Count;
        $CombinedCounts += $Count;
        
    };
    
    
};
close(FILE1);

open(FILE2, "<$File2") || die "$File2: $!\n";
while(<FILE2>)
{
    # Skip comments
    next if /^\s*\# /;
    # Skip empty lines
    next unless /\S/;
    
    if(/^\s*(\d+)\t/)
    {
        my $Count = $1;
        my $Value = $';
        $Value =~ s/\s+//ig;
        
        $CountFile2{$Value} += $Count;
        $Counts2 += $Count;
        
        $CombinedCountFile{$Value} += $Count;
        $CombinedCounts += $Count;
        
    };
    
    
};
close(FILE2);

#######################################################################

# Calculate the Kullback-Leibler distances
# KL(1,2)
my $KL12 = 0;
foreach $Word (keys(%CountFile1))
{
    my ($p, $q) = (0, 0);
    if(exists($CountFile2{$Word}))
    {
        $p = $CountFile1{$Word} / $Counts1;
        $q = $CountFile2{$Word} / $Counts2;
    };
    
    if($q > 0 && $p > 0)
    {
        $KL12 += $p * log($p / $q) / log(2);
    };
};

my $KL21 = 0;
foreach $Word (keys(%CountFile2))
{
    my ($p, $q) = (0, 0);
    if(exists($CountFile1{$Word}))
    {
        $p = $CountFile2{$Word} / $Counts2;
        $q = $CountFile1{$Word} / $Counts1;
    };
    
    if($q > 0 && $p > 0)
    {
        $KL21 += $p * log($p / $q) / log(2);
    };
};

# Calculate type counts
my $TypeCount1 = scalar(keys(%CountFile1));
my $TypeCount2 = scalar(keys(%CountFile2));

my $MinFileCount = $TypeCount1 <= $TypeCount2 ? $TypeCount1 : $TypeCount2;
my $MaxFileCount = $TypeCount1 >= $TypeCount2 ? $TypeCount1 : $TypeCount2;
$MaxFileCount = 1 unless $MaxFileCount > 0;

my $CombinedTypeCount = scalar(keys(%CombinedCountFile));

my $NCD = ($CombinedTypeCount - $MinFileCount) / $MaxFileCount;

# Print results
print "# File     \tTypes  \tCounts  \n";
print "# $File1:\t$TypeCount1\t$Counts1\n";
print "# $File2:\t$TypeCount2\t$Counts2\n";
print "# Combined:\t$CombinedTypeCount\t$CombinedCounts\n";
print "# NCD = $NCD\n";
print "# \n";
print "# KL($File1, $File2) = $KL12\n";
print "# KL($File2, $File1) = $KL21\n";
print "# \n";




