#!/usr/bin/perl
use strict;
use warnings;

#
# Automata are represented as hash tables containing three sets :
# (see the example of the SNP/Indel below)
#
#  ( SIZE => [integer that gives the number of states]
#    INITPROB => [hash table that gives the probability to be at a given state (key) for any state]
#    TRANSITION => [ hash (key "from") of hash (key "to") of hash (key "char") that gives the probability
#    of reaching the state "to" and  generating the char "char" starting from state "from"]
# )

# SNP/Indel model : probability of SNP and IndelRate
my $pr_snp_or_indel  = 0.01;
my $indel_rate       = 0.15;

# SNP/Indel automaton : three states that generate Matches (0),  SNP (1) and Indel (2)
my %SNP_AA =
    (
      SIZE => 3, # nbstates
      INITPROB => {  # initial dist
          0 => 1.0-$pr_snp_or_indel,
          1 => $pr_snp_or_indel*(1.0-$indel_rate),
          2 => $pr_snp_or_indel*$indel_rate,
      },
      TRANSITIONS => {
          0 => # from
          {
              0 => # to
              {
                  "M" => 1.0-$pr_snp_or_indel,
                  #"E" => 0.0,
                  #"I" => 0.0,
              },
              1 => # to
              {
                  #"M" => 0.0,
                  "E" => $pr_snp_or_indel*(1.0 - $indel_rate),
                  #"I" => 0.0,
              },
              2 => # to
              {
                  #"M" => 0.0,
                  "E" => $pr_snp_or_indel*$indel_rate,
                  #"I" => 0.0,
              },
          },
          1 => # from
          {
              0 => # to
              {
                  #"M" => 0.0,
                  "E" => (1.0-$pr_snp_or_indel),
                  #"I" => 0.0,
              },
              1 => # to
              {
                  "M" => 0.25*$pr_snp_or_indel*(1.0-$indel_rate),
                  "E" => 0.75*$pr_snp_or_indel*(1.0-$indel_rate),
                  #"I" => 0.0,
              },
              2 => # to
              {
                  "M" => 0.25*$pr_snp_or_indel*$indel_rate,
                  "E" => 0.75*$pr_snp_or_indel*$indel_rate,
                  #"I" => 0.0,
              },
          },
          2 => # from
          {
              0 => # to
              {
                  #"M" => 0.0,
                  #"E" => 0.0,
                  "I" => (1.0-$pr_snp_or_indel),
              },
              1 => # to
              {
                  #"M" => 0.0,
                  #"E" => 0.0,
                  "I" => $pr_snp_or_indel*(1.0-$indel_rate),
              },
              2 => # to
              {
                  #"M" => 0.0,
                  #"E" => 0.0,
                  "I" => $pr_snp_or_indel*$indel_rate,
              },
          },
      },
    );

# Output the automaton in simple form for debugging
sub print_aa($) {
    my %Automaton = @_;
    my $size  = $Automaton{SIZE};
    print "size : ".$size."\n";
    print "initprob : (\n";
    while( my ($s,$p) = each (%{$Automaton{INITPROB}})) {
        print "\tstate : $s\tpr : $p\n";
    }
    print ")\n";
    print "(\n";
    while( my ($s,$qlp) = each (%{$Automaton{TRANSITIONS}})) {
        print "\t".$s.",\n";
        while( my ($q,$lp) = each (%{$qlp})) {
            print "\t\t".$q.",\n";
            while( my ($l,$p) = each (%{$lp})) {
                print "\t\t\t".$l." : ".$p.",\n";
            }
        }
    }
    print ")\n";
}

