dtmc_model_checking

against(P_phi,ReferenceP,equal) :- !,
    (ground_number(ReferenceP,P) ->
    P_phi =< P + 0.00000000000000023,
    P =< P_phi + 0.00000000000000023
    ; ReferenceP = P_phi % ReferenceP is an open variable
).

against(P_phi,ReferenceP,less) :-
    ground_number(ReferenceP,P), !,
    P_phi =< P + 0.00000000000000023. % less

against(P_phi,ReferenceP,greater) :-
    ground_number(ReferenceP,P), !,
    P_phi >= P - 0.00000000000000023. % greater

against(P_phi,ReferenceP,strictlyless) :-
    ground_number(ReferenceP,P), !,
    P_phi + 0.00000000000000023 < P . % strictly less

against(P_phi,ReferenceP,strictlygreater) :-
    ground_number(ReferenceP,P), !,
    P_phi - 0.00000000000000023 > P . % strictlygreater

against(P_phi,ReferenceP,not(equal)) :- !,
    (ground_number(ReferenceP,P) ->
    (P_phi > P + 0.00000000000000023 ;
    P =< P_phi + 0.00000000000000023
    )
    ; % ReferenceP is an open variable
    add_warning(dtmc_model_checking,'Using symbolic variable for inequality: ',P_phi),
    (P_phi < 1-0.00000000000000023
    -> ReferenceP = P_phi + 0.00000000000000023
    ; ReferenceP = P_phi - 0.00000000000000023
)
).

against(P_phi,P,not(less)) :- !, against(P_phi,P,strictlygreater).

against(P_phi,P,not(greater)) :- !, against(P_phi,P,strictlyless).

against(P_phi,P,not(strictlyless)) :- !, against(P_phi,P,greater).

against(P_phi,P,not(strictlygreater)) :- !, against(P_phi,P,less).

against(P1,P2,Op) :- add_internal_error('Illegal comparison operator: ',against(P1,P2,Op)),fail.

find_prob_u([E],[P],E,P).

find_prob_u([S1,S2|List_S],[Prob1,Prob2|P_vect],E,P) :-
    (E=S1,P=Prob1)
    ; find_prob_u([S2|List_S],[Prob2|P_vect],E,P).

ground_number(P,Res) :- number(P),!,Res=P.

ground_number(P,Res) :- atom(P), atom_codes(P,Codes), % convert atom from parser into number
    safe_number_codes(Nr,Codes),!,Res=Nr.

negate_operator(not(Op),R) :- !, R=Op.

negate_operator(Op,not(Op)).

pctl_model_check(Formula,_MaxNodes,Mode,Res) :-
    (temporal_parser(Formula,pctl,PCtlFormula) -> true ; add_error(ctl,'PCTL Parser failed: ',Formula),fail),
    %set_max_nr_of_new_impl_trans_nodes(MaxNodes),
    % TODO: provide better feedback, extract probability for P={p} formulas
    (pre_process_formula(PCtlFormula,ProcessedFormula,Bindings,[]) -> true
    ; add_error(dtmc_model_checking,'Could not pre-process: ',PCtlFormula),
    ProcessedFormula = PCtlFormula
    ),
    (Mode= specific_node(ID) -> Start=ID
    ; Mode = starthere -> current_state_id(Start)
    ; Mode = init -> find_initialised_states(SN), member(Start,SN) % TODO try out all start nodes??
    ; add_error(dtmc_model_checking,'Illegal starting mode (init, starthere supported):',Mode),
    Start = root
    ),
    format('Checking PCTL formula from state ~w : ~w~n AST: ~w~n Open Probability Variables: ~w~n',[Start,Formula,ProcessedFormula,Bindings]),
    reset_tables,
    start_ms_timer(T1),
    (sat(ProcessedFormula,Start)
    -> (Bindings=[] -> Res=true ; Res = solution(Bindings))
    ; Res = false),
    stop_ms_timer_with_msg(T1,'PTCTL model checking').

pre_process_formula(true,true) --> [].

pre_process_formula(false,false) --> [].

pre_process_formula(ap(P),ap(P)) --> [].

pre_process_formula(p(P),p(P)) --> [].

