:- module(qdecproc,[qdecproc/2,qproject/4]).

/* SP
:- use_module(library(terms)).
*/

/** <module> qdecproc

qdecproc/2 implements a decision procedure for first order formulas over
<Q,{0,1,+},{=, >}>, where Q denotes the set of rational numbers.
A recent (>= 5.6.x) version of SWI-Prolog is needed.

Syntax:

    * a _|linear expression|_ is any clp(Q) linear expression.

    * an _|atomic constraint|_ is a term =|ac(LinExp1 OpRel LinExp2)|=
    where =LinExp1= and =LinExp2= are linear expressions
    and =OpRel= is a relational operator from {=, >, >=, < , =<, =\=}.
    For instance, =|ac(0<1)|= denotes _true_ and =|ac(0>1)|= _false_. 
    
    * _formulas_ are defined inductively:
        * an atomic constraint is a formula;
	* if =F= is a formula, its negation =|-F|= is a formula;
	* if =F1= and =F2= are formulas, then
	the conjunction =|(F1,F2)|=, 
	the disjunction =|(F1;F2)|=, and
	the implication =|(F1 -> F2)|= are formulas;
	* if =F= is a formula and =X= is a Prolog variable, 
	=|forall(X,F)|= and =|exists(X,F)|= are formulas.


Use:

    * the first argument of qdecproc/2 is a formula F.

    * the second argument, computed by the procedure,
    is a formula G in disjunctive normal form without quantifier
    which is equivalent to F. In particular, if F is a _closed_
    formula, i.e. a formula without free variable,
    G is either =|ac(0<1)|= or =|ac(0>1)|=.

At the end of this file, some examples are given in
the SWI-Prolog Unit Tests format.

@author Fred Mesnard 
@version 1.0
@license GPL
@tbd 
    * remove the cuts in disjunctive_normal_form2/2
*/

:- use_module(library(clpq)).
:- use_module(library(lists)).

%%	qdecproc(+F:formula,-SimplifiedF:formula) is det.
%
%	True if SimplifiedF is a quantifier-free dnf equivalent to the
%	formula F.
%	
qdecproc(F0,F) :-
	prenex_dnf(F0,F1),
	quantif_elim(F1,F).
	
%%	quantif_elim(+F:prenex_dnf,-SimplifiedF:dnf) is det.
%
%	True if SimplifiedF is a dnf equivalent to the prenex_dnf F but
%	without quantifiers.
%	
/* the initial version was too slow: we have to factorize quantifiers
quantif_elim(exists(X,F0),F) :-
	quantif_elim(F0,F1),
	elim_exs(F1,[X],F).
quantif_elim(forall(X,F0),F) :-
	quantif_elim(F0,F1),
	prenex_dnf(-F1,F2), 
	elim_exs(F2,[X],F3), 
	prenex_dnf(-F3,F).
*/
quantif_elim(exists(X,F0),F) :-
	fact_quantif_ex(F0,[X],Xs,Fe),
	quantif_elim(Fe,F1),
	elim_exs(F1,Xs,F).
quantif_elim(forall(X,F0),F) :-
	fact_quantif_fa(F0,[X],Xs,Fa),
	quantif_elim(Fa,F1),
	prenex_dnf(-F1,F2), 
	elim_exs(F2,Xs,F3), 
	prenex_dnf(-F3,F).
quantif_elim((C;Cs),(C;Cs)). quantif_elim((A,As),(A,As)).
quantif_elim(ac(A),ac(A)).

fact_quantif_fa(ac(A),Xs,Xs,ac(A)).
fact_quantif_fa((A,As),Xs,Xs,(A,As)).
fact_quantif_fa((C;Cs),Xs,Xs,(C;Cs)).
fact_quantif_fa(exists(X,F),Xs,Xs,exists(X,F)).
fact_quantif_fa(forall(X,F0),Xs,Ys,F) :-
	fact_quantif_fa(F0,[X|Xs],Ys,F).