# Output the automaton for the iederalinear "C" program (generate "C" tables)
sub print_aa_c($) {
    my %Automaton = @_;
    my $size  = $Automaton{SIZE};
    print "#define NB_STATES ".$size."\n";
    print "#define ALIGN_ALPHABET_SIZE 3\n";
    print "double automaton_init_prob[NB_STATES] = {";
    my $comma = 0;
    my $i = 0;
    for(my $s = 0; $s < $size ; $s++) {
        if ($comma) {
            print ",";
        }
        if (defined $Automaton{INITPROB}{$s}) {
            print "/*".$s."*/ ".$Automaton{INITPROB}{$s};
        } else {
            print "/*".$s."*/ 0.0";
        }
        $comma = 1;
    }
    print "};\n";

    print "int automaton_size[NB_STATES] =  {";
    $comma = 0;
    for(my $s = 0; $s < $size; $s++) {
        my $size_t =  keys(%{$Automaton{TRANSITIONS}{$s}});
        if ($comma) {
            print ",";
        }
        print $size_t;
        $comma = 1;
    }
    print "};\n";

    print "int automaton_states[NB_STATES][NB_STATES] = {\n";
    $comma = 0;
    for(my $s = 0; $s < $size ; $s++) {
        my $comma_s = 0;
        if ($comma) {
            print ",\n\t{";
        } else {
            $comma = 1;
            print "\t{";
        }

        foreach my $q (sort {$a <=> $b} (keys(%{$Automaton{TRANSITIONS}{$s}}))) {
            if ($comma_s) {
                print ",";
            }
            $comma_s = 1;
            print $q;
        }
        print "}";
    }
    print "};\n";

    print "double automaton_prob[NB_STATES][ALIGN_ALPHABET_SIZE][NB_STATES] = {\n";
    $comma = 0;
    for (my $s=0; $s < $size ;$s++) {
        my $comma_q = 0;
        if ($comma) {
            print "\t,\n";
        }
        print "\t{";
        $comma = 1;
        print " /* ".$s." */\n";
        foreach my $l ("I","E","M") {
            my $comma_s = 0;
            if ($comma_q) {
                print ",";
            }
            print "\n\t\t{";
            $comma_q = 1;
            print " /* ".$l." */ ";
            foreach my $q (sort {$a <=> $b} (keys(%{$Automaton{TRANSITIONS}{$s}}))) {
                if ($comma_s) {
                    print ",";
                }
                $comma_s = 1;
                if (defined $Automaton{TRANSITIONS}{$s}{$q}{$l}) {
                    print " /* ".$q." */ ".$Automaton{TRANSITIONS}{$s}{$q}{$l}."";
                } else {
                    print " /* ".$q." */ 0.0";
                }
            }
            print "}";
        }
        print "\n\t}\n";
    }
    print "};\n";
 }



# Output the automaton for the iedera program [lossless case]
sub print_aa_iedera_raw($) {
    my @lookupletter = ("I","E","M");
    my %Automaton = @_;
    # give the size
    print "".($Automaton{SIZE})."\n";
    # set the states
    for (my $s=0; $s < $Automaton{SIZE} ; $s++) {
        if (defined $Automaton{FINAL}{$s} &&  $Automaton{FINAL}{$s} > 0) {
            print "".($s)."\t1\n"; # state number + final
        } else {
            print "".($s)."\t0\n"; # state number + final
        }
        for (my $c = 0 ; $c < 3 ; $c++) {
            my %set = ();
            foreach my $q (sort {$a <=> $b} (keys(%{$Automaton{TRANSITIONS}{$s}}))) {
                if (defined $Automaton{TRANSITIONS}{$s}{$q}{$lookupletter[$c]} && $Automaton{TRANSITIONS}{$s}{$q}{$lookupletter[$c]} > 0) {
                    my $pr = $Automaton{TRANSITIONS}{$s}{$q}{$lookupletter[$c]};
                    $set{$q} = $pr;
                }
            }
            print "\t".$c."\t".(scalar (keys %set))."\n"; # a + alistsize
            for my $q (sort {$a <=> $b} (keys(%set))) {
                print "\t\t".($q)."\t".($set{$q})."\n"; # q + pr
            }
        }
    }
 }