pre_process_formula(probformula(Operator,RefProbAtom,Ctl_formula),
    probformula(Operator,RefProbNr,Ctl_formula)) -->
    ({ground_number(RefProbAtom,RefProbNr)} -> []
    ; {Operator=equal} -> [RefProbAtom/RefProbNr]
    ; {add_error(dtmc_model_checking,'Illegal reference probability for operator (only equal is supported for symbolic values):',RefProbAtom/Operator)},
    [RefProbAtom/RefProbNr]
). % TODO: pre_process_path_formula

pre_process_formula(and(F,G),and(PF,PG)) -->
    pre_process_formula(F,PF),
    pre_process_formula(G,PG).

pre_process_formula(or(F,G),and(PF,PG)) -->
    pre_process_formula(F,PF),
    pre_process_formula(G,PG).

pre_process_formula(implies(F,G),and(PF,PG)) -->
    pre_process_formula(F,PF),
    pre_process_formula(G,PG).

pre_process_formula(not(F),not(PF)) -->
    pre_process_formula(F,PF).

prob0(_F,_G,E,Node) :-
    table_prob0(E,true,Node),!.

prob0(_F,G,E,Node) :-
    sat(G,E),
    !,
    asserta(table_prob0(E,true,Node)).

prob0(F,G,E1,Node) :-
    sat(F,E1),
    assertz(table_prob0(E1,computing,Node)),
    trans(E1,E2,_P),
    \+ table_prob0(E2,computing,Node),
    (prob0(F,G,E2,Node) ->
    asserta(table_prob0(E1,true,Node)),
    retract(table_prob0(E1,computing,Node))
    ),!.

prob1(_F,_G,E,Node) :-
    table_prob1(E,true,Node),!.

prob1(_F,_G,E,Node) :-
    \+(table_prob0(E,true,Node)),
    !,
    asserta(table_prob1(E,true,Node)).

prob1(F,G,E1,Node) :-
    sat(F,E1),
    \+(sat(G,E1)),
    assert(table_prob1(E1,computing,Node)),
    trans(E1,E2,_P),
    \+ table_prob1(E2,computing,Node),
    (prob1(F,G,E2,Node) ->
    asserta(table_prob1(E1,true,Node)),
    retract(table_prob1(E1,computing,Node))
    ),!.

prob_calc(u(F,G),E,P_phi,Node) :-
    retractall(table_prob0(_,_,Node)),
    retractall(table_prob1(_,_,Node)),
    prob_calc_u1(F,G,List_S,P_vect,Node),
    state(E),
    (find_prob_u(List_S,P_vect,E,P) -> P_phi=P
    ; P_phi=0.0
).

prob_calc(x(F),E,P_phi,Operator,Node) :-
    retractall(prob_current(Node,_)),
    assert(prob_current(Node,0.0)),
    (prob_calc_sub(x(F),E,P_phi,Operator,Node) ->
    true
    ; against(0.0,P_phi,Operator)),!.

prob_calc(uk(F,K,G),E,P_phi,Operator,Node) :-
    retractall(prob_current(Node,_)),
    assert(prob_current(Node,0.0)),
    (sat(G,E) ->
    retractall(prob_current(Node,_)),
    assert(prob_current(Node,1.0)),
    against(1.0,P_phi,Operator)
    ; sat(F,E) ->
    (prob_calc_sub(uk(F,K,G),E,P_phi,1.0,Operator,Node) ->
    true
    ; against(0.0,P_phi,Operator))
    ; against(0.0,P_phi,Operator)),!.

prob_calc_sub(x(F),E,P_phi,Operator,Node) :-
    transition_probability(E,S,P),
    sat(F,S),
    prob_current(Node,Previous_P),
    Current_prob is Previous_P +P,
    retract(prob_current(Node,Previous_P)),
    assert(prob_current(Node,Current_prob)),
    against(Current_prob,P_phi,Operator).

prob_calc_sub(uk(F,K_new,G),E,P_phi,P_trace,Operator,Node) :-
    (sat(G,E) ->
    prob_current(Node,P),
    P_new is P+P_trace,
    retract(prob_current(Node,P)),
    assert(prob_current(Node,P_new)),
    against(P_new,P_phi,Operator)
    ; K_new > 0,
    transition_probability(E,S,P_trans),
    sat(F,S),
    K is K_new -1,
    P_trace_new is P_trace*P_trans,
    prob_calc_sub(uk(F,K,G),S,P_phi,P_trace_new,Operator,Node)
).

