######################################
# we need numeric
# 
import Numeric
######################################
def viterbi_posterior(hmm, apos):
    ''' viterbi algorithm on the posterior probability
        apos[L][N]= matrix of posterior values (a Numeric.Float64 array)
                    L length of the sequence + 1
                    N the number of the hmm states   
        hmm = a hidden markov model object, 
              hmm.num_states = number of state 
              hmm.emits = list of the emitting states
              hmm.nulls = list of the silent states
              hmm.topo_order = list of all state including begin
                               in topological order
              hmm.in_s = list of inlinks.  hmm.in_s[s] is the list of the s inlinks
              hmm.in_s_e = list of emitting inlinks for each state (like hmm.in_s)
              hmm.in_s_n = list of silet inlinks for each state (like hmm.in_s)
              hmm.end_s = list of the states that can end
        -> 
           return(best_state_path,val)
           best_state_path = the best path
           val = the score of the best path
    '''
    # settings
    big_negative=-10000.0
    L=len(seq)+1 # the emissions start with index 1
    # lf keeps the values of each sequence/state position
    # as in the viterbi
    lf=Numeric.array([[big_negative]*hmm.num_states]*L,Numeric.Float64)
    # bkt is the backtrace matrix
    bkt=Numeric.array([[-1]*hmm.num_states]*L,Numeric.Int)
    # nicknames
    E=hmm.emits # the list of the emitting states
    N=hmm.nulls # the list of the silent states 
    S=hmm.topo_order # the list of all state including B
    IS=hmm.in_s # inlinks 
    IE=hmm.in_s_e # inlinks from emittings only
    IN=hmm.in_s_n # inlinks from silent only
    ENDS=hmm.end_s # end states
    B=0 # begin state
    #
    ###### START PHASE
    lf[0][B]=0.0
    bkt[0][B]=-1 # no state to backtrace
    for s in E:
        lf[0][s]=big_negative
        bkt[0][s]=B
    # from N+B -> N | ARE all equivalent to Begin
    # we are assuming a null can have only a null inlink
    for s in N:
        for t in IN[s]: # IN[s] inludes also B (if this is possible)  
            if lf[0][t] > lf[0][s]: # there exitsts some inlinks or B
                lf[0][s]=lf[0][t]
                bkt[0][s]=t
    ###### RECURRENCE PHASE for 1 to L-1
    for i in range(1,L):
        # S -> E
        for s in E:
            # update done only if labels are not used
            # or labels corresponds (labels[i-1] is position i !!)
            for t in IS[s]:
                ltmp=lf[i-1][t] 
                if(lf[i][s]<ltmp):
                    lf[i][s]=ltmp
                    bkt[i][s]=t
            lf[i][s]+=Numeric.log(apos[i,s]) # 
        # E -> N . Note the i-level is the same
        for s in N:
            for t in IE[s]:
                ltmp=lf[i][t] 
                if(lf[i][s]<ltmp):
                    lf[i][s]=ltmp
                    bkt[i][s]=t
        # N -> N . Note the i-level is the same
        for s in N:
            for t in IN[s]:
                ltmp=lf[i][t] # ln_a(t,s)+lf[i][t]
                if(lf[i][s]<ltmp):
                    lf[i][s]=ltmp
                    bkt[i][s]=t
    ###### END PHASE 
    bestval=big_negative
    pback=None # back pointer
    for s in ENDS:
        if(bestval<lf[L-1][s]):
            bestval=lf[L-1][s]
            pback=s
    if not pback:
        print "ERRORRRRR!!!"
    ###### BACKTRACE
    best_path=[]
    i=L-1
    while i >= 0 and pback != B:
        best_path.insert(0,pback)
        ptmp=bkt[i][pback]
        if pback in hmm.emits:
           i-=1
        pback=ptmp

    return(best_path,bestval)
######################################