# Output the automaton for the iedera program [lossy case]
sub print_aa_iedera_prob($) {
    my %Automaton = @_;
    my $size  = $Automaton{SIZE};
    print "".($size+2)."\n"; # + final state 0 and init state 1

    # final state 0
    print "0\t1\n"; # nb state + final
    for (my $c = 0 ; $c < 3 ; $c++) {
        print "\t".$c."\t1\n"; # c + alistsize
        print "\t\t0\t0.333333\n"; # q + pr
    }

    # init state 1
    my $nb_init_probs = 0;
    for (my $s=0; $s < $size ; $s++) {
        if (defined $Automaton{INITPROB}{$s}) {
            my $nb_states_transitions_from_s = keys(%{$Automaton{TRANSITIONS}{$s}});
            $nb_init_probs += $nb_states_transitions_from_s;
        }
    }

    # set init state 1
    print "1\t0\n"; # nb state + final
    my @lookupletter = ("I","E","M");
    for (my $c = 0 ; $c < 3 ; $c++) {
        my %set = ();
        for (my $s=0; $s < $size ; $s++) {
            if (defined $Automaton{INITPROB}{$s} && $Automaton{INITPROB}{$s} > 0) {
                foreach my $q (sort {$a <=> $b} (keys(%{$Automaton{TRANSITIONS}{$s}}))) {
                    if (defined $Automaton{TRANSITIONS}{$s}{$q}{($lookupletter[$c])} && $Automaton{TRANSITIONS}{$s}{$q}{($lookupletter[$c])} > 0) {
                        my $pr =
                            ($Automaton{TRANSITIONS}{$s}{$q}{($lookupletter[$c])})
                            *
                            ($Automaton{INITPROB}{$s});
                        if (defined ($set{$q})) {
                            $set{$q} += $pr;
                        } else {
                            $set{$q} = $pr;
                        }
                    }
                }
            }
        }
        print "\t".$c."\t".(scalar (keys %set))."\n"; # a + alistsize
        for my $q (sort {$a <=> $b} (keys(%set))) {
            print "\t\t".($q + 2)."\t".($set{$q})."\n"; # q + pr
        }
    }

    # set the other states
    for (my $s=0; $s < $size ; $s++) {
        print "".($s + 2)."\t0\n"; # nb state + final
        for (my $c = 0 ; $c < 3 ; $c++) {
            my %set = ();
            foreach my $q (sort {$a <=> $b} (keys(%{$Automaton{TRANSITIONS}{$s}}))) {
                if (defined $Automaton{TRANSITIONS}{$s}{$q}{$lookupletter[$c]} && $Automaton{TRANSITIONS}{$s}{$q}{$lookupletter[$c]} > 0) {
                    my $pr = $Automaton{TRANSITIONS}{$s}{$q}{$lookupletter[$c]};
                    $set{$q} = $pr;
                }
            }
            print "\t".$c."\t".(scalar (keys %set))."\n"; # a + alistsize
            for my $q (sort {$a <=> $b} (keys(%set))) {
                print "\t\t".($q + 2)."\t".($set{$q})."\n"; # q + pr
            }
        }
    }
}



# Build a simple (non correlated) model
#   $len      : length of the alignment
#   $pr_start : probability of read error at the start of the alignment (in the non correlated case)
#   $pr_end   : probability of read error at the end of the alignment (in the non correlated case)
#   $indel_rate : should be 0 unless they are indel during the reading process (??)

sub build_error_simple_aa($$$$) {
    my ($len,$pr_start,$pr_end,$indel_rate) = @_;
    # size
    my %A = (SIZE => ($len+1));
    # init prob
    $A{INITPROB}{0} = 1.0;
    for(my $i=1;$i<=$len;$i++) {
        $A{INITPROB}{$i} = 0.0;
    }
    # transitions
    for(my $l=0;$l<=$len;$l++) {
        my $pr_nM = ($pr_start*($len-$l) + $pr_end*($l)) / ($len);
        my $pr_M  = 1.0 - $pr_nM;
        my $new   = $l+1;
        if ($new > $len) {
            $new  = $len;
        }
        $A{TRANSITIONS}{$l}{$new}{"M"} = $pr_M;
        $A{TRANSITIONS}{$l}{$new}{"E"} = $pr_nM*(1.0-$indel_rate);
        $A{TRANSITIONS}{$l}{$new}{"I"} = $pr_nM*$indel_rate;
    }
    return %A;
}