prob_calc_u1(F,G,List_S,P_vect,Node) :-
    findall(S,(state(S),search_prob0(F,G,S,Node)),List_S),
    findall(E,(state(E),search_prob1(F,G,E,Node)),_List),
    prob_calc_u2(F,G,List_S,List_S,P_vect,P_vect,Node).

prob_calc_u2(_F,_G,[],_List_S,[],_P_vect,_Node).

prob_calc_u2(F,G,[S|List_S_explored],List_S,[P_phi|P_vect_explored],P_vect,Node) :-
    (table_prob1(S,true,Node) -> prob_calc_u3(S,List_S,P_phi,P_vect)
    ; P_phi=1.0
    ),
    prob_calc_u2(F,G,List_S_explored,List_S,P_vect_explored,P_vect,Node).

prob_calc_u3(_S,[],0.0,[]).

prob_calc_u3(S,[E|List_S],P_phi_new,[P|P_vect]) :-
    prob_calc_u3(S,List_S,P_phi,P_vect),
    (transition_probability(S,E,P_trans) -> {P_phi_new = P*P_trans + P_phi}
    ; P_phi_new=P_phi
).

reset_tables :- retractall(table_prob0(_,_,_)),retractall(node(_)), retractall(prob_current(_,_)).

sat(Formula) :- start(E), sat(Formula,E).

sat(true,_E).

sat(false,_E):-fail.

sat(ap(Property),E) :- check_ap(Property,E).

sat(p(Property),E) :- check_ap(Property,E).

sat(and(F,G),E) :- sat(F,E), sat(G,E).

sat(or(F,_G),E) :- sat(F,E).

sat(or(_F,G),E) :- sat(G,E).

sat(implies(F,G),E) :- sat(or(not(F),G),E).

sat(equivalence(F,G),E) :- sat(and(implies(F,G),implies(G,F)),E).

sat(not(F),E) :- probformula(_,_,_)\=F, sat_not(F,E).

sat(probformula(Operator,P,Ctl_formula),E) :-
    Ctl_formula = gk(_,_) ->
    sat_gk(probformula(Operator,P,Ctl_formula),E)
    ; Ctl_formula = g(_) ->
    sat_g(probformula(Operator,P,Ctl_formula),E)
    ; (node(Node) ->
    New_Node is Node +1,
    retract(node(Node)),
    assert(node(New_Node)),
    sat_node(probformula(Operator,P,Ctl_formula),E,New_Node)
    ; assert(node(0)),
    sat_node(probformula(Operator,P,Ctl_formula),E,0),
    retractall(node(_)) /*reinitialize nodes*/
).

sat(not(probformula(Operator,P,Ctl_formula)),E) :-
    negate_operator(Operator,NotOp),
    sat(probformula(NotOp,P,Ctl_formula),E).

sat(probformula(Operator,P,gk(K,F)),E) :-
    ground(P) ->
    Q is 1-P,
    (Operator = equal ->
    sat(probformula(Operator,Q,fk(K,not(F))),E)
    ; sat(not(probformula(Operator,Q,fk(K,not(F)))),E))
    ; ((Operator = equal ->
    sat(probformula(Operator,Q,fk(K,not(F))),E)
    ; sat(not(probformula(Operator,Q,fk(K,not(F)))),E)),
    P is 1-Q)
.

sat_dynamic(probformula(Operator,P,Ctl_formula),E,Node) :-
    ((Operator = greater ; Operator= strictlygreater) ->
    ground(P),
    prob_calc(Ctl_formula,E,P,Operator,Node)
   
    ; Operator = equal ->
    retractall(prob_current(Node,_)),
    (ground(P) ->
    prob_calc(Ctl_formula,E,P,equal,Node)
   
    ; (prob_calc(Ctl_formula,E,1.0,equal,Node) ->
    P=1.0
    ; prob_current(Node,P)))
   
    ; Operator = less ->
    ground(P),
    \+(prob_calc(Ctl_formula,E,P,strictlygreater,Node))
   
    ; Operator = strictlyless ->
    ground(P),
    \+(prob_calc(Ctl_formula,E,P,greater,Node))
).

