1		% (c) 2025-2026 Lehrstuhl fuer Softwaretechnik und Programmiersprachen,
2		% Heinrich Heine Universitaet Duesseldorf
3		% This software is licenced under EPL 1.0 (http://www.eclipse.org/org/documents/epl-v10.html)
4
5		%###############################################
6
7		% This is a model-checker verifying Probabilistic Computation
8		% Tree Logic (PCTL) formulas over Discrete-Time Markov Chains (DTMC)
9		% written by Honore Marmion
10
11		% Working on SICStus prolog 4.10
12
13		%###############################################
14
15		:- module(dtmc_model_checking,[pctl_model_check/4
16		% sat/1, sat/2,
17		% search_prob0/4,prob0/4,prob1/4,table_prob0/3,state/1
18		]).
19		:- use_module(library(clpr),[{}/1]).
20
21
22		:- use_module(probsrc(error_manager),[add_error/3, add_internal_error/2, add_warning/3]).
23		:- use_module(probltlsrc(ltl_tools),[temporal_parser/3]).
24		:- use_module(probsrc(state_space),[find_initialised_states/1, current_state_id/1]).
25		:- use_module(probsrc(tools),[start_ms_timer/1, stop_ms_timer_with_msg/2]).
26
27		pctl_model_check(Formula,_MaxNodes,Mode,Res) :-
28		(temporal_parser(Formula,pctl,PCtlFormula) -> true ; add_error(ctl,'PCTL Parser failed: ',Formula),fail),
29		%set_max_nr_of_new_impl_trans_nodes(MaxNodes),
30		% TODO: provide better feedback, extract probability for P={p} formulas
31		(pre_process_formula(PCtlFormula,ProcessedFormula,Bindings,[]) -> true
32		; add_error(dtmc_model_checking,'Could not pre-process: ',PCtlFormula),
33		ProcessedFormula = PCtlFormula
34		),
35		(Mode= specific_node(ID) -> Start=ID
36		; Mode = starthere -> current_state_id(Start)
37		; Mode = init -> find_initialised_states(SN), member(Start,SN) % TODO try out all start nodes??
38		; add_error(dtmc_model_checking,'Illegal starting mode (init, starthere supported):',Mode),
39		Start = root
40		),
41		format('Checking PCTL formula from state ~w : ~w~n AST: ~w~n Open Probability Variables: ~w~n',[Start,Formula,ProcessedFormula,Bindings]),
42		reset_tables,
43		start_ms_timer(T1),
44		(sat(ProcessedFormula,Start)
45		-> (Bindings=[] -> Res=true ; Res = solution(Bindings))
46		; Res = false),
47		stop_ms_timer_with_msg(T1,'PTCTL model checking').
48
49		reset_tables :- retractall(table_prob0(_,_,_)),retractall(node(_)), retractall(prob_current(_,_)).
50
51		pre_process_formula(true,true) --> [].
52		pre_process_formula(false,false) --> [].
53		pre_process_formula(ap(P),ap(P)) --> [].
54		pre_process_formula(p(P),p(P)) --> [].
55		pre_process_formula(probformula(Operator,RefProbAtom,Ctl_formula),
56		probformula(Operator,RefProbNr,Ctl_formula)) -->
57		({ground_number(RefProbAtom,RefProbNr)} -> []
58		; {Operator=equal} -> [RefProbAtom/RefProbNr]
59		; {add_error(dtmc_model_checking,'Illegal reference probability for operator (only equal is supported for symbolic values):',RefProbAtom/Operator)},
60		[RefProbAtom/RefProbNr]
61		). % TODO: pre_process_path_formula
62		pre_process_formula(and(F,G),and(PF,PG)) -->
63		pre_process_formula(F,PF),
64		pre_process_formula(G,PG).
65		pre_process_formula(or(F,G),and(PF,PG)) -->
66		pre_process_formula(F,PF),
67		pre_process_formula(G,PG).
68		pre_process_formula(implies(F,G),and(PF,PG)) -->
69		pre_process_formula(F,PF),
70		pre_process_formula(G,PG).
71		pre_process_formula(not(F),not(PF)) -->
72		pre_process_formula(F,PF).
73
74		%pre_process_path_formula(u(F,G),u(PF,PG)) -->
75		% pre_process_formula(F,PF),
76		% pre_process_formula(G,PG).
77
78
79		% Interface predicates defining the system to be checked
80
81		:- use_module(probsrc(xtl_interface),[xtl_transition/4]).
82		:- use_module(probltlsrc(ltl_propositions), [trans/4, check_ap/2]). % check_transition_pred/5
83		:- use_module(probsrc(state_space),[visited_expression/2, visited_expression_id/1]).
84
85		% TODO: move to specfile and extend to work with B pragmas,...
86		transition_probability(From,To,Prob) :- trans(From,To,_,_),
87		visited_expression(From,FS),
88		visited_expression(To,TS),
89		xtl_transition(FS,_,TS,Infos),
90		(member(probability/X,Infos) -> Prob=X ; write(no_probability(From,To,Infos)),nl,fail).
91
92		state(X) :- visited_expression_id(X).
93
94		% example public_examples/XTL/markov/SimpleMarkov.P
95		% ?- use_module(extension('markov/dtmc_model_checking.pl')).
96		% ?- sat(prob_formula(eq,P,f(p(xtl_predicate_check(finished))))).
97
98		%PCTL Model-checking
99
100		% This term allow to handle nested formulas in case of dynamic model-checking
101		:- dynamic node/1.
102
103		%start(root).
104		start(ID) :- find_initialised_states(I), member(ID,I).
105
106		:- public sat/1.
107		sat(Formula) :- start(E), sat(Formula,E).
108
109
110		% Classic cases
111
112		sat_not(false,_E).
113		sat_not(true,_E) :- fail.
114		sat_not(p(Property),E) :- \+(check_ap(Property,E)).
115		sat_not(ap(Property),E) :- \+(check_ap(Property,E)).
116		sat_not(and(F,_G),E) :- sat_not(F,E).
117		sat_not(and(_F,G),E) :- sat_not(G,E).
118		sat_not(or(F,G),E) :- sat_not(F,E),sat_not(G,E).
119		sat_not(implies(F,G),E) :- sat(not(or(not(F),G)),E).
120		sat_not(equivalence(F,G),E) :- sat_not(and(implies(F,G),implies(G,F)),E).
121		sat_not(not(F),E) :- sat(F,E).
122
123		%sat(Formula,State) :- write(sat(Formula,State)),nl,fail.
124		sat(true,_E).
125		sat(false,_E):-fail.
126		sat(ap(Property),E) :- check_ap(Property,E).
127		sat(p(Property),E) :- check_ap(Property,E).
128		sat(and(F,G),E) :- sat(F,E), sat(G,E).
129		sat(or(F,_G),E) :- sat(F,E).
130		sat(or(_F,G),E) :- sat(G,E).
131		sat(implies(F,G),E) :- sat(or(not(F),G),E).
132		sat(equivalence(F,G),E) :- sat(and(implies(F,G),implies(G,F)),E).
133		sat(not(F),E) :- probformula(_,_,_)\=F, sat_not(F,E).
134
135		% Probabilistic-formula cases. Operator is =, < >,<= or >=. P is a number between 0 and 1 (or a Variable),
136		% E is a state.
137		% Check if the formula is nested
138
139		sat(probformula(Operator,P,Ctl_formula),E) :-
140		Ctl_formula = gk(_,_) ->
141		sat_gk(probformula(Operator,P,Ctl_formula),E)
142		; Ctl_formula = g(_) ->
143		sat_g(probformula(Operator,P,Ctl_formula),E)
144		; (node(Node) ->
145		New_Node is Node +1,
146		retract(node(Node)),
147		assert(node(New_Node)),
148		sat_node(probformula(Operator,P,Ctl_formula),E,New_Node)
149		; assert(node(0)),
150		sat_node(probformula(Operator,P,Ctl_formula),E,0),
151		retractall(node(_)) /*reinitialize nodes*/
152		).
153
154		sat(not(probformula(Operator,P,Ctl_formula)),E) :-
155		negate_operator(Operator,NotOp),
156		sat(probformula(NotOp,P,Ctl_formula),E).
157
158		% Always bounded formula, we use the dual probabilistic event of fk(K,not(F))
159		sat(probformula(Operator,P,gk(K,F)),E) :-
160		ground(P) ->
161		Q is 1-P,
162		(Operator = equal ->
163		sat(probformula(Operator,Q,fk(K,not(F))),E)
164		; sat(not(probformula(Operator,Q,fk(K,not(F)))),E))
165		; ((Operator = equal ->
166		sat(probformula(Operator,Q,fk(K,not(F))),E)
167		; sat(not(probformula(Operator,Q,fk(K,not(F)))),E)),
168		P is 1-Q)
169		.
170
171		% Always bounded formula, we use the dual probabilistic event of fk(K,not(F))
172		sat_gk(probformula(Operator,P,gk(K,F)),E) :-
173		ground(P) ->
174		Q is 1-P,
175		(Operator = equal ->
176		sat(probformula(Operator,Q,fk(K,not(F))),E)
177		; Operator = greater ->
178		sat(probformula(less,Q,fk(K,not(F))),E)
179		; Operator = less ->
180		sat(probformula(greater,Q,fk(K,not(F))),E)
181		; Operator = strictlygreater ->
182		sat(probformula(strictlyless,Q,fk(K,not(F))),E)
183		; Operator = strictlyless ->
184		sat(probformula(strictlygreater,Q,fk(K,not(F))),E))
185		; ((Operator = equal ->
186		sat(probformula(Operator,Q,fk(K,not(F))),E)
187		; Operator = greater ->
188		sat(probformula(less,Q,fk(K,not(F))),E)
189		; Operator = less ->
190		sat(probformula(greater,Q,fk(K,not(F))),E)
191		; Operator = strictlygreater ->
192		sat(probformula(strictlyless,Q,fk(K,not(F))),E)
193		; Operator = strictlyless ->
194		sat(probformula(strictlygreater,Q,fk(K,not(F))),E)),
195		P is 1-Q)
196		.
197
198		% Always formula
199		sat_g(probformula(Operator,P,g(F)),E) :-
200		ground(P) ->
201		Q is 1-P,
202		(Operator = equal ->
203		sat(probformula(Operator,Q,f(not(F))),E)
204		; Operator = greater ->
205		sat(probformula(less,Q,f(not(F))),E)
206		; Operator = less ->
207		sat(probformula(greater,Q,f(not(F))),E)
208		; Operator = strictlygreater ->
209		sat(probformula(strictlyless,Q,f(not(F))),E)
210		; Operator = strictlyless ->
211		sat(probformula(strictlygreater,Q,f(not(F))),E))
212		; ((Operator = equal ->
213		sat(probformula(Operator,Q,f(not(F))),E)
214		; Operator = greater ->
215		sat(probformula(less,Q,f(not(F))),E)
216		; Operator = less ->
217		sat(probformula(greater,Q,f(not(F))),E)
218		; Operator = strictlygreater ->
219		sat(probformula(strictlyless,Q,f(not(F))),E)
220		; Operator = strictlyless ->
221		sat(probformula(strictlygreater,Q,f(not(F))),E)),
222		P is 1-Q)
223		.
224
225		% Check the type of the formula
226		sat_node(probformula(Operator,P,Ctl_formula),E,Node) :-
227		Ctl_formula = u(F,G) ->
228		prob_calc(u(F,G),E,P_phi,Node),
229		against(P_phi,P,Operator)
230		; Ctl_formula = f(G) ->
231		prob_calc(u(true,G),E,P_phi,Node),
232		against(P_phi,P,Operator)
233		; Ctl_formula = fk(K,G) ->
234		sat_dynamic(probformula(Operator,P,uk(true,K,G)),E,Node)
235		; sat_dynamic(probformula(Operator,P,Ctl_formula),E,Node).
236
237		% Use a different technic depending on the operator
238		% This allow notably the calculation of a probability for the equal operator
239		sat_dynamic(probformula(Operator,P,Ctl_formula),E,Node) :-
240		((Operator = greater ; Operator= strictlygreater) ->
241		ground(P),
242		prob_calc(Ctl_formula,E,P,Operator,Node)
243
244		; Operator = equal ->
245		retractall(prob_current(Node,_)),
246		(ground(P) ->
247		prob_calc(Ctl_formula,E,P,equal,Node)
248
249		; (prob_calc(Ctl_formula,E,1.0,equal,Node) ->
250		P=1.0
251		; prob_current(Node,P)))
252
253		; Operator = less ->
254		ground(P),
255		\+(prob_calc(Ctl_formula,E,P,strictlygreater,Node))
256
257		; Operator = strictlyless ->
258		ground(P),
259		\+(prob_calc(Ctl_formula,E,P,greater,Node))
260		).
261
262		% Compare different formulas using a specific operator
263
264		% For the equal comparison, we compare the results using epsilon precision in case
265		% of a given probability
266		against(P_phi,ReferenceP,equal) :- !,
267		(ground_number(ReferenceP,P) ->
268		P_phi =< P + 0.00000000000000023,
269		P =< P_phi + 0.00000000000000023
270		; ReferenceP = P_phi % ReferenceP is an open variable
271		).
272		against(P_phi,ReferenceP,less) :-
273		ground_number(ReferenceP,P), !,
274		P_phi =< P + 0.00000000000000023. % less
275		against(P_phi,ReferenceP,greater) :-
276		ground_number(ReferenceP,P), !,
277		P_phi >= P - 0.00000000000000023. % greater
278		against(P_phi,ReferenceP,strictlyless) :-
279		ground_number(ReferenceP,P), !,
280		P_phi + 0.00000000000000023 < P . % strictly less
281		against(P_phi,ReferenceP,strictlygreater) :-
282		ground_number(ReferenceP,P), !,
283		P_phi - 0.00000000000000023 > P . % strictlygreater
284		against(P_phi,ReferenceP,not(equal)) :- !,
285		(ground_number(ReferenceP,P) ->
286		(P_phi > P + 0.00000000000000023 ;
287		P =< P_phi + 0.00000000000000023
288		)
289		; % ReferenceP is an open variable
290		add_warning(dtmc_model_checking,'Using symbolic variable for inequality: ',P_phi),
291		(P_phi < 1-0.00000000000000023
292		-> ReferenceP = P_phi + 0.00000000000000023
293		; ReferenceP = P_phi - 0.00000000000000023
294		)
295		).
296		against(P_phi,P,not(less)) :- !, against(P_phi,P,strictlygreater).
297		against(P_phi,P,not(greater)) :- !, against(P_phi,P,strictlyless).
298		against(P_phi,P,not(strictlyless)) :- !, against(P_phi,P,greater).
299		against(P_phi,P,not(strictlygreater)) :- !, against(P_phi,P,less).
300		against(P1,P2,Op) :- add_internal_error('Illegal comparison operator: ',against(P1,P2,Op)),fail.
301
302		negate_operator(not(Op),R) :- !, R=Op.
303		negate_operator(Op,not(Op)).
304
305		:- use_module(probsrc(tools),[safe_number_codes/2]).
306		% TODO: in future fully pre-process AST of formula once before launching model checker
307		ground_number(P,Res) :- number(P),!,Res=P.
308		ground_number(P,Res) :- atom(P), atom_codes(P,Codes), % convert atom from parser into number
309		safe_number_codes(Nr,Codes),!,Res=Nr.
310
311
312		% Next formula
313		:- dynamic prob_current/2.
314
315		prob_calc(x(F),E,P_phi,Operator,Node) :-
316		retractall(prob_current(Node,_)),
317		assert(prob_current(Node,0.0)),
318		(prob_calc_sub(x(F),E,P_phi,Operator,Node) ->
319		true
320		; against(0.0,P_phi,Operator)),!.
321
322		% Until Bounded formula
323		prob_calc(uk(F,K,G),E,P_phi,Operator,Node) :-
324		retractall(prob_current(Node,_)),
325		assert(prob_current(Node,0.0)),
326		(sat(G,E) ->
327		retractall(prob_current(Node,_)),
328		assert(prob_current(Node,1.0)),
329		against(1.0,P_phi,Operator)
330		; sat(F,E) ->
331		(prob_calc_sub(uk(F,K,G),E,P_phi,1.0,Operator,Node) ->
332		true
333		; against(0.0,P_phi,Operator))
334		; against(0.0,P_phi,Operator)),!.
335
336		% Until formula
337		% For this formula we have to calculate the probability for
338		% all states
339		prob_calc(u(F,G),E,P_phi,Node) :-
340		retractall(table_prob0(_,_,Node)),
341		retractall(table_prob1(_,_,Node)),
342		prob_calc_u1(F,G,List_S,P_vect,Node),
343		state(E),
344		(find_prob_u(List_S,P_vect,E,P) -> P_phi=P
345		; P_phi=0.0
346		).
347
348
349		%*******************************************************
350
351		% recursion for the next formula
352		prob_calc_sub(x(F),E,P_phi,Operator,Node) :-
353		transition_probability(E,S,P),
354		sat(F,S),
355		prob_current(Node,Previous_P),
356		Current_prob is Previous_P +P,
357		retract(prob_current(Node,Previous_P)),
358		assert(prob_current(Node,Current_prob)),
359		against(Current_prob,P_phi,Operator).
360
361		% recursion for the bounded until formula
362		prob_calc_sub(uk(F,K_new,G),E,P_phi,P_trace,Operator,Node) :-
363		(sat(G,E) ->
364		prob_current(Node,P),
365		P_new is P+P_trace,
366		retract(prob_current(Node,P)),
367		assert(prob_current(Node,P_new)),
368		against(P_new,P_phi,Operator)
369		; K_new > 0,
370		transition_probability(E,S,P_trans),
371		sat(F,S),
372		K is K_new -1,
373		P_trace_new is P_trace*P_trans,
374		prob_calc_sub(uk(F,K,G),S,P_phi,P_trace_new,Operator,Node)
375		).
376
377		% 1rst precomputation for the until formula
378
379		:- dynamic table_prob0/3.
380
381		prob0(_F,_G,E,Node) :-
382		table_prob0(E,true,Node),!.
383
384		prob0(_F,G,E,Node) :-
385		sat(G,E),
386		!,
387		asserta(table_prob0(E,true,Node)).
388
389		prob0(F,G,E1,Node) :-
390		sat(F,E1),
391		assertz(table_prob0(E1,computing,Node)),
392		trans(E1,E2,_P,_TransId),
393		\+ table_prob0(E2,computing,Node),
394		(prob0(F,G,E2,Node) ->
395		asserta(table_prob0(E1,true,Node)),
396		retract(table_prob0(E1,computing,Node))
397		),!.
398
399		search_prob0(F,G,E,Node):-
400		retractall(table_prob0(_,computing,Node)),
401		prob0(F,G,E,Node).
402
403		% 2nd precomputation for the until formula
404		:- dynamic table_prob1/3.
405
406		prob1(_F,_G,E,Node) :-
407		table_prob1(E,true,Node),!.
408
409		prob1(_F,_G,E,Node) :-
410		\+(table_prob0(E,true,Node)),
411		!,
412		asserta(table_prob1(E,true,Node)).
413
414		prob1(F,G,E1,Node) :-
415		sat(F,E1),
416		\+(sat(G,E1)),
417		assert(table_prob1(E1,computing,Node)),
418		trans(E1,E2,_P,_),
419		\+ table_prob1(E2,computing,Node),
420		(prob1(F,G,E2,Node) ->
421		asserta(table_prob1(E1,true,Node)),
422		retract(table_prob1(E1,computing,Node))
423		),!.
424
425		search_prob1(F,G,E,Node):-
426		retractall(table_prob1(_,computing,Node)),
427		prob1(F,G,E,Node).
428
429		% P_vect is the vector of non-null probabilities
430		% for the until formula
431		prob_calc_u1(F,G,List_S,P_vect,Node) :-
432		findall(S,(state(S),search_prob0(F,G,S,Node)),List_S),
433		findall(E,(state(E),search_prob1(F,G,E,Node)),_List),
434		prob_calc_u2(F,G,List_S,List_S,P_vect,P_vect,Node).
435
436		prob_calc_u2(_F,_G,[],_List_S,[],_P_vect,_Node).
437		prob_calc_u2(F,G,[S|List_S_explored],List_S,[P_phi|P_vect_explored],P_vect,Node) :-
438		(table_prob1(S,true,Node) -> prob_calc_u3(S,List_S,P_phi,P_vect)
439		; P_phi=1.0
440		),
441		prob_calc_u2(F,G,List_S_explored,List_S,P_vect_explored,P_vect,Node).
442
443		prob_calc_u3(_S,[],0.0,[]).
444		prob_calc_u3(S,[E|List_S],P_phi_new,[P|P_vect]) :-
445		prob_calc_u3(S,List_S,P_phi,P_vect),
446		(transition_probability(S,E,P_trans) -> {P_phi_new = P*P_trans + P_phi}
447		; P_phi_new=P_phi
448		).
449
450		find_prob_u([E],[P],E,P).
451		find_prob_u([S1,S2|List_S],[Prob1,Prob2|P_vect],E,P) :-
452		(E=S1,P=Prob1)
453		; find_prob_u([S2|List_S],[Prob2|P_vect],E,P).