# Build the complex model used for SOLiD read errors
#   $len      : length of the alignment
#   $pr_start : probability of read error at the start of the alignment (in the non correlated case)
#   $pr_end   : probability of read error at the end of the alignment (in the non correlated case)
#   $indel_rate : should be 0 unless they are indel during the reading process (??)
#   $e_freq : gives the frequency to jump from a "non correlated" black arrow to a "correlated" red arrow (should be low).
#   $e_err : probability of read error on "corrrelated" red arrow (should be higher than $pr_end)


sub build_error_complex_aa($$$$$) {
    my ($len,$pr_start,$pr_end,$indel_rate,$e_freq,$e_err) = @_;
    my $nbstates = (2**5)*($len+1);
    # size
    my %A = (SIZE => ($nbstates));
    # init prob
    my $currentstate = 0;

    for(my $pattern = 0 ; $pattern < 2**5 ; $pattern++) {
        my $nberr = 0;
        for(my $i=0;$i<5;$i++) {
            if ($pattern & (1<<$i)) {
                $nberr++;
            }
        }
        my $pr = ($e_freq**$nberr) * (1-$e_freq)**(5-$nberr);
        $A{INITPROB}{$currentstate} = $pr;
        $currentstate++;
    }

    for(my $i = 0 ; $i < (2**5)*$len ; $i++) {
        $A{INITPROB}{$currentstate} = 0.0;
        $currentstate++;
    }

    # transitions
    $currentstate = 0;
    for(my $l = 0; $l <= $len; $l++ ) {
        for(my $pattern = 0 ; $pattern < 2**5 ; $pattern++) {
            my $newstate = $currentstate + 2**5;
            if ($newstate >= $nbstates) {
                $newstate = $currentstate;
            }

            # necklaces hard coded for "5"
            if ($pattern & (1 << ($l % 5))) {
                # red arrow : fixed "$e_err" error rate (usually high eg 0.25)
                my $pr_nM = $e_err;
                my $pr_M  = 1.0 - $pr_nM;
                $A{TRANSITIONS}{$currentstate}{$newstate}{"M"} = $pr_M;
                $A{TRANSITIONS}{$currentstate}{$newstate}{"E"} = $pr_nM*(1.0-$indel_rate);
                $A{TRANSITIONS}{$currentstate}{$newstate}{"I"} = $pr_nM*$indel_rate;
            } else {
                # black arrow : normal (and increasing to the end) error rate
                # two possibilities : the arrow stays black (1-$e_freq) or goes red ($e_freq)
                my $newstate_red = $currentstate + 2**5 + (1 << ($l%5));
                if ($newstate_red >= $nbstates) {
                    $newstate_red = $currentstate;
                }

                my $pr_nM_red = $e_err;
                my $pr_M_red  = 1.0 - $pr_nM_red;

                $A{TRANSITIONS}{$currentstate}{$newstate_red}{"M"} = $pr_M_red*$e_freq;
                $A{TRANSITIONS}{$currentstate}{$newstate_red}{"E"} = $pr_nM_red*(1.0-$indel_rate)*$e_freq;
                $A{TRANSITIONS}{$currentstate}{$newstate_red}{"I"} = $pr_nM_red*$indel_rate*$e_freq;

                my $pr_nM = ($pr_start*($len-$l) + $pr_end*($l)) / ($len);
                my $pr_M  = 1.0 - $pr_nM;

                if (defined $A{TRANSITIONS}{$currentstate}{$newstate}) {
                    $A{TRANSITIONS}{$currentstate}{$newstate}{"M"} += $pr_M*(1.0-$e_freq);
                    $A{TRANSITIONS}{$currentstate}{$newstate}{"E"} += $pr_nM*(1.0-$indel_rate)*(1-$e_freq);
                    $A{TRANSITIONS}{$currentstate}{$newstate}{"I"} += $pr_nM*$indel_rate*(1-$e_freq);
                } else {
                    $A{TRANSITIONS}{$currentstate}{$newstate}{"M"} = $pr_M*(1.0-$e_freq);
                    $A{TRANSITIONS}{$currentstate}{$newstate}{"E"} = $pr_nM*(1.0-$indel_rate)*(1-$e_freq);
                    $A{TRANSITIONS}{$currentstate}{$newstate}{"I"} = $pr_nM*$indel_rate*(1-$e_freq);
                }
            }
            $currentstate++;
        }
    }
    return %A;
}


