1/*
    2Hidden Markov model for modeling DNA sequences.
    3The model has three states, q1, q2 and end, and four output symbols,
    4a, c, g, and t, corresponding to the four nucleotides (letter).
    5From
    6Christiansen, H. and Gallagher, J. P. 2009. Non-discriminating arguments and
    7their uses. In International Conference on Logic Programming. LNCS,
    8vol. 5649. Springer, 55-69.
    9
   10*/
   11:- use_module(library(viterbi)).   12
   13:- if(current_predicate(use_rendering/1)).   14:- use_rendering(c3).   15:- use_rendering(graphviz).   16:- endif.   17
   18:- viterbi.   19
   20:- begin_lpad.   21
   22% hmm(O): O is the output sequence
   23% hmm1(S,O): O is the output sequence and S is the sequence of states
   24% hmm(Q,S0,S,O):  from state Q and previous state S0, generates output O and
   25% sequence of states S
   26
   27hmm(O):-hmm1(_,O).
   28% O is an output sequence if there is a state seuqnece S such that hmm1(S,O)
   29% holds
   30
   31hmm1(S,O):-hmm(q1,[],S,O).
   32% O is an output sequence and S a state sequence if the chain stats at state
   33% q1 and ends generating state sequence S and output sequence O
   34
   35hmm(end,S,S,[]).
   36% an HMM in state end terminates the sequence without emitting any symbol
   37
   38hmm(Q,S0,S,[L|O]):-
   39	Q\= end,
   40	next_state(Q,Q1,S0),
   41	letter(Q,L,S0),
   42	hmm(Q1,[Q|S0],S,O).
   43% an HMM in state Q different from end goes in state Q1, emits the letter L
   44% and continues the chain
   45
   46next_state(q1,q1,S):1/2;next_state(q1,q2,S):0.45;next_state(q1,end,S):0.05.
   47% from state q1 the HMM can go to q1, q2 or end with equal probability
   48
   49next_state(q2,q1,S):0.45;next_state(q2,q2,S):1/2;next_state(q2,end,S):0.05.
   50% from state q2 the HMM can go to q1, q2 or end with equal probability
   51
   52
   53letter(q1,a,S):0.4;letter(q1,c,S):0.3;letter(q1,g,S):0.2;letter(q1,t,S):0.1.
   54% from state q1 the HMM emits one of the letters with equal probability
   55
   56letter(q2,a,S):0.1;letter(q2,c,S):0.2;letter(q2,g,S):0.3;letter(q2,t,S):0.4.
   57% from state q1 the HMM emits one of the letters with equal probability
   58
   59:- end_lpad.   60
   61
   62state_diagram(digraph(G)):-
   63    findall(edge(A -> B,[label=P]),
   64      (clause(next_state(A,B,_,_,_),
   65        (get_var_n(_,_,_,_,Probs,_),equality(_,_,N,_))),
   66        nth0(N,Probs,P)),
   67      G).