fact_quantif_ex(ac(A),Xs,Xs,ac(A)).
fact_quantif_ex((A,As),Xs,Xs,(A,As)).
fact_quantif_ex((C;Cs),Xs,Xs,(C;Cs)).
fact_quantif_ex(forall(X,F),Xs,Xs,forall(X,F)).
fact_quantif_ex(exists(X,F0),Xs,Ys,F) :-
	fact_quantif_ex(F0,[X|Xs],Ys,F).

%%	elim_exs(+F:dnf,+Xs:list,-G:dnf) is det.
%
%	True if G is a dnf equivalent to the dnf F where
%	existentially quantified variables belonging to Xs have been
%	eliminated from F.

elim_exs(ac(A),Xs,F) :-
	elim_exs((ac(A),ac(0<1)),Xs,F). 
elim_exs((A,As),Xs,F) :-
	round_list((A,As),F1),
	term_variables(F1,Vars),
	remove(Vars,Xs,Vars1),
	qproject(Vars1,F1,Vars1,F2),
	list_round(F2,F).
elim_exs((C; Cs),Xs,F) :-
	elim_exs(C,Xs,F1),
	elim_exs(Cs,Xs,F2),
	simplify_dnf((F1;F2),F).

round_list(ac(A),[A]).
round_list((ac(A),B), [A|Cs]) :- round_list(B,Cs).

list_round([],ac(0<1)).
list_round([A|As],Rs) :- list_round(As,A,Rs).

list_round([],A,ac(A)).
list_round([B|Cs], A, (ac(A),Ds)):-list_round(Cs,B,Ds).

remove([],_,[]).
remove([X|Xs],Ys,Zs) :- 
	( notin(Ys,X) ->
	     Zs=[X|Ts],
	     remove(Xs,Ys,Ts)
	;
	     remove(Xs,Ys,Zs)).

notin([],_).
notin([Y|Ys],X) :-
	X\==Y,
	notin(Ys,X).

%%%%
simplify_dnf(ac(A),ac(A)).
simplify_dnf((A,As),Cs) :-
	simplify_constraints((A,As),Cs).
simplify_dnf((D0;D1), D) :-
	simplify_constraints(D0,D2),
	simplify_dnf(D1,D3),
	or_combine(D2,D3,D).

simplify_constraints(ac(A),Cs) :-
	simplify_constraints((ac(A),ac(0<1)),Cs).
simplify_constraints((A,As),Cs) :-
	round_list((A,As),Ds),
	term_variables(Ds,Vars),
	qproject(Vars,Ds,Vars,Es),
	list_round(Es,Cs).

or_combine(D1,D2,D3) :-
	(ground(D1) ->
	    (clpq_call(D1) ->
		D3=ac(0<1)
	    ;
		(ground(D2) ->
		    (clpq_call(D2) ->
			D3=ac(0<1)
		    ;
			D3=ac(0>1))
		;
		    D3=D2))
	;
	    (ground(D2) ->
		(clpq_call(D2) ->
		    D3=ac(0<1)
		;
		    D3=D1)
	    ;
		D3=(D1;D2))).
		
clpq_call((ac(D),Ds)) :- {D},clpq_call(Ds).
clpq_call(ac(D)) :- {D}.
clpq_call((A;B)) :- 
	(clpq_call(A) ->
	    true
	;   
	    clpq_call(B)).
	

%%	prenex_dnf(+F:formula,-G:prenex_dnf) is det.
%
%	True if G is a prenex dnf equivalent to the formula F.
%
% adapted from the Prolog programs of Appendix B of:
% SIMPLY LOGICAL: Intelligent reasoning by example,
% Peter A. Flach, Ed. John Wiley & Sons

prenex_dnf(Formula,F5) :-
	rewrite_implication_dif(Formula,F1),
	negations_inside(F1,F2),
	prenex_normal_form(F2,F3),
	disjunctive_normal_form(F3,F5).