# Dominant symbol ("I"ndel > "E"rror > "M"atch)
sub dominant($$) {
    my ($a,$b) = @_;
    if (($a eq "I") || ($b eq "I")) {
        return "I";
    } elsif (($a eq "E") || ($b eq "E")) {
        return "E";
    } else {
        return "M";
    }
}


# Does the product of two models $A_ref and $B_ref : generates new automaton
sub product_aa($$) {
    my ($A_ref,$B_ref) = @_;
    my %A = %$A_ref;
    my %B = %$B_ref;

    my $A_size = $A{SIZE};
    my $B_size = $B{SIZE};

    # size
    my $C_size = $A_size * $B_size;
    my %C      = (SIZE => $C_size);

    my @ABstates_stack = ();

    # init prob (on the first reachable states ...)
    my %instack = ();
    for (my $currentstate_a = 0; $currentstate_a < $A_size; $currentstate_a++) {
        if ($A{INITPROB}{$currentstate_a} > 0) {
            for (my $currentstate_b = 0; $currentstate_b < $B_size; $currentstate_b++) {
                my $currentstate   = $currentstate_a*$B_size+$currentstate_b;
                if ($B{INITPROB}{$currentstate_b} > 0) {
                    $C{INITPROB}{$currentstate} =
                        $A{INITPROB}{$currentstate_a} * $B{INITPROB}{$currentstate_b};
                    push @ABstates_stack,$currentstate_a;
                    push @ABstates_stack,$currentstate_b;
                    $instack{$currentstate} = 1;
                }
            }
        }
    }

    # transitions
    while ($#ABstates_stack > 0) {
        my $currentstate_b = pop @ABstates_stack;
        my $currentstate_a = pop @ABstates_stack;
        my $currentstate   = $currentstate_a*$B_size+$currentstate_b;
        while( my ($nextstate_a,$lp_a) = each (%{$A{TRANSITIONS}{$currentstate_a}})) {
            while( my ($nextstate_b,$lp_b) = each (%{$B{TRANSITIONS}{$currentstate_b}})) {
                my $nextstate = $nextstate_a*$B_size+$nextstate_b;
                # pair of transitions on letters
                while( my ($l_a,$p_a) = each (%{$lp_a})) {
                    while( my ($l_b,$p_b) = each (%{$lp_b})) {
                        if ($p_a > 0 && $p_b > 0) {
                            my $l = &dominant($l_a,$l_b);
                            if (defined $C{TRANSITIONS}{$currentstate}{$nextstate}{$l}) {
                                $C{TRANSITIONS}{$currentstate}{$nextstate}{$l} += $p_a*$p_b;
                            } else {
                                $C{TRANSITIONS}{$currentstate}{$nextstate}{$l} = $p_a*$p_b;
                            }
                        }
                    }
                }
                if (!defined $instack{$nextstate}) {
                    push @ABstates_stack,$nextstate_a;
                    push @ABstates_stack,$nextstate_b;
                    $instack{$nextstate} = 1;
                }
            }
        }
    }
    return %C;
}



# This funtion generate an automaton with at most total_cost_threshold
# provided that Indels, Errors and SNP have a given cost
#


