#!/usr/local/bin/perl

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


$usage .= "@{[nameonly($0)]} -- convert HMMER2.0 format Plan7 model files into ModelFile 1.0 format\n";
$usage .= "\n";
$usage .= "Usage: @{[nameonly($0)]} [options] filename\n";
$usage .= "\n";
$usage .= "Options:\n";
$usage .= " -debug   Print debugging information (state labels)\n";
$usage .= " -ratio   Do not factor out null model (yields an un-normalised model!)\n";
$usage .= " -corebug Replicate bug in core_algorithms.c for transition from begin to final match state\n";$
usage .= "\n";

while (@ARGV) {
    last unless ($ARGV[0] =~ /^-/);
    $opt = shift;
    if ($opt eq "-debug") { $debug = 1 }
    elsif ($opt eq "-ratio") { $ratio = 1 }
    elsif ($opt eq "-corebug") { $corebug = 1 }
    else { die "Unknown option: $opt\n\n$usage" }
}

(@ARGV==1) or die $usage;

# first read in the model

$_ = <>;
unless ($_ eq "HMMER2.0\n") { die "Format error: HMMER2.0 identifier expected"; }

while (<>) {
    ($tag,$val) = split /\s+/,$_,2;
    chop $val;
    $val{$tag} = $val;
    last if /^HMM/;
}

$len = $val{LENG};

$alphabet = $val{HMM};
$alphabet =~ s/\s//g;

@nult = split /\s+/,$val{NULT};
foreach (@nult) { $_ = ascii2prob($_) }
($p1,$p2) = @nult;

@nule = split /\s+/,$val{NULE};
foreach (@nule) { $_ = ascii2prob($_) / length($alphabet) }
unless (@nule==length($alphabet)) { die "Expected ".length($alphabet)." null emission probabilities; found ".scalar(@nule); }

@xt = split /\s+/,$val{XT};
foreach (@xt) { $_ = ascii2prob($_) }
($nb,$nn,$ec,$ej,$ct,$cc,$jb,$jj) = @xt;

$dummy = <>;

$_ = <>;
@p = split;
foreach (@p) { $_ = ascii2prob($_) }
($bm1,$bi1,$bd1) = @p;
($bi1==0) or die "Expected b->i1 transition to have probability zero";

for ($s=1;$s<=$len;$s++) {
    $_ = <>;
    ($n,@p) = split;
    for ($i=0;$i<@p;$i++) { $p[$i] = ascii2prob($p[$i]) }
    unless ($ratio) { for ($i=0;$i<@p;$i++) { $p[$i] *= $nule[$i] } }
    @{$matchEmit[$s]} = @p;
    
    unless ($n==$s) { die "Expected state $s, found state $n"; }
    unless (@p==length($alphabet)) { die "Expected ".length($alphabet)." match emission probabilities; found ".scalar(@p); }
    
    $_ = <>;
    ($dummy,@p) = split;
    for ($i=0;$i<@p;$i++) { $p[$i] = ascii2prob($p[$i]) }
    unless ($ratio) { for ($i=0;$i<@p;$i++) { $p[$i] *= $nule[$i] } }
    @{$insertEmit[$s]} = @p;

    unless (@p==length($alphabet)) { die "Expected ".length($alphabet)." insert emission probabilities; found ".scalar(@p); }

    $_ = <>;
    ($dummy,@p) = split;

    unless (@p==9) { die "Expected 9 insert emission probabilities; found ".scalar(@p); }

    foreach (@p) { $_ = ascii2prob($_) }
    ($mm[$s],$mi[$s],$md[$s],$im[$s],$ii[$s],$dm[$s],$dd[$s],$bm[$s],$me[$s]) = @p;

}

if ($ratio) {
    $nn /= $p1;
    $cc /= $p1;
    $jj /= $p1;
    $ct /= $p2;
    for ($s=1;$s<=$len;$s++) {
	$bm[$s] /= $p1;
	$mm[$s] /= $p1;
	$im[$s] /= $p1;
	$dm[$s] /= $p1;
	$mi[$s] /= $p1;
	$ii[$s] /= $p1;
    }
    foreach (@nule) { $_ = 1 }
}