% rewrite_implication_dif(F1,F2) 
%    F2 is a predicate logic formula without implication
%    and without =\=, logically equivalent to F1
rewrite_implication_dif(ac(A),F) :- rewrite_implication_dif2(A,F).
rewrite_implication_dif((A -> B),         (-C ; D)) :-	
	rewrite_implication_dif(A,C),
	rewrite_implication_dif(B,D).
rewrite_implication_dif((A,B),            (C,D)) :-
	rewrite_implication_dif(A,C),	
	rewrite_implication_dif(B,D).
rewrite_implication_dif((A ; B),          (C ; D)) :-
	rewrite_implication_dif(A,C),
	rewrite_implication_dif(B,D).
rewrite_implication_dif(-A,               -C) :-
	rewrite_implication_dif(A,C).
rewrite_implication_dif(forall(X,B),      forall(X,D)) :-
	rewrite_implication_dif(B,D).
rewrite_implication_dif(exists(X,B),      exists(X,D)) :-
	rewrite_implication_dif(B,D).

rewrite_implication_dif2(X = Y,        ac(X = Y)).
rewrite_implication_dif2(X =\= Y,      (ac(X < Y) ; ac(X > Y))).
rewrite_implication_dif2(X < Y,        ac(X < Y)).
rewrite_implication_dif2(X =< Y,       ac(X =< Y)).
rewrite_implication_dif2(X > Y,        ac(X > Y)).
rewrite_implication_dif2(X >= Y,       ac(X >= Y)).

% negations_inside(F1,F2) 
%  F2 is a predicate logic formula without negation
%  logically equivalent to F1
negations_inside(ac(A),ac(A)).
negations_inside(- F, G) :- negations_inside2(F,G).
negations_inside((A , B), (C , D)) :-
	negations_inside(A,C),
	negations_inside(B,D).
negations_inside((A ; B), (C ; D)) :-
	negations_inside(A,C),
	negations_inside(B,D).
negations_inside(exists(X,A),exists(X,B)) :-
	negations_inside(A,B).
negations_inside(forall(X,A),forall(X,B)) :-
	negations_inside(A,B).

negations_inside2(ac(A),B) :- 
	neg_atomic_constraint(A,B).
negations_inside2((A,B), (C;D)) :- 
	negations_inside(-A,C),
	negations_inside(-B,D).
negations_inside2((A ; B), (C , D)) :- 
	negations_inside(-A,C),
	negations_inside(-B,D).
negations_inside2(- A,B) :-
	negations_inside(A,B).
negations_inside2(exists(X,A),forall(X,B)) :- 
	negations_inside(-A,B).
negations_inside2(forall(X,A),exists(X,B)) :- 
	negations_inside(-A,B).

neg_atomic_constraint( X = Y, (ac(X < Y) ; ac(X > Y))).
neg_atomic_constraint( X =\= Y, ac(X = Y)).
neg_atomic_constraint( X < Y, ac(X >= Y)).
neg_atomic_constraint( X =< Y, ac(X > Y)).
neg_atomic_constraint( X > Y, ac(X =< Y)).
neg_atomic_constraint( X >= Y, ac(X < Y)).

% prenex_normal_form(F1,F2) 
%  F2 is a predicate logic formula with all 
%  quantifiers in front, logically equivalent to F1

prenex_normal_form(F,PNF) :- 
	rename(F,[],G),
	pnf(G,PNF,B,B).

pnf(forall(X,F),forall(X,Quants),V,Body) :-
	pnf(F,Quants,V,Body).
pnf(exists(X,F),exists(X,Quants),V,Body) :-
	pnf(F,Quants,V,Body).
pnf(ac(A),V,V,ac(A)).
pnf((A , B),Quants,V,(BodyA , BodyB)) :-
	pnf(A,Quants,QB,BodyA),
	pnf(B,QB,V,BodyB).
pnf((A ; B),Quants,V,(BodyA ; BodyB)) :-
	pnf(A,Quants,QB,BodyA),
	pnf(B,QB,V,BodyB).