sub lossless_aa($$$$) {
    my ($snp_cost,$err_cost,$indel_cost,$total_cost_threshold) = @_;


    # 1) check if the given costs are ok
    if ($snp_cost < 0 || $err_cost < 0 || $indel_cost < 0 || $total_cost_threshold < 0) {
        print STDERR "all costs must be >= 0";
        exit(0);
    }

    if ($snp_cost > 2*$err_cost) {
        print STDERR "invalid SNP cost : must be \"<= 2 x Error cost\"";
        exit(0);
    }

    # 2) build the automaton
    my %A = (
        SIZE => 2,
        INITPROB => {  # initial dist
            1 => 1.0
        },
        TRANSITIONS => {
         # "0" = Fail state
            0 => # from
            {
                0 => # to
                {
                    "M" => 0.333333,
                    "E" => 0.333333,
                    "I" => 0.333333,
                }
            },
        },

        FINAL => {
            # "0" = Fail state
            0 => 1
        },
     );

    # stack to set states that have to be processed
    my @ScoreStates_stack = ();
    push @ScoreStates_stack,1; #state
    push @ScoreStates_stack,0; #score

    # stack to store states on given score;
    my %MatchStateOfScore = ( 0 => 1 );
    my %IndelStateOfScore = ();

    while ($#ScoreStates_stack > 0) {

        my $score = pop @ScoreStates_stack;
        my $state = pop @ScoreStates_stack;

        # A) score is >= 0 -> "Match" (otherwise "indel" encoded)
        if ($score >= 0) {

            # add 1st Error if needed (check before that SNP,Error and Indel do no give a fail state)

            # selfloop on a match
            $A{TRANSITIONS}{$state}{$state}{"M"} = 0.33333;

            # error and indel only (no SNP)
            if ($score + $snp_cost > $total_cost_threshold) {

                # error? (no SNP)
                if ($score + $err_cost > $total_cost_threshold) {
                    # no error, (no SNP)
                    $A{TRANSITIONS}{$state}{0}{"E"} = 0.33333;
                } else {
                    # error, (no SNP)
                    my $next_state;
                    if (defined $MatchStateOfScore{($score + $err_cost)}) {
                        $next_state = $MatchStateOfScore{($score + $err_cost)};
                    } else {
                        $next_state = $A{SIZE};
                        $MatchStateOfScore{($score + $err_cost)} = $next_state;
                        $A{SIZE} = $next_state + 1;
                        push @ScoreStates_stack,$next_state; #state
                        push @ScoreStates_stack,($score + $err_cost); #score
                    }
                    $A{TRANSITIONS}{$state}{$next_state}{"E"} = 0.33333;
                }

                # indel? (no SNP)
                if ($score + $indel_cost > $total_cost_threshold) {
                    # no indel, (no SNP)
                    $A{TRANSITIONS}{$state}{0}{"I"} = 0.33333;
                } else {
                    # indel, (no SNP)
                    my $next_state;
                    if (defined $IndelStateOfScore{($score + $indel_cost)}) {
                        $next_state = $IndelStateOfScore{($score + $indel_cost)};
                    } else {
                        $next_state = $A{SIZE};
                        $IndelStateOfScore{($score + $indel_cost)} = $next_state;
                        $A{SIZE} = $next_state + 1;
                        push @ScoreStates_stack,$next_state; #state
                        push @ScoreStates_stack,(-($score + $indel_cost)); #score
                    }
                    $A{TRANSITIONS}{$state}{$next_state}{"I"} = 0.33333;
                }

            } else {

                # (SNP)

                # (state) to (SNP_x_state)
                my $snp_x_state = $A{SIZE};
                $A{TRANSITIONS}{$state}{$snp_x_state}{"E"} = 0.33333;
                $A{SIZE} = $snp_x_state + 1;

                # (state) to (indel)
                my $next_state;
                if ($score + $indel_cost > $total_cost_threshold) {
                    # no indel, (SNP)
                    $next_state = 0;
                } else {
                    # indel, (SNP)
                    if (defined $IndelStateOfScore{($score + $indel_cost)}) {
                        $next_state = $IndelStateOfScore{($score + $indel_cost)};
                    } else {
                        $next_state = $A{SIZE};
                        $IndelStateOfScore{($score + $indel_cost)} = $next_state;
                        $A{SIZE} = $next_state + 1;
                        push @ScoreStates_stack,$next_state; #state
                        push @ScoreStates_stack,(-($score + $indel_cost)); #score
                    }
                }
                $A{TRANSITIONS}{$state}{$next_state}{"I"} = 0.33333;

                # from (SNP_x_state)
                # with "M"
                if ($score + $err_cost > $total_cost_threshold) {
                    # no error
                    $next_state = 0;
                } else {
                    # error
                    if (defined $MatchStateOfScore{($score + $err_cost)}) {
                        $next_state = $MatchStateOfScore{($score + $err_cost)};
                    } else {
                        $next_state = $A{SIZE};
                        $MatchStateOfScore{($score + $err_cost)} = $next_state;
                        $A{SIZE} = $next_state + 1;
                        push @ScoreStates_stack,$next_state; #state
                        push @ScoreStates_stack,($score + $err_cost); #score
                    }
                }
                $A{TRANSITIONS}{$snp_x_state}{$next_state}{"M"} = 0.33333;

                # with "I"
                if ($score + $err_cost + $indel_cost > $total_cost_threshold) {
                    # no indel
                    $next_state = 0;
                } else {
                    # indel
                    if (defined $IndelStateOfScore{($score + $err_cost + $indel_cost)}) {
                        $next_state = $IndelStateOfScore{($score + $err_cost + $indel_cost)};
                    } else {
                        $next_state = $A{SIZE};
                        $IndelStateOfScore{($score + $err_cost + $indel_cost)} = $next_state;
                        $A{SIZE} = $next_state + 1;
                        push @ScoreStates_stack,$next_state; #state
                        push @ScoreStates_stack,(-($score + $err_cost + $indel_cost)); #score
                    }
                }
                $A{TRANSITIONS}{$snp_x_state}{$next_state}{"I"} = 0.33333;

                # with "E"
                if (defined $MatchStateOfScore{($score + $snp_cost)}) {
                    $next_state = $MatchStateOfScore{($score + $snp_cost)};
                } else {
                    $next_state = $A{SIZE};
                    $MatchStateOfScore{($score + $snp_cost)} = $next_state;
                    $A{SIZE} = $next_state + 1;
                    push @ScoreStates_stack,$next_state; #state
                    push @ScoreStates_stack,($score + $snp_cost); #score
                }
                $A{TRANSITIONS}{$snp_x_state}{$next_state}{"E"} = 0.33333;
            }
        } else {

            # B) score is < 0 -> "Indel" state
            $score = -$score;

            my $next_state;
            if (defined $MatchStateOfScore{($score)}) {
                $next_state = $MatchStateOfScore{($score)};
            } else {
                $next_state = $A{SIZE};
                $MatchStateOfScore{($score)} = $next_state;
                $A{SIZE} = $next_state + 1;
                push @ScoreStates_stack,$next_state; #state
                push @ScoreStates_stack,($score); #score
            }

            $A{TRANSITIONS}{$state}{$next_state}{"M"} = 0.33333;
            $A{TRANSITIONS}{$state}{$next_state}{"E"} = 0.33333;

            # indel?
            if ($score + $indel_cost > $total_cost_threshold) {
                # no indel
                $A{TRANSITIONS}{$state}{0}{"I"} = 0.33333;
            } else {
                # indel
                my $next_state;
                if (defined $IndelStateOfScore{($score + $indel_cost)}) {
                    $next_state = $IndelStateOfScore{($score + $indel_cost)};
                } else {
                    $next_state = $A{SIZE};
                    $IndelStateOfScore{($score + $indel_cost)} = $next_state;
                    $A{SIZE} = $next_state + 1;
                    push @ScoreStates_stack,$next_state; #state
                    push @ScoreStates_stack,(-($score + $indel_cost)); #score
                }
                $A{TRANSITIONS}{$state}{$next_state}{"I"} = 0.33333;
            }
        }
    } # while not (stack empty)
    return %A;
}




