#!/usr/local/bin/perl

use Math::MatrixReal;

require glob("~ihh/perl/standard.pl");

$identifier = "MODELFILE_1.0";

$usage .= "@{[nameonly($0)]} -- calculate entropy of models in ModelFile 1.0 format\n";
$usage .= "\n";
$usage .= "Usage: @{[nameonly($0)]} [files...]\n";
$usage .= "\n";

while (@ARGV) {
    last unless ($ARGV[0] =~ /^-/);
    $opt = shift;
    die "Unknown option $opt\n$usage";
}

if (@ARGV==0) { @ARGV = ("-") }

foreach $file (@ARGV) {
    open file or die "Couldn't open $file: $!";
    while (<file>) {
	next unless /\S/;   # skip whitespace
	
	$mfile = ($file eq "-") ? "<STDIN>" : "File $file";
	$line = $.;

	($id,$alphabet,$states,$links) = split;
	($id eq $identifier) or die "$mfile line $line: expected $identifier\n";
	
	$alen = length($alphabet);

	warn "Initializing matrices...\n";

	$A = Math::MatrixReal->new($states,$states);
	$b = Math::MatrixReal->new($states,1);

	$A->one;
	$A->negate($A);
	$b->zero;
	
	warn "Reading file...\n";

	%p = ();
	for ($i=0;$i<$links;$i++) {
	    $_ = <file>;
	    ($from,$to,$dx,$dy,@t) = split;
	    $num = $alen ** ($dx+$dy);
	    (@t == ($alen**($dx+$dy))) or warn "$mfile line $line: expected $num transition probabilities\n";
	    foreach $prob (@t) {
		if ($prob>0) {
		    $A->assign($from+1, $to+1, $A->element($from+1,$to+1) + $prob);
		    $b->assign($from+1, 1, $b->element($from+1,1) - $prob*log($prob));
		}
	    }
	}
    }

    warn "Normalizing matrices...\n";
    ($norm_A,$norm_b) = $A->normalize($b);
    warn "LR-decomposing...\n";
    $LR = $norm_A->decompose_LR();
    warn "LR-solving...\n";
    ($dim,$x,$B) = $LR->solve_LR($norm_b);

    $entropy = - $x->element(1,1);
    print "File $file\tEntropy $entropy\n";
}


