1/*
    2Model checking of the Synchronous Leader Election Protocol expressed in 
    3Probabilistic Computation Tree Logic (PCTL).
    4
    5Given a synchronous ring of N processes the protocol is used 
    6to elect a leader (a uniquely designated processor) by sending messages around 
    7the ring.
    8
    9The protocol proceeds in rounds and is parametrised by a constant K. Each round 
   10begins by all processors (independently) choosing a random number (uniformly) 
   11from {1,...,K} as an id. The processors then pass their ids around the ring. 
   12If there is a unique id, then the processor with the maximum unique id is 
   13elected the leader, and otherwise the processors begin a new round.
   14
   15With this program you can 
   16- check that the probability of eventually electing a leader is 1
   17- compute the probability of electing a leader within a certain 
   18  number of rounds
   19- compute the expected number of rounds to elect a leader
   20- graph the probability of electing a leader within L rounds as a function of L
   21- graph the expected number of rounds to elect a leader as a function of the 
   22  number of processes or of K
   23From
   24Gorlin, 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.
   25This program was kindly provided by Andrey Gorlin and adapted to MCINTYRE
   26by Fabrizio Riguzzi.
   27See also http://www.prismmodelchecker.org/casestudies/synchronous_leader.php
   28*/

% see http://www.prismmodelchecker.org/casestudies/synchronous_leader.php % What is the probability that eventually a leader is elected? ?- mc_sample(eventually(elect),100,P). % expected result 1

% What is the probability of electing a leader within 3 rounds? ?- mc_sample(bounded_eventually(elect,3),1000,P). % expected result 0.97

% What is the expected number of rounds to elect a leader? ?- mc_expectation(eventually(elect,T),1000,T,E). % expected result 1.2

?- graph_prob(G). % graph the probability of electing a leader within L rounds as % a function of L

?- graph_exp_rounds_n(G). % graph the expected number of rounds to elect a leader as a % funtion of the number of processes when K=3

?- graph_exp_rounds_k(G). % graph the expected number of rounds to elect a leader as a % funtion of K when N=3