# This function does a light cleaning of the automaton :
#  - it removes unreachable states (FIXME not done yet)
#  - it remove states with probability < $pr_threshold
sub cleaning_aa($$$) {
    my ($A_ref,$len,$pr_threshold) = @_;
    my %A = %$A_ref;
    my $A_size = $A{SIZE};
    my @P_old = ();
    my @P_new = ();
    my @P_max = ();

    # init P_old
    for(my $q=0; $q<$A_size; $q++) {
        if (defined $A{INITPROB}{$q}) {
            $P_old[$q] = $A{INITPROB}{$q};
        } else {
            $P_old[$q] = 0.0;
        }
        $P_max[$q] = $P_old[$q];
    }

    # compute P_new according to P_old
    for(my $l=0; $l<$len; $l++) {
        # erase P_new
        for(my $i=0; $i<$A_size; $i++) {
            $P_new[$i] = 0;
        }

        # recompute P_new
        for(my $q=0; $q<$A_size; $q++) {
            if ($P_old[$q] > 0.0) {
                foreach my $s_prim (keys(%{$A{TRANSITIONS}{$q}})) {
                    my $s = $s_prim;
                    foreach my $l_prim ("I","E","M") {
                        my $l = $l_prim;
                        if ((defined $A{TRANSITIONS}{$q}{$s}{$l}) && ($A{TRANSITIONS}{$q}{$s}{$l} > 0)) {
                            $P_new[$s] += $P_old[$q] * $A{TRANSITIONS}{$q}{$s}{$l};
                        }
                    }
                }
            }
        }
        # and keep the max in the P_max (and copy new in old)
        for(my $q=0; $q<$A_size; $q++) {
            if ($P_new[$q] > $P_max[$q]) {
                $P_max[$q] = $P_new[$q];
            }
            $P_old[$q] = $P_new[$q];
        }
    }


    #
    for(my $q=0; $q<$A_size; $q++) {
        foreach my $s_prim (keys(%{$A{TRANSITIONS}{$q}})) {
            my $s = $s_prim;
            if ($P_max[$s] <= $pr_threshold) {
                foreach my $l_prim ("I","E","M") {
                    my $l = $l_prim;
                    if (defined $A{TRANSITIONS}{$q}{$s}{$l}) {
                        my $pr = $A{TRANSITIONS}{$q}{$s}{$l};
                        delete $A{TRANSITIONS}{$q}{$s}{$l};
                        #FIXME (set to "0" now)
                        if (defined $A{TRANSITIONS}{$q}{($A_size-1)}{$l}) {
                            $A{TRANSITIONS}{$q}{($A_size-1)}{$l} += $pr;
                        } else {
                            $A{TRANSITIONS}{$q}{($A_size-1)}{$l} = $pr;
                        }
                    }
                }
            }
        }
    }
}