# now add the links

for ($s=1;$s<=$len;$s++) {

# do the match->match, insert->match, delete->match links first
# to avoid pointless null cycles

    @bme = @mme = @ime = @dme = @{$matchEmit[$s]};
    @mie = @iie = @{$insertEmit[$s]};

    foreach (@dme) { $_ *= $dm[$s-1] }
    addLink(D($s-1),M($s),1,@dme) if ($s>2);

    foreach (@mme) { $_ *= $mm[$s-1] }
    addLink(M($s-1),M($s),1,@mme) if ($s>1);

    foreach (@ime) { $_ *= $im[$s-1] }
    addLink(I($s-1),M($s),1,@ime) if ($s>1);

# next do the begin->match links

    unless ($s==$len && $corebug) {
	$bm = $bm[$s];
	if ($s>1) {
	    $bd = $bd1;
	    for ($i=1;$i<$s-1;$i++) { $bd *= $dd[$i] }
	    $bm += $bd * $dm[$s-1];
	}
    }
    foreach (@bme) { $_ *= $bm }
    addLink(B(),M($s),1,@bme);
    foreach (@mie) { $_ *= $mi[$s] }
    addLink(M($s),I($s),1,@mie) if ($s<$len);

# do all the other internal links

    foreach (@iie) { $_ *= $ii[$s] }
    addLink(I($s),I($s),1,@iie) if ($s<$len);
    
    addLink(M($s),D($s+1),0,$md[$s]) if ($s<$len-1);

    addLink(D($s),D($s+1),0,$dd[$s]) if ($s>1 && $s<$len-1);

    $me = $me[$s];
    if ($s<$len) {
	$de = $md[$s];
	for ($i=$s+1;$i<$len;$i++) { $de *= $dd[$i] }
	$me += $de;
    }
    addLink(M($s),E(),0,$me);
}

for ($s=1;$s<=$len;$s++) {
}

# finally, the special links

@nne = @cce = @jje = @nule;
foreach (@nne) { $_ *= $nn }
foreach (@cce) { $_ *= $cc }
foreach (@jje) { $_ *= $jj }

addLink(N(),N(),1,@nne);
addLink(C(),C(),1,@cce);
addLink(J(),J(),1,@jje);

addLink(S(),N(),0,1);
addLink(N(),B(),0,$nb);
addLink(E(),C(),0,$ec);
addLink(E(),J(),0,$ej);
addLink(J(),B(),0,$jb);
addLink(C(),T(),0,$ct);

# the following link ensures a probabilistically valid model but represents a null cycle & so is omitted:
#
# $be = $bd1;
# for ($i=1;$i<=$len-1;$i++) { $be *= $dd[$i] }
# addLink(B(),E(),0,$be);

print "MODELFILE_1.0 $alphabet @{[states()]} @{[scalar(@links)]}\n";
print join "\n",@links;

if ($debug) {
    warn "N=@{[N()]} B=@{[B()]} M1=@{[M(1)]} I1=@{[I(1)]} D2=@{[D(2)]} E=@{[E()]} J=@{[J()]} C=@{[C()]} len=$len\n";
}

sub addLink {
    my ($from,$to,$emitting,@emit) = @_;
    my $link = "$from $to ";
    $link .= $emitting ? "1 " : "0 ";
    $link .= "0 ";
    $link .= join " ",@emit;
    push @links,$link;
}

sub ascii2prob {
    my ($a) = @_;
    if ($a =~ /\*/) { return 0 }
    return exp(log(2)*$a/1000);
}

# for an interpretation of the following state labels, see the Plan 7 architecture

sub S { 0 }
sub T { 1 }
sub N { 2 }
sub M { N()+1 + ($_[0]-1) }
sub I { M($len)+1 + ($_[0]-1) }
sub E { I($len-1)+1 }
sub C { E()+1 }
sub J { C()+1 }
sub B { J()+1 }
sub D { B()+1 + ($_[0]-2) }
sub states { D($len-1)+1 }