rename(ac(A),_Vars,ac(A)).
rename((A0,B0),Vars,(A,B)) :-
	rename(A0,Vars,A),
	rename(B0,Vars,B).
rename((A0;B0),Vars,(A;B)) :-
	rename(A0,Vars,A),
	rename(B0,Vars,B).
rename(forall(X,F0),Vars,forall(Y,F)) :-
	(notin(Vars,X) ->
	    Y=X,
	    rename(F0,[X|Vars],F)
	;
	    rename(F0,Vars,F1),
	    term_variables(F1,VarsF1),
	    exchange(VarsF1,X,Y,VarsF2),
	    copy_term(VarsF1-F1,VarsF2-F)).
rename(exists(X,F0),Vars,exists(Y,F)) :-
	(notin(Vars,X) ->
	    Y=X,
	    rename(F0,[X|Vars],F)
	;
	    rename(F0,Vars,F1),
	    term_variables(F1,VarsF1),
	    exchange(VarsF1,X,Y,VarsF2),
	    copy_term(VarsF1-F1,VarsF2-F)).

exchange([],_,_,[]).
exchange([A|As],X,Y,Zs) :-
	(X==A ->
	    Zs=[Y|As]
	;
	    Zs=[A|Ts],
	    exchange(As,X,Y,Ts)).

%%%%
disjunctive_normal_form(forall(X,F0),forall(X,F))  :-
	disjunctive_normal_form(F0,F).
disjunctive_normal_form(exists(X,F0),exists(X,F))  :-
	disjunctive_normal_form(F0,F).
disjunctive_normal_form(ac(A),ac(B)) :-
	simplify_constraints(ac(A),ac(B)).
disjunctive_normal_form((A,As),DNF) :-
	disjunctive_normal_form2((A,As),DNF).
disjunctive_normal_form((A;As),DNF) :-
	disjunctive_normal_form2((A;As),DNF).

disjunctive_normal_form2(A,C) :-	% base case
	conjunction_atoms(A,B),!,
	simplify_constraints(B,C).
disjunctive_normal_form2((A , (B ; C)), G ) :- !,	% distribution
	disjunctive_normal_form2((A , B),D),
	disjunctive_normal_form2((A , C),E),
	disjunction_conjunctions(D,E,F),
	simplify_dnf(F,G).
disjunctive_normal_form2(((A ; B) , C), G) :-!,	% distribution
	disjunctive_normal_form2((A , C),D),
	disjunctive_normal_form2((B , C),E),
	disjunction_conjunctions(D,E,F),
	simplify_dnf(F,G).
disjunctive_normal_form2((A ; B),F) :-	% disjunction
	disjunctive_normal_form2(A,C),
	disjunctive_normal_form2(B,D),
	disjunction_conjunctions(C,D,E),
	simplify_dnf(E,F).
disjunctive_normal_form2((A , B),E) :-	% other cases
	disjunctive_normal_form2(A,C),
	disjunctive_normal_form2(B,D),
	disjunctive_normal_form2((C , D),E).
	  
conjunction_atoms(ac(A),ac(A)).
conjunction_atoms((C,Cs),Ds) :-
	conjunction_atoms(C,Cs,Ds).

conjunction_atoms(ac(A),Cs,(ac(A),Ds)) :-
	conjunction_atoms(Cs,Ds).
conjunction_atoms((A,As),Cs,Ds) :-
	conjunction_atoms(A,(As,Cs),Ds).

disjunction_conjunctions(ac(A),B,(ac(A);B)).
disjunction_conjunctions((A,As),B,((A,As);B)).
disjunction_conjunctions((A;B),C,(A;E)) :-
	disjunction_conjunctions(B,C,E).


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  qproject(+Xs,S,-Ys,-T) if
%%		env(Ys,T) is a projected copy of env(Xs,S)

