1		% (c) 2025-2025 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		sat(Formula) :- start(E), sat(Formula,E).
107
108
109		% Classic cases
110
111		sat_not(false,_E).
112		sat_not(true,_E) :- fail.
113		sat_not(p(Property),E) :- \+(check_ap(Property,E)).
114		sat_not(ap(Property),E) :- \+(check_ap(Property,E)).
115		sat_not(and(F,_G),E) :- sat_not(F,E).
116		sat_not(and(_F,G),E) :- sat_not(G,E).
117		sat_not(or(F,G),E) :- sat_not(F,E),sat_not(G,E).
118		sat_not(implies(F,G),E) :- sat(not(or(not(F),G)),E).
119		sat_not(equivalence(F,G),E) :- sat_not(and(implies(F,G),implies(G,F)),E).
120		sat_not(not(F),E) :- sat(F,E).
121
122		%sat(Formula,State) :- write(sat(Formula,State)),nl,fail.
123		sat(true,_E).
124		sat(false,_E):-fail.
125		sat(ap(Property),E) :- check_ap(Property,E).
126		sat(p(Property),E) :- check_ap(Property,E).
127		sat(and(F,G),E) :- sat(F,E), sat(G,E).
128		sat(or(F,_G),E) :- sat(F,E).
129		sat(or(_F,G),E) :- sat(G,E).
130		sat(implies(F,G),E) :- sat(or(not(F),G),E).
131		sat(equivalence(F,G),E) :- sat(and(implies(F,G),implies(G,F)),E).
132		sat(not(F),E) :- probformula(_,_,_)\=F, sat_not(F,E).
133
134		% Probabilistic-formula cases. Operator is =, < >,<= or >=. P is a number between 0 and 1 (or a Variable),
135		% E is a state.
136		% Check if the formula is nested
137
138		sat(probformula(Operator,P,Ctl_formula),E) :-
139		Ctl_formula = gk(_,_) ->
140		sat_gk(probformula(Operator,P,Ctl_formula),E)
141		; Ctl_formula = g(_) ->
142		sat_g(probformula(Operator,P,Ctl_formula),E)
143		; (node(Node) ->
144		New_Node is Node +1,
145		retract(node(Node)),
146		assert(node(New_Node)),
147		sat_node(probformula(Operator,P,Ctl_formula),E,New_Node)
148		; assert(node(0)),
149		sat_node(probformula(Operator,P,Ctl_formula),E,0),
150		retractall(node(_)) /*reinitialize nodes*/
151		).
152
153		sat(not(probformula(Operator,P,Ctl_formula)),E) :-
154		negate_operator(Operator,NotOp),
155		sat(probformula(NotOp,P,Ctl_formula),E).
156
157		% Always bounded formula, we use the dual probabilistic event of fk(K,not(F))
158		sat(probformula(Operator,P,gk(K,F)),E) :-
159		ground(P) ->
160		Q is 1-P,
161		(Operator = equal ->
162		sat(probformula(Operator,Q,fk(K,not(F))),E)
163		; sat(not(probformula(Operator,Q,fk(K,not(F)))),E))
164		; ((Operator = equal ->
165		sat(probformula(Operator,Q,fk(K,not(F))),E)
166		; sat(not(probformula(Operator,Q,fk(K,not(F)))),E)),
167		P is 1-Q)
168		.
169
170		% Always bounded formula, we use the dual probabilistic event of fk(K,not(F))
171		sat_gk(probformula(Operator,P,gk(K,F)),E) :-
172		ground(P) ->
173		Q is 1-P,
174		(Operator = equal ->
175		sat(probformula(Operator,Q,fk(K,not(F))),E)
176		; Operator = greater ->
177		sat(probformula(less,Q,fk(K,not(F))),E)
178		; Operator = less ->
179		sat(probformula(greater,Q,fk(K,not(F))),E)
180		; Operator = strictlygreater ->
181		sat(probformula(strictlyless,Q,fk(K,not(F))),E)
182		; Operator = strictlyless ->
183		sat(probformula(strictlygreater,Q,fk(K,not(F))),E))
184		; ((Operator = equal ->
185		sat(probformula(Operator,Q,fk(K,not(F))),E)
186		; Operator = greater ->
187		sat(probformula(less,Q,fk(K,not(F))),E)
188		; Operator = less ->
189		sat(probformula(greater,Q,fk(K,not(F))),E)
190		; Operator = strictlygreater ->
191		sat(probformula(strictlyless,Q,fk(K,not(F))),E)
192		; Operator = strictlyless ->
193		sat(probformula(strictlygreater,Q,fk(K,not(F))),E)),
194		P is 1-Q)
195		.
196
197		% Always formula
198		sat_g(probformula(Operator,P,g(F)),E) :-
199		ground(P) ->
200		Q is 1-P,
201		(Operator = equal ->
202		sat(probformula(Operator,Q,f(not(F))),E)
203		; Operator = greater ->
204		sat(probformula(less,Q,f(not(F))),E)
205		; Operator = less ->
206		sat(probformula(greater,Q,f(not(F))),E)
207		; Operator = strictlygreater ->
208		sat(probformula(strictlyless,Q,f(not(F))),E)
209		; Operator = strictlyless ->
210		sat(probformula(strictlygreater,Q,f(not(F))),E))
211		; ((Operator = equal ->
212		sat(probformula(Operator,Q,f(not(F))),E)
213		; Operator = greater ->
214		sat(probformula(less,Q,f(not(F))),E)
215		; Operator = less ->
216		sat(probformula(greater,Q,f(not(F))),E)
217		; Operator = strictlygreater ->
218		sat(probformula(strictlyless,Q,f(not(F))),E)
219		; Operator = strictlyless ->
220		sat(probformula(strictlygreater,Q,f(not(F))),E)),
221		P is 1-Q)
222		.
223
224		% Check the type of the formula
225		sat_node(probformula(Operator,P,Ctl_formula),E,Node) :-
226		Ctl_formula = u(F,G) ->
227		prob_calc(u(F,G),E,P_phi,Node),
228		against(P_phi,P,Operator)
229		; Ctl_formula = f(G) ->
230		prob_calc(u(true,G),E,P_phi,Node),
231		against(P_phi,P,Operator)
232		; Ctl_formula = fk(K,G) ->
233		sat_dynamic(probformula(Operator,P,uk(true,K,G)),E,Node)
234		; sat_dynamic(probformula(Operator,P,Ctl_formula),E,Node).
235
236		% Use a different technic depending on the operator
237		% This allow notably the calculation of a probability for the equal operator
238		sat_dynamic(probformula(Operator,P,Ctl_formula),E,Node) :-
239		((Operator = greater ; Operator= strictlygreater) ->
240		ground(P),
241		prob_calc(Ctl_formula,E,P,Operator,Node)
242
243		; Operator = equal ->
244		retractall(prob_current(Node,_)),
245		(ground(P) ->
246		prob_calc(Ctl_formula,E,P,equal,Node)
247
248		; (prob_calc(Ctl_formula,E,1.0,equal,Node) ->
249		P=1.0
250		; prob_current(Node,P)))
251
252		; Operator = less ->
253		ground(P),
254		\+(prob_calc(Ctl_formula,E,P,strictlygreater,Node))
255
256		; Operator = strictlyless ->
257		ground(P),
258		\+(prob_calc(Ctl_formula,E,P,greater,Node))
259		).
260
261		% Compare different formulas using a specific operator
262
263		% For the equal comparison, we compare the results using epsilon precision in case
264		% of a given probability
265		against(P_phi,ReferenceP,equal) :- !,
266		(ground_number(ReferenceP,P) ->
267		P_phi =< P + 0.00000000000000023,
268		P =< P_phi + 0.00000000000000023
269		; ReferenceP = P_phi % ReferenceP is an open variable
270		).
271		against(P_phi,ReferenceP,less) :-
272		ground_number(ReferenceP,P), !,
273		P_phi =< P + 0.00000000000000023. % less
274		against(P_phi,ReferenceP,greater) :-
275		ground_number(ReferenceP,P), !,
276		P_phi >= P - 0.00000000000000023. % greater
277		against(P_phi,ReferenceP,strictlyless) :-
278		ground_number(ReferenceP,P), !,
279		P_phi + 0.00000000000000023 < P . % strictly less
280		against(P_phi,ReferenceP,strictlygreater) :-
281		ground_number(ReferenceP,P), !,
282		P_phi - 0.00000000000000023 > P . % strictlygreater
283		against(P_phi,ReferenceP,not(equal)) :- !,
284		(ground_number(ReferenceP,P) ->
285		(P_phi > P + 0.00000000000000023 ;
286		P =< P_phi + 0.00000000000000023
287		)
288		; % ReferenceP is an open variable
289		add_warning(dtmc_model_checking,'Using symbolic variable for inequality: ',P_phi),
290		(P_phi < 1-0.00000000000000023
291		-> ReferenceP = P_phi + 0.00000000000000023
292		; ReferenceP = P_phi - 0.00000000000000023
293		)
294		).
295		against(P_phi,P,not(less)) :- !, against(P_phi,P,strictlygreater).
296		against(P_phi,P,not(greater)) :- !, against(P_phi,P,strictlyless).
297		against(P_phi,P,not(strictlyless)) :- !, against(P_phi,P,greater).
298		against(P_phi,P,not(strictlygreater)) :- !, against(P_phi,P,less).
299		against(P1,P2,Op) :- add_internal_error('Illegal comparison operator: ',against(P1,P2,Op)),fail.
300
301		negate_operator(not(Op),R) :- !, R=Op.
302		negate_operator(Op,not(Op)).
303
304		:- use_module(probsrc(tools),[safe_number_codes/2]).
305		% TODO: in future fully pre-process AST of formula once before launching model checker
306		ground_number(P,Res) :- number(P),!,Res=P.
307		ground_number(P,Res) :- atom(P), atom_codes(P,Codes), % convert atom from parser into number
308		safe_number_codes(Nr,Codes),!,Res=Nr.
309
310
311		% Next formula
312		:- dynamic prob_current/2.
313
314		prob_calc(x(F),E,P_phi,Operator,Node) :-
315		retractall(prob_current(Node,_)),
316		assert(prob_current(Node,0.0)),
317		(prob_calc_sub(x(F),E,P_phi,Operator,Node) ->
318		true
319		; against(0.0,P_phi,Operator)),!.
320
321		% Until Bounded formula
322		prob_calc(uk(F,K,G),E,P_phi,Operator,Node) :-
323		retractall(prob_current(Node,_)),
324		assert(prob_current(Node,0.0)),
325		(sat(G,E) ->
326		retractall(prob_current(Node,_)),
327		assert(prob_current(Node,1.0)),
328		against(1.0,P_phi,Operator)
329		; sat(F,E) ->
330		(prob_calc_sub(uk(F,K,G),E,P_phi,1.0,Operator,Node) ->
331		true
332		; against(0.0,P_phi,Operator))
333		; against(0.0,P_phi,Operator)),!.
334
335		% Until formula
336		% For this formula we have to calculate the probability for
337		% all states
338		prob_calc(u(F,G),E,P_phi,Node) :-
339		retractall(table_prob0(_,_,Node)),
340		retractall(table_prob1(_,_,Node)),
341		prob_calc_u1(F,G,List_S,P_vect,Node),
342		state(E),
343		(find_prob_u(List_S,P_vect,E,P) -> P_phi=P
344		; P_phi=0.0
345		).
346
347
348		%*******************************************************
349
350		% recursion for the next formula
351		prob_calc_sub(x(F),E,P_phi,Operator,Node) :-
352		transition_probability(E,S,P),
353		sat(F,S),
354		prob_current(Node,Previous_P),
355		Current_prob is Previous_P +P,
356		retract(prob_current(Node,Previous_P)),
357		assert(prob_current(Node,Current_prob)),
358		against(Current_prob,P_phi,Operator).
359
360		% recursion for the bounded until formula
361		prob_calc_sub(uk(F,K_new,G),E,P_phi,P_trace,Operator,Node) :-
362		(sat(G,E) ->
363		prob_current(Node,P),
364		P_new is P+P_trace,
365		retract(prob_current(Node,P)),
366		assert(prob_current(Node,P_new)),
367		against(P_new,P_phi,Operator)
368		; K_new > 0,
369		transition_probability(E,S,P_trans),
370		sat(F,S),
371		K is K_new -1,
372		P_trace_new is P_trace*P_trans,
373		prob_calc_sub(uk(F,K,G),S,P_phi,P_trace_new,Operator,Node)
374		).
375
376		% 1rst precomputation for the until formula
377
378		:- dynamic table_prob0/3.
379
380		prob0(_F,_G,E,Node) :-
381		table_prob0(E,true,Node),!.
382
383		prob0(_F,G,E,Node) :-
384		sat(G,E),
385		!,
386		asserta(table_prob0(E,true,Node)).
387
388		prob0(F,G,E1,Node) :-
389		sat(F,E1),
390		assertz(table_prob0(E1,computing,Node)),
391		trans(E1,E2,_P,_TransId),
392		\+ table_prob0(E2,computing,Node),
393		(prob0(F,G,E2,Node) ->
394		asserta(table_prob0(E1,true,Node)),
395		retract(table_prob0(E1,computing,Node))
396		),!.
397
398		search_prob0(F,G,E,Node):-
399		retractall(table_prob0(_,computing,Node)),
400		prob0(F,G,E,Node).
401
402		% 2nd precomputation for the until formula
403		:- dynamic table_prob1/3.
404
405		prob1(_F,_G,E,Node) :-
406		table_prob1(E,true,Node),!.
407
408		prob1(_F,_G,E,Node) :-
409		\+(table_prob0(E,true,Node)),
410		!,
411		asserta(table_prob1(E,true,Node)).
412
413		prob1(F,G,E1,Node) :-
414		sat(F,E1),
415		\+(sat(G,E1)),
416		assert(table_prob1(E1,computing,Node)),
417		trans(E1,E2,_P,_),
418		\+ table_prob1(E2,computing,Node),
419		(prob1(F,G,E2,Node) ->
420		asserta(table_prob1(E1,true,Node)),
421		retract(table_prob1(E1,computing,Node))
422		),!.
423
424		search_prob1(F,G,E,Node):-
425		retractall(table_prob1(_,computing,Node)),
426		prob1(F,G,E,Node).
427
428		% P_vect is the vector of non-null probabilities
429		% for the until formula
430		prob_calc_u1(F,G,List_S,P_vect,Node) :-
431		findall(S,(state(S),search_prob0(F,G,S,Node)),List_S),
432		findall(E,(state(E),search_prob1(F,G,E,Node)),_List),
433		prob_calc_u2(F,G,List_S,List_S,P_vect,P_vect,Node).
434
435		prob_calc_u2(_F,_G,[],_List_S,[],_P_vect,_Node).
436		prob_calc_u2(F,G,[S|List_S_explored],List_S,[P_phi|P_vect_explored],P_vect,Node) :-
437		(table_prob1(S,true,Node) -> prob_calc_u3(S,List_S,P_phi,P_vect)
438		; P_phi=1.0
439		),
440		prob_calc_u2(F,G,List_S_explored,List_S,P_vect_explored,P_vect,Node).
441
442		prob_calc_u3(_S,[],0.0,[]).
443		prob_calc_u3(S,[E|List_S],P_phi_new,[P|P_vect]) :-
444		prob_calc_u3(S,List_S,P_phi,P_vect),
445		(transition_probability(S,E,P_trans) -> {P_phi_new = P*P_trans + P_phi}
446		; P_phi_new=P_phi
447		).
448
449		find_prob_u([E],[P],E,P).
450		find_prob_u([S1,S2|List_S],[Prob1,Prob2|P_vect],E,P) :-
451		(E=S1,P=Prob1)
452		; find_prob_u([S2|List_S],[Prob2|P_vect],E,P).