1/*
    2Model checking of a Markov chain: we want to know what is the likelihood 
    3that on an execution of the chain from a start state s, a final state t 
    4will be reached?
    5From
    6Gorlin, Andrey, C. R. Ramakrishnan, and Scott A. Smolka. "Model checking with probabilistic tabled logic programming." Theory and Practice of Logic Programming 12.4-5 (2012): 681-700.
    7*/
    8
    9
   10:- use_module(library(mcintyre)).   11
   12:- if(current_predicate(use_rendering/1)).   13:- use_rendering(c3).   14:- use_rendering(graphviz).   15:- endif.   16
   17:- mc.   18
   19:- begin_lpad.   20% reach(S, I, T) starting at state S at instance I,
   21%   state T is reachable.
   22reach(S, I, T) :-
   23  trans(S, I, U),
   24  reach(U, next(I), T).
   25
   26reach(S, _, S).
   27
   28% trans(S,I,T) transition from S at instance I goes to T
   29
   30trans(s0,S,s0):0.5; trans(s0,S,s1):0.3; trans(s0,S,s2):0.2.
   31
   32trans(s1,S,s1):0.4; trans(s1,S,s3):0.1; trans(s1,S,s4):0.5.
   33
   34trans(s4,_,s3).
   35:- end_lpad.   36
   37markov_chain(digraph(G)):-
   38    findall(edge(A -> B,[label=P]),
   39      (clause(trans(A,_,B),
   40        (sample_head(_,_,_,Probs,N))),
   41        nth0(N,Probs,_:P)),
   42      G0),
   43    findall(edge(A -> B,[label=1.0]),
   44      clause(trans(A,_,B),true),
   45      G1),
   46    append(G0,G1,G).

?- mc_prob(reach(s0,0,s0),P). % expected result ~ 1.

?- mc_prob(reach(s0,0,s1),P). % expected result ~ 0.5984054054054054.

?- mc_prob(reach(s0,0,s2),P). % expected result ~ 0.4025135135135135.

?- mc_prob(reach(s0,0,s3),P). % expected result ~ 0.5998378378378378.

?- mc_prob(reach(s0,0,s4),P). % expected result ~ 0.49948717948717947.

?- mc_prob(reach(s1,0,s0),P). % expected result ~ 0.

?- mc_sample(reach(s0,0,s1),1000,P,[successes(T),failures(F)]). % expected result ~ 0.5984054054054054.

?- mc_sample(reach(s0,0,s1),1000,Prob),bar(Prob,C).

?- mc_sample_arg(reach(s0,0,S),50,S,Values). % take 50 samples of L in findall(S,reach(s0,0,S),L)

?- mc_sample_arg(reach(s0,0,S),50,S,O),argbar(O,C). % take 50 samples of L in findall(S,reach(s0,0,S),L)

?- mc_sample_arg_first(reach(s0,0,S),50,S,Values). % take 50 samples of the first value returned for S in reach(s0,0,S)

?- mc_sample_arg_first(reach(s0,0,S),50,S,O),argbar(O,C). % take 50 samples of the first value returned for S in reach(s0,0,S)

?- markov_chain(G). % draw the Markov chain */