# This function print an alignment of size $len generated by the given automaton $A_ref
sub gen_seq($$) {
    my ($len,$A_ref) = @_;
    my %A = %$A_ref;

    my $p = rand();
    my $psum = 0.0;
    my $q = 0;
    for ($q=0; $q<$A{SIZE}; $q++) {
        $psum += $A{INITPROB}{$q};
        if ($psum >= $p) {
            goto eol1;
        }
    }
  eol1:;
    my $s = 0;
    my $l = "X";
    for (my $i=0;$i<$len;$i++) {
        $psum = 0.0;
        $p = rand();
        foreach my $s_prim (keys(%{$A{TRANSITIONS}{$q}})) {
            $s = $s_prim;
            foreach my $l_prim ("I","E","M") {
                $l = $l_prim;
                if ((defined $A{TRANSITIONS}{$q}{$s}{$l})) {
                    $psum += $A{TRANSITIONS}{$q}{$s}{$l};
                    if ($psum >= $p) {
                        goto eol2;
                    }
                }
            }
        }
      eol2:
        print $l;
        $q = $s;
    }
    print "\n";
}


# 1) lossless automaton
# my %L = &lossless_aa(3,#$snp_cost
#                      2,#$err_cost
#                      4,#$indel_cost,
#                      7);#,$total_cost_threshold
#
# &print_aa_iedera_raw(%L);
# exit();


# 2) lossy automaton
my %A = &build_error_complex_aa(34,
                                0.01,#prstart
                                0.10,#prend
                                0.00,#indelrate
                                0.02,#e_freq
                                0.25,#e_err
    );
my %C = &product_aa(\%A,\%SNP_AA);
#&print_aa_c(%C);
&print_aa_iedera_prob(%C);
exit();

# # 3) generate sequences
# while (1) {
#     &gen_seq(34,\%C);
# }