sat_g(probformula(Operator,P,g(F)),E) :-
    ground(P) ->
    Q is 1-P,
    (Operator = equal ->
    sat(probformula(Operator,Q,f(not(F))),E)
    ; Operator = greater ->
    sat(probformula(less,Q,f(not(F))),E)
    ; Operator = less ->
    sat(probformula(greater,Q,f(not(F))),E)
    ; Operator = strictlygreater ->
    sat(probformula(strictlyless,Q,f(not(F))),E)
    ; Operator = strictlyless ->
    sat(probformula(strictlygreater,Q,f(not(F))),E))
    ; ((Operator = equal ->
    sat(probformula(Operator,Q,f(not(F))),E)
    ; Operator = greater ->
    sat(probformula(less,Q,f(not(F))),E)
    ; Operator = less ->
    sat(probformula(greater,Q,f(not(F))),E)
    ; Operator = strictlygreater ->
    sat(probformula(strictlyless,Q,f(not(F))),E)
    ; Operator = strictlyless ->
    sat(probformula(strictlygreater,Q,f(not(F))),E)),
    P is 1-Q)
.

sat_gk(probformula(Operator,P,gk(K,F)),E) :-
    ground(P) ->
    Q is 1-P,
    (Operator = equal ->
    sat(probformula(Operator,Q,fk(K,not(F))),E)
    ; Operator = greater ->
    sat(probformula(less,Q,fk(K,not(F))),E)
    ; Operator = less ->
    sat(probformula(greater,Q,fk(K,not(F))),E)
    ; Operator = strictlygreater ->
    sat(probformula(strictlyless,Q,fk(K,not(F))),E)
    ; Operator = strictlyless ->
    sat(probformula(strictlygreater,Q,fk(K,not(F))),E))
    ; ((Operator = equal ->
    sat(probformula(Operator,Q,fk(K,not(F))),E)
    ; Operator = greater ->
    sat(probformula(less,Q,fk(K,not(F))),E)
    ; Operator = less ->
    sat(probformula(greater,Q,fk(K,not(F))),E)
    ; Operator = strictlygreater ->
    sat(probformula(strictlyless,Q,fk(K,not(F))),E)
    ; Operator = strictlyless ->
    sat(probformula(strictlygreater,Q,fk(K,not(F))),E)),
    P is 1-Q)
.

sat_node(probformula(Operator,P,Ctl_formula),E,Node) :-
    Ctl_formula = u(F,G) ->
    prob_calc(u(F,G),E,P_phi,Node),
    against(P_phi,P,Operator)
    ; Ctl_formula = f(G) ->
    prob_calc(u(true,G),E,P_phi,Node),
    against(P_phi,P,Operator)
    ; Ctl_formula = fk(K,G) ->
    sat_dynamic(probformula(Operator,P,uk(true,K,G)),E,Node)
    ; sat_dynamic(probformula(Operator,P,Ctl_formula),E,Node).

sat_not(false,_E).

sat_not(true,_E) :- fail.

sat_not(p(Property),E) :- \+(check_ap(Property,E)).

sat_not(ap(Property),E) :- \+(check_ap(Property,E)).

sat_not(and(F,_G),E) :- sat_not(F,E).

sat_not(and(_F,G),E) :- sat_not(G,E).

sat_not(or(F,G),E) :- sat_not(F,E),sat_not(G,E).

sat_not(implies(F,G),E) :- sat(not(or(not(F),G)),E).

sat_not(equivalence(F,G),E) :- sat_not(and(implies(F,G),implies(G,F)),E).

sat_not(not(F),E) :- sat(F,E).

search_prob0(F,G,E,Node):-
    retractall(table_prob0(_,computing,Node)),
    prob0(F,G,E,Node).

search_prob1(F,G,E,Node):-
    retractall(table_prob1(_,computing,Node)),
    prob1(F,G,E,Node).

start(ID) :- find_initialised_states(I), member(ID,I).

state(X) :- visited_expression_id(X).

transition_probability(From,To,Prob) :- trans(From,To,_),
    visited_expression(From,FS),
    visited_expression(To,TS),
    xtl_transition(FS,_,TS,Infos),
    (member(probability/X,Infos) -> Prob=X ; write(no_probability(From,To,Infos)),nl,fail).