qproject(Xs, Cxs, Ys, Cys) :-
    copy_term(Xs-Cxs, CopyXs-CopyCxs),
    (tell_cs(CopyCxs) -> 
	prepare_dump(CopyXs, Zs, Ys, Eqns),
	dump(Zs, Ys, Ds),
	append(Eqns, Ds, Cys)
    ;
	Ys=CopyXs,
	Cys=[0>1]).

tell_cs([]).
tell_cs([C|Cs]) :-  {C}, tell_cs(Cs).

prepare_dump([], [], [], []).
prepare_dump([X|Xs], [Z|Zs],[Y|Ys], Cs) :-
	(var(X) ->
	    Z=X,
	    Cs=Ds,
	    prepare_dump(Xs,Zs,Ys,Ds)
	;
	    Cs=[Y=X|Ds],
	    prepare_dump(Xs,Zs,Ys,Ds)).



:- begin_tests(tests).

test(a0) :- qdecproc(ac(1+1=2),ac(0<1)).
test(a1) :- qdecproc(ac(1+1=0),ac(0>1)).

test(b0) :- qdecproc(exists(X,exists(Y,ac(Y > X))),ac(0<1)).
test(b1) :- qdecproc(forall(X,exists(Y,ac(Y > X))),ac(0<1)).
test(b2) :- qdecproc(exists(X,forall(Y,ac(Y > X))),ac(0>1)).
test(b3) :- qdecproc(forall(X,forall(Y,ac(Y > X))),ac(0>1)).

test(c0) :- qdecproc(exists(Y,(ac(X=2*Y+1),ac(Y>1))),ac(X>3)).
test(c1) :- qdecproc(exists(X,exists(Y,(ac(X=2*Y+1),ac(Y>1)))),ac(0<1)).
test(c2) :- qdecproc(forall(X,forall(Y,(ac(X=Y) ; ac(X =\=Y)))),ac(0<1)).
test(c3) :- qdecproc(forall(X,(ac(X=Y) ; ac(X =\=Y))),ac(0<1)).

test(d0) :- qdecproc(forall(X,forall(Y, (ac(X > Y) -> exists(Z, (ac(X > Z), ac(Z > Y)))))),ac(0<1)).
test(d1) :- qdecproc(forall(X,forall(Y, (ac(X >= Y) -> exists(Z, (ac(X > Z), ac(Z > Y)))))),ac(0>1)).
test(d2) :- qdecproc(forall(X,forall(Y, (ac(X >= Y) -> exists(Z, (ac(X >= Z), ac(Z >= Y)))))),ac(0<1)).
test(d3) :- qdecproc(forall(X,forall(Y, (ac(X =\= Y) -> 
			exists(Z, ((ac(X > Z), ac(Z > Y));(ac(Y > Z),ac(Z > X))))))),ac(0<1)).
test(d4) :- qdecproc(forall(X,forall(Y, 
		       (ac(X = Y) ;exists(Z, ((ac(X > Z), ac(Z > Y));(ac(Y > Z),ac(Z > X))))))),ac(0<1)).
test(d5) :- qdecproc(forall(X,ac(X=_Y)),ac(0>1)).
test(d6) :- qdecproc(exists(X,ac(X=_Y)),ac(0<1)).
test(d7) :- qdecproc(exists(X,ac(X=\=_Y)),ac(0<1)).

test(e0) :- qdecproc(forall(X,(ac(X=1) -> exists(_Y,ac(X=1)))),ac(0<1)).

test(f0) :- qdecproc(forall(_454,forall(_474,forall(_494,forall(_548,forall(_584,forall(_642,(
	      ac(_454=_548+1),ac(_584=2*_474),ac(_642=_494),ac(_474>=1),ac(_494>=2)->
		forall(_454,forall(_474,forall(_494,(ac(0<1)->
			exists(_548,exists(_584,exists(_642,
			        (ac(_454=_548+1),ac(_584=2*_474),ac(_642=_494),
				 ac(_474>=1),ac(_494>=2)))))))))))))))),ac(0>1)).
:- end_tests(tests).

%:- run_tests.