?- network_topology(G). % draw a graph representing the network topology */

   60:- use_module(library(mcintyre)).   61
   62:- if(current_predicate(use_rendering/1)).   63:- use_rendering(c3).   64:- use_rendering(graphviz).   65:- endif.   66:- dynamic kr/1,num/1.   67:- mc.   68
   69:- begin_lpad.   70
   71% State Formulae 
   72models(_S, tt,_Hist,_Limit,_Time).
   73models(S, prop(P),_Hist,_Limit,_Time) :-
   74	proposition(P, S).
   75models(S, and(F1, F2),Hist,Limit,Time) :-
   76	models(S, F1,Hist,Limit,Time), models(S, F2,Hist,Limit,Time).
   77models(S, or(F1, _F2),Hist,Limit,Time) :-
   78	models(S, F1,Hist,Limit,Time).
   79models(S, or(F1, F2),Hist,Limit,Time) :-
   80	\+ models(S, F1,Hist,Limit,Time),
   81	models(S, F2,Hist,Limit,Time).
   82models(S, not(F), Hist,Limit,Time) :-
   83	\+ models(S, F,Hist,Limit,Time).
   84models(S, prob_until(comp(Op, P), F1, F2),Hist,Limit,Time) :-
   85	mc_sample(pmodels(S, until(F1, F2),Hist,Limit,Time),20, Q),
   86	comp(Q, Op, P).
   87models(S, prob_next(comp(Op, P), F),Hist,Limit,Time) :-
   88	mc_sample(pmodels(S, next(F),Hist,Limit,Time),20, Q),
   89	comp(Q, Op, P).
   90
   91comp(Q,>,P):-
   92  Q>P.
   93
   94comp(Q,>=,P):-
   95  Q>=P.
   96
   97comp(Q,<,P):-
   98  Q<P.
   99
  100comp(Q,=<,P):-
  101  Q=<P.
  102
  103
  104% Path Formulae
  105pmodels(S,F):-
  106  pmodels(S,F,[],nolimit,0,_Time).
  107
  108pmodels(S,F,Hist,Limit,Time):-
  109  pmodels(S,F,Hist,Limit,Time,_Time).
  110
  111pmodels(S, until(_F1, F2),Hist,Limit,Time,Time) :-
  112	models(S, F2,Hist,Limit,Time),!.
  113	
  114pmodels(S, until(F1, F2),Hist0,Limit,Time0,Time) :-
  115	within_limit(Time0,Limit),
  116	models(S, F1,Hist0,Limit,Time0),
  117	ctrans(S, _, T, Hist0, Hist),!,
  118	Time1 is Time0+1,
  119	pmodels(T, until(F1,F2),Hist,Limit,Time1,Time).
  120
  121pmodels(S, next(F), Hist,Limit,Time0,Time) :-
  122	within_limit(Time0,Limit),
  123	ctrans(S, _, T, Hist, _),!,
  124	Time is Time0+1,
  125	models(T, F,Hist,Limit,Time).
  126
  127within_limit(_Time,nolimit):-!.
  128
  129within_limit(Time,Limit):-
  130  Time<Limit.
  131
  132bounded_eventually(Prop,Rounds):-
  133  num(N),
  134  B is Rounds*(N+1),
  135  eventually(Prop,B,_T).
  136
  137eventually(Prop):-
  138  eventually(Prop,_T).
  139
  140eventually(Prop,Rounds):-
  141  eventually(Prop,nolimit,T),
  142  num(N),
  143  Rounds is T/(N+1).
  144
  145eventually(Prop,Limit,T) :-
  146	init(S),
  147	pctlspec(Prop, F),
  148	pmodels(S, F,[],Limit,0,T).
  149
  150
  151pctlspec(X, until(tt, prop(X))).
  152proposition(P, S) :- final(P, S).
  153
  154final(elect, [_|L]) :-
  155	num(N),
  156	gen_elected_state(N, L).
  157
  158gen_elected_state(J, L) :-
  159	(J==0
  160	->    L=[]
  161	;     L = [state(3,_,_,_)|Rest],
  162	      J1 is J-1,
  163	      gen_elected_state(J1,Rest)
  164	).
  165	
  166
  167% transitions
  168% module counter
  169% [read] c<N-1 -> (c'=c+1);
  170% reading
  171trans(counter, counter(C), read, counter(D),_S,H,H) :-
  172  num(N),
  173  C < N-1,
  174  D is C+1.
  175
  176% [read] c=N-1 -> (c'=c);
  177% finished reading
  178trans(counter, counter(C), read, counter(C),_S,H,H) :-
  179  num(N),
  180  C =:= N-1.
  181
  182% [done] u1 | u2 | u3 | u4 -> (c'=c);
  183% done
  184trans(counter, counter(C), done, counter(C),S,H,H) :-
  185  get_processid(P), 
  186  nonlocal(process(P,_), uniqueid, 1,S).
  187   
  188
  189% [retry] !(u1 | u2 | u3 | u4) -> (c'=1);
  190% pick again reset counter 
  191trans(counter, counter(_C), retry, counter(1),S,H,H) :-
  192        findall(P,get_processid(P),PL),
  193	maplist(nl(S),PL).
  194
  195% [loop] s1=3 -> (c'=c);
  196% loop (when finished to avoid deadlocks)
  197trans(counter, counter(C), loop, counter(C),S,H,H) :-
  198  nonlocal(process(1,_), state, 3,S).
  199
  200% module process
  201% local state
  202% s1=0 make random choice
  203% s1=1 reading
  204% s1=2 deciding
  205% s1=3 finished
  206
  207% [pick] s1=0 -> 1/K : (s1'=1) & (p1'=0) & (v1'=0) & (u1'=true) + ...;
  208% pick value
  209trans(process(_N,_Next), state(0,_,_,_), pick, state(1,1,R,R),_S,H,[pick(R)|H]) :-
  210  pick(H,R).
  211
  212%read 
  213% [read] s1=1 &  u1 & !p1=v2 & c<N-1 -> (u1'=true) & (v1'=v2);
  214trans(process(_N,Next), state(1,1,_,P), read, state(1,1,V,P),S,H,H) :-
  215  nonlocal(counter, counter, C,S),
  216  num(CN),
  217  C < CN - 1,
  218  nonlocal(process(Next,_), value, V,S),
  219  P \== V.
  220
  221% [read] s1=1 &  u1 &  p1=v2 & c<N-1 -> (u1'=false) & (v1'=v2) & (p1'=0);
  222trans(process(_N,Next), state(1,1,_,P), read, state(1,0,V,0),S,H,H) :-
  223  nonlocal(counter, counter, C,S),
  224  num(CN),
  225  C < CN - 1,
  226  nonlocal(process(Next,_), value, V,S),
  227  P == V.
  228
  229% [read] s1=1 & !u1 &  c<N-1 -> (u1'=false) & (v1'=v2);
  230trans(process(_N,Next), state(1,0,_,P), read, state(1,0,V,P),S,H,H) :-
  231  nonlocal(counter, counter, C,S),
  232  num(CN),
  233  C < CN - 1,
  234  nonlocal(process(Next,_), value, V,S).
  235 
  236% read and move to decide 
  237% [read] s1=1 &  u1 & !p1=v2 & c=N-1 -> (s1'=2) & (u1'=true) & (v1'=0) & (p1'=0);
  238trans(process(_N,Next), state(1,1,_,P), read, state(2,1,0,0),S,H,H) :-
  239  nonlocal(counter, counter, C,S),
  240  num(CN),
  241  C =:= CN - 1,
  242  nonlocal(process(Next,_), value, V,S),
  243  P \== V.
  244
  245% [read] s1=1 &  u1 &  p1=v2 & c=N-1 -> (u1'=false) & (v1'=0) & (p1'=0);
  246trans(process(_N,Next), state(1,1,_,P), read, state(2,0,0,0),S,H,H) :-
  247  nonlocal(counter, counter, C,S),
  248  num(CN),
  249  C =:= CN - 1,
  250  nonlocal(process(Next,_), value, V,S),
  251  P == V.
  252
  253% [read] s1=1 & !u1 &  c=N-1 -> (s1'=2) & (u1'=false) & (v1'=0);
  254trans(process(_N,_Next), state(1,0,_,P), read, state(2,0,0,P),S,H,H) :-
  255  nonlocal(counter, counter, C,S),
  256  num(CN),
  257  C =:= CN - 1.
  258
  259% done
  260% [done] s1=2 -> (s1'=3) & (u1'=false) & (v1'=0) & (p1'=0);
  261trans(process(_N,_Next), state(2,_,_,_), done, state(3,0,0,0),_S,H,H).
  262
  263% retry
  264% [retry] s1=2 -> (s1'=0) & (u1'=false) & (v1'=0) & (p1'=0);
  265trans(process(_N,_Next), state(2,_,_,_), retry, state(0,0,0,0),_S,H,H).
  266
  267% loop (when finished to avoid deadlocks)
  268% [loop] s1=3 -> (s1'=3);
  269trans(process(_N,_Next), state(3,U,V,P), loop, state(3,U,V,P),_S,H,H).
  270
  271pick(H,V):-
  272  kr(K),
  273  K1 is K-1,
  274  PH is 1/K,
  275  findall(I,between(0,K1,I),L),
  276  foldl(pick_value(H,PH),L,(1,_),(_,V)).
  277
  278pick_value(_H,_PH,_I,(P0,V0),(P0,V0)):-
  279  nonvar(V0).
  280
  281pick_value(H,PH,I,(P0,V0),(P1,V1)):-
  282  var(V0),
  283  PF is PH/P0,
  284  (pick_fact(H,V0,PF)->
  285    P1=PF,
  286    V1=I
  287  ;
  288    P1 is P0*(1-PF),
  289    V1=V0
  290  ).
  291
  292pick_fact(_,_,P):P.
  293
  294%pick(H,0):0.5; pick(H,1):0.5.
  295
  296ctrans(S, A, T, Hi, Ho) :-
  297	config(P),
  298	find_matching_trans(P,S,S,[],A,T,Hi,Ho).
  299
  300find_matching_trans([], [], _CS, _PA, A, [], H,H) :- nonvar(A).
  301find_matching_trans([P|Ps], [S|Ss], CS, PA, A, [T|Ts],Hi,Ho) :-
  302	pick_trans(P, S, CS, PA, A, T, Hi,H1),
  303	find_matching_trans(Ps, Ss, CS, PA, A, Ts,H1,Ho).
  304find_matching_trans([P|Ps], [S|Ss], CS, PA, A, [S|Ts], Hi,Ho) :-
  305	% skip current process; but then all transitions enabled in the current
  306	% process will be prohibited for selection in later processes.
  307	enabled_trans(P,L),
  308	(nonvar(A) -> \+ member(A,L); true),
  309	append(L, PA, PA1),
  310	find_matching_trans(Ps, Ss, CS, PA1, A, Ts, Hi, Ho).
  311
  312pick_trans(P, S, CS, PA, A, T, H0,H) :-
  313	etrans(P, S, PA, A, T,CS, H0,H).
  314
  315etrans(P, S, PA, A, T, CS,H0,H) :-
  316	trans(P, S, A, T,CS,H0,H),
  317	A \= epsilon,
  318	\+ member(A, PA).
  319
  320enabled_trans(P, L) :-
  321	setof(A, enabled_trans_in_process(P,A), L).
  322
  323enabled_trans_in_process(P,A) :-
  324	clause(trans(P,_,A,_,_,_,_),_),
  325	A \= epsilon.
  326
  327nonlocal(Proc, Var, Val,CS) :-
  328	getpid(Proc, Var, Pid, Idx),
  329	nth1(Pid, CS, State),
  330	arg(Idx, State, Val).
  331%	writeln(nonlocal_read(Proc, State, Var, Val)).
  332
  333getpid(Proc, Var, Pid, Idx) :-
  334	config(Config),
  335	nth1(Pid, Config, Proc),
  336	nonlocal_access(Proc, Var, Idx).
  337
  338get_processid(P):-
  339  num(N),
  340  between(1,N,P).
  341
  342config([counter| L]) :-
  343	findall(P,get_processid(P),PL),
  344	maplist(neighbor,PL,L).
  345
  346neighbor(I,process(I,J)) :-
  347	num(N),
  348	(I<N
  349	->  J is I+1
  350	;   J = 1
  351	).
  352
  353%config([counter, process(1,2), process(2,3), process(3,4), process(4,1)]).
  354
  355init(S) :-
  356	config(P),
  357	maplist(init_state,P,S).
  358
  359init_state(counter, counter(1)).
  360init_state(process(_,_), state(0,0,0,0)).
  361
  362nonlocal_access(counter, counter, 1).
  363nonlocal_access(process(_,_), state, 1).
  364nonlocal_access(process(_,_), uniqueid, 2).
  365nonlocal_access(process(_,_), value, 3).
  366
  367nl(S,P):-
  368  nonlocal(process(P, _), uniqueid, 0,S).
  369
  370num(4).
  371kr(4).
  372
  373
  374:- end_lpad.  375
  376network_topology(digraph(G)):-
  377  config([_|L]),
  378  maplist(to_edge,L,G).
  379
  380to_edge(process(A,B),edge(A -> B,[])).
  381
  382graph_prob(G):-
  383  retract(num(N)),
  384  retract(kr(K)),
  385  assert(num(4)),
  386  assert(kr(2)),
  387  findall(L-P,
  388    (between(1,6,L),mc_sample(bounded_eventually(elect,L),100,P)),LV),
  389  G=c3{data:_{x:x, rows:[x-'Probability of leader elected within L rounds (N=4, K=2)'|LV]},%legend:_{show: false},
  390    axis:_{x:_{min:1,max:6,label:'L',padding:_{bottom:0.0,top:0.0}},
  391%        tick:_{values:[0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]}},
  392           y:_{label:'Probability',padding:_{bottom:0.0,top:0.0}}}},
  393%        tick:_{values:[0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]}}}},
  394  retract(num(4)),
  395  retract(kr(2)),
  396  assert(num(N)),
  397  assert(kr(K)).
  398
  399graph_exp_rounds_n(G):-
  400  retract(num(NI)),
  401  retract(kr(KI)),
  402  assert(kr(3)),
  403  findall(N-E,
  404    (between(3,8,N),
  405     assert(num(N)),
  406     mc_expectation(eventually(elect,T),100,T,E),
  407     retract(num(N))),
  408    LV),
  409  G=c3{data:_{x:x, rows:[x-'Expected rounds to elect a leader (K=3)'|LV]},%legend:_{show: false},
  410    axis:_{x:_{min:3,max:8,label:'N',padding:_{bottom:0.0,top:0.0}},
  411           y:_{label:'Expected rounds',padding:_{bottom:0.0,top:0.0}}}},
  412  retract(kr(3)),
  413  assert(num(NI)),
  414  assert(kr(KI)).
  415
  416graph_exp_rounds_k(G):-
  417  retract(num(NI)),
  418  retract(kr(KI)),
  419  assert(num(3)),
  420  findall(K-E,
  421    (between(3,8,K),
  422     assert(kr(K)),
  423     mc_expectation(eventually(elect,T),500,T,E),
  424     retract(kr(K))),
  425    LV),
  426  G=c3{data:_{x:x, rows:[x-'Expected rounds to elect a leader (N=3)'|LV]},%legend:_{show: false},
  427    axis:_{x:_{min:3,max:8,label:'K',padding:_{bottom:0.0,top:0.0}},
  428           y:_{label:'Expected rounds',padding:_{bottom:0.0,top:0.0}}}},
  429  retract(num(3)),
  430  assert(num(NI)),
  431  assert(kr(KI))