1% FINITE and INFINITE DOMAINS		
    2% 910527 ECRC thom fruehwirth
    3% 910913 modified 
    4% 920409 element/3 added
    5% 920616 more CHIP predicates added
    6% 930726 started porting to CHR release
    7% 931014 mult/3 added for CHIC user meeting
    8% 931201 ported to CHR release
    9% 931208 removed special case of integer domain
   10% 940304 element/3 constraint loop fixed
   11% 961017 Christian Holzbaur SICStus mods
   12% 980714 Thom Fruehwirth, some updates reagrding alread_in*
   13% 060527 Marco Gavanelli, use uf unnamed variables to remove "singleton variables" warning
   14
   15% just quick port from Eclipse CHR library version
   16% does not take advantage of Sicstus CHR library features!
   17		
   18% Simplifies domains together with inequalities and some more CHIP predicates:
   19% 	element/3, atmost/3, alldistinct/1, circuit/1 and mult/3
   20% It also includes paired (!) domains (see element constraint)
   21
   22
   23
   24:-module(domains,[(::)/2,
   25		  le/2,
   26		  lt/2,
   27		  plus/3,
   28		  ne/2,
   29		  ge/2,
   30		  gt/2,
   31		  alldistinct/1]).   32
   33operator(700,xfx,'::').
   34
   35:- use_module(library(chr)).   36%:- use_module(library('chr/getval')).
   37:- use_module(library(lists), [member/2,last/2]).   38
   39:- use_module( library(ordsets),
   40	[
   41          list_to_ord_set/2,
   42	  ord_intersection/3
   43        ]).   44
   45handler domain.
   46
   47option(already_in_store, on).   
   48option(already_in_heads, off).   % see pragma already_in_heads
   49option(check_guard_bindings, off).
   50
   51
   52% for domain constraints
   53
   54operator(600,xfx,'..').
   55operator(600,xfx,':').  % clash with module operator?
   56
   57% for inequality constraints
   58operator(700,xfx,lt).
   59operator(700,xfx,le).
   60operator(700,xfx,gt).
   61operator(700,xfx,ge).
   62operator(700,xfx,ne).
   63
   64% X::Dom - X must be element of the finite or infinite domain Dom
   65
   66% Domains can be either numbers (including arithemtic expressions)
   67% or arbitrary ground terms (!), the domain is set with setval(domain,Kind),
   68% where Kind is either number or term. Default for Kind is term.
   69
   70%:- setval(domain,term). 	% set default
   71
   72
   73% INEQUALITIES ===============================================================
   74% inequalities over numbers (including arithmetic expressions) or terms
   75
   76constraints lt/2,le/2,ne/2.
   77
   78gt(A,B) :- lt(B,A).				% constraints gt/2,ge/2
   79ge(A,B) :- le(B,A).
   80% some basic simplifications
   81lt(A,A) <=> fail.
   82le(A,A) <=> true.
   83ne(A,A) <=> fail.
   84lt(A,B),lt(B,A) <=> fail.
   85le(A,B), le(B,A) <=> A=B.
   86ne(A,B)
   87\
   88ne(B,A) <=> true.
   89% for number domain, allow arithmetic expressions in the arguments
   90lt(A,B) <=> domain(number),ground(A),\+ number(A) | A1 is A, lt(A1,B).
   91lt(B,A) <=> domain(number),ground(A),\+ number(A) | A1 is A, lt(B,A1).
   92le(A,B) <=> domain(number),ground(A),\+ number(A) | A1 is A, le(A1,B).
   93le(B,A) <=> domain(number),ground(A),\+ number(A) | A1 is A, le(B,A1).
   94ne(A,B) <=> domain(number),ground(A),\+ number(A) | A1 is A, ne(A1,B).
   95ne(B,A) <=> domain(number),ground(A),\+ number(A) | A1 is A, ne(B,A1).
   96% use built-ins to solve the predicates if arguments are known
   97lt(A,B) <=> ground(A),ground(B) | (domain(number) -> A < B ; A @< B).
   98le(A,B) <=> ground(A),ground(B) | (domain(number) -> A =< B ; A @=< B).
   99ne(A,B) <=> ground(A),ground(B) | (domain(number) -> A =\= B ; A \== B).
  100
  101
  102
  103% FINITE and INFINITE DOMAINS ================================================
  104
  105constraints (::)/2.
  106
  107% enforce groundness of domain expression
  108 ::(X,Dom) <=> nonground(Dom) | 
  109        raise_exception( instantiation_error(::(X,Dom),2)).
  110
  111constraints labeling/0.
  112
  113labeling, (X::[Y|L]) # Ph <=> 
  114	member(X,[Y|L]), labeling
  115    pragma passive(Ph).
  116
  117% binary search by splitting domain in halves
  118labeling, (X::Min:Max) # Ph <=> domain(number),Min+0.5<Max |  % ensure termination
  119	(integer(Min),integer(Max) ->  % assume we have integer domain
  120	Mid is (Min+Max)//2, Next is Mid+1
  121	;
  122	Mid is (Min+Max)/2, Next=Mid   % splitted domains overlap at Mid for floats
  123	),
  124	(
  125	X::Min:Mid
  126	;
  127	X::Next:Max
  128	% ;
  129	% Min+1>Max,	% for floats only, to get X also bound
  130	% X=Min		% or X=Max etc.
  131	),
  132	labeling
  133    pragma passive(Ph).
  134
  135	nonground(X) :- ground(X), !, fail.
  136        nonground(_).
  137
  138	domain(Kind) :- getval(domain,Kind).
  139
  140% CHIP list shorthand for domain variables
  141% list must be known (end in the empty list)
  142
  143 [X|L]::Dom <=> makedom([X|L],Dom).
  144
  145	makedom([],_) :- true.
  146	makedom([X|L],D) :- 
  147                nonvar(L),
  148		X::D,
  149		makedom(L,D).
  150
  151
  152% Consecutive integer domain ---------------------------------------------
  153% X::Min..Max - X is an integer between the numbers Min and Max (included)
  154% constraint is mapped to enumeration domain constraint
  155 X::Min..Max <=> 
  156        Min0 is Min, 
  157        (Min0=:=round(float(Min0)) -> Min1 is integer(Min0) ; Min1 is integer(Min0+1)),
  158	Max1 is integer(Max),
  159	interval(Min1,Max1,L), 
  160	X::L.
  161
  162 	interval(M,N,[M|Ns]):- 
  163		M<N, 
  164		!, 
  165		M1 is M+1, 
  166		interval(M1,N,Ns).
  167	interval(N,N,[N]).
  168
  169
  170% Enumeration domain -----------------------------------------------------
  171
  172% X::Dom - X must be a ground term in the ascending sorted ground list Dom
  173 X::[A|L] <=> list_to_ord_set([A|L],SL), SL\==[A|L] | X::SL.
  174% for number domain, allow arithmetic expressions in domain
  175 X::[A|L] <=> domain(number), member(X,[A|L]), \+ number(X) |
  176		eval_list([A|L],L1),list_to_ord_set(L1,L2), X::L2.
  177
  178	eval_list([],[]).
  179	eval_list([X|L1],[Y|L2]):-
  180		Y is X,
  181		eval_list(L1,L2).
  182
  183% special cases
  184 X::[] <=> fail.				
  185 X::[Y] <=> X=Y.
  186 X::[A|L] <=> ground(X) | (member(X,[A|L]) -> true).
  187
  188% intersection of domains for the same variable
  189% without pragma already_in_heads, needs already_in_store
  190 (X::[A1|L1]) \ (X::[A2|L2]) <=> 
  191    ord_intersection([A1|L1],[A2|L2],L),
  192    L \== [A2|L2]
  193    | 
  194    X::L.
  195
  196% interaction with inequalities
  197 (X::[A|L]) \ (X ne Y) <=> integer(Y), remove(Y,[A|L],L1) | X::L1.
  198 (X::[A|L]) \ (Y ne X) <=> integer(Y), remove(Y,[A|L],L1) | X::L1.
  199
  200 X::[A|L], Y le X ==> ground(Y), remove_lower(Y,[A|L],L1) | X::L1.
  201 X::[A|L], X le Y ==> ground(Y), remove_higher(Y,[A|L],L1) | X::L1.
  202 X::[A|L], Y lt X ==> ground(Y), remove_lower(Y,[A|L],L1),remove(Y,L1,L2) | X::L2.
  203 X::[A|L], X lt Y ==> ground(Y), remove_higher(Y,[A|L],L1),remove(Y,L1,L2) | X::L2.
  204
  205% interaction with interval domain
  206 X::[A|L], X::Min:Max ==> remove_lower(Min,[A|L],L1),remove_higher(Max,L1,L2) | X::L2.
  207
  208% propagation of bounds
  209 X le Y, Y::[A|L]   ==> var(X) | last([A|L],Max), X le Max.
  210 X le Y, X::[Min|_] ==> var(Y) | Min le Y.
  211 X lt Y, Y::[A|L]   ==> var(X) | last([A|L],Max), X lt Max.
  212 X lt Y, X::[Min|_] ==> var(Y) | Min lt Y.
  213
  214% Interval domain ---------------------------------------------------------
  215% X::Min:Max - X must be a ground term between Min and Max (included)
  216% for number domain, allow for arithmetic expressions ind omain
  217% for integer domains, X::Min..Max should be used
  218 X::Min:Max <=> domain(number), \+ (number(Min),number(Max)) |
  219		Min1 is Min, Max1 is Max, X::Min1:Max1.
  220% special cases
  221 X::Min:Min <=> X=Min.
  222 X::Min:Max <=> (domain(number) -> Min>Max ; Min@>Max) | fail.
  223 X::Min:Max <=> ground(X) | 
  224		(domain(number) -> Min=<X,X=<Max ; Min@=<X,X@=<Max).
  225% intersection of domains for the same variable
  226% without pragma already_in_heads, needs already_in_store
  227 (X::Min1:Max1) \ (X::Min2:Max2) <=> maximum(Min1,Min2,Min),
  228	                         minimum(Max1,Max2,Max),
  229                                 (Min \== Min2  ; Max \== Max2 ) |
  230		X::Min:Max.
  231
  232	minimum(A,B,C):- (domain(number) -> A<B ; A@<B) -> A=C ; B=C.
  233	maximum(A,B,C):- (domain(number) -> A<B ; A@<B) -> B=C ; A=C.
  234
  235% interaction with inequalities
  236 (X::Min:Max) \ (X ne Y) <=> ground(Y),
  237	(domain(number) -> (Y<Min;Y>Max) ; (Y@<Min;Y@>Max)) | true.
  238 (X::Min:Max) \ (Y ne X) <=> ground(Y),
  239	(domain(number) -> (Y<Min;Y>Max) ; (Y@<Min;Y@>Max)) | true.
  240 (X::Min1:Max) \ (Min2 le X) <=> ground(Min2) , maximum(Min1,Min2,Min) | X::Min:Max.
  241 (X::Min:Max1) \ (X le Max2) <=> ground(Max2) , minimum(Max1,Max2,Max) | X::Min:Max.
  242 (X::Min1:Max) \ (Min2 lt X) <=> ground(Min2) , maximum(Min1,Min2,Min) |
  243		X::Min:Max, X ne Min.
  244 (X::Min:Max1) \ (X lt Max2) <=> ground(Max2) , minimum(Max1,Max2,Max) |
  245		X::Min:Max, X ne Max.
  246% propagation of bounds
  247 X le Y, Y::Min:Max ==> var(X) | X le Max.
  248 X le Y, X::Min:Max ==> var(Y) | Min le Y.
  249 X lt Y, Y::Min:Max ==> var(X) | X lt Max.
  250 X lt Y, X::Min:Max ==> var(Y) | Min lt Y.
  251
  252
  253
  254% MULT/3 EXAMPLE EXTENSION ==================================================
  255% mult(X,Y,C) - integer X multiplied by integer Y gives the integer constant C.
  256
  257constraints mult/3.
  258
  259mult(X,Y,C) <=> ground(X) | (X=:=0 -> C=:=0 ; 0=:=C mod X, Y is C//X).
  260mult(Y,X,C) <=> ground(X) | (X=:=0 -> C=:=0 ; 0=:=C mod X, Y is C//X).
  261mult(X,Y,C), X::MinX:MaxX ==> 
  262	%(Dom=MinX:MaxX -> true ; Dom=[MinX|L],last(L,MaxX)),
  263	MinY is (C-1)//MaxX+1,
  264        MaxY is C//MinX,
  265	Y::MinY:MaxY.
  266mult(Y,X,C), X::MinX:MaxX ==>
  267	%(Dom=MinX:MaxX -> true ; Dom=[MinX|L],last(L,MaxX)),
  268	MinY is (C-1)//MaxX+1,
  269        MaxY is C//MinX,
  270	Y::MinY:MaxY.
  271
  272/*
  273:- mult(X,Y,156),[X,Y]::2:156,X le Y.
  274
  275X = X_g307
  276Y = Y_g331
  277 
  278Constraints:
  279(1) mult(X_g307, Y_g331, 156)
  280(7) Y_g331 :: 2 : 78
  281(8) X_g307 :: 2 : 78
  282(10) X_g307 le Y_g331
  283
  284yes.
  285:- mult(X,Y,156),[X,Y]::2:156,X le Y,labeling.
  286
  287X = 12
  288Y = 13     More? (;) 
  289
  290X = 6
  291Y = 26     More? (;) 
  292
  293X = 4
  294Y = 39     More? (;) 
  295
  296X = 2
  297Y = 78     More? (;) 
  298
  299X = 3
  300Y = 52     More? (;) 
  301
  302no (more) solution.
  303*/
  304
  305
  306
  307
  308% CHIP ELEMENT/3 ============================================================
  309% translated to "pair domains", a very powerful extension of usual domains
  310% this version does not work with arithmetic expressions!
  311
  312element(I,VL,V):- length(VL,N),interval(1,N,IL),gen_pair(IL,VL,BL), I-V::BL.
  313
  314	gen_pair([],[],[]).
  315	gen_pair([A|L1],[B|L2],[A-B|L3]):-
  316		gen_pair(L1,L2,L3).
  317
  318% special cases
  319 I-I::L <=> setof(X,member(X-X,L),L1), I::L1.
  320 I-V::L <=> ground(I) | setof(X,member(I-X,L),L1), V::L1.
  321 I-V::L <=> ground(V) | setof(X,member(X-V,L),L1), I::L1.
  322% intersections
  323 X::[A|L1], X-Y::L2 <=> intersect(I::[A|L1],I-V::L2,I-V::L3),
  324			length(L2,N2),length(L3,N3),N2>N3 | X-Y::L3.
  325 Y::[A|L1], X-Y::L2 <=> intersect(V::[A|L1],I-V::L2,I-V::L3),
  326			length(L2,N2),length(L3,N3),N2>N3 | X-Y::L3.
  327 X-Y::L1, Y-X::L2 <=> intersect(I-V::L1,V-I::L2,I-V::L3) | X-Y::L3.
  328 X-Y::L1, X-Y::L2 <=> intersect(I-V::L1,I-V::L2,I-V::L3) | X-Y::L3 pragma already_in_heads.
  329
  330    intersect(A::L1,B::L2,C::L3):- setof(C,A^B^(member(A,L1),member(B,L2)),L3).
  331
  332% inequalties with two common variables
  333 Y lt X, X-Y::L <=> A=R-S,setof(A,(member(A,L),R@< S),L1) | X-Y::L1.
  334 X lt Y, X-Y::L <=> A=R-S,setof(A,(member(A,L),S@< R),L1) | X-Y::L1.
  335 Y le X, X-Y::L <=> A=R-S,setof(A,(member(A,L),R@=<S),L1) | X-Y::L1.
  336 X le Y, X-Y::L <=> A=R-S,setof(A,(member(A,L),S@=<R),L1) | X-Y::L1.
  337 Y ne X, X-Y::L <=> A=R-S,setof(A,(member(A,L),R\==S),L1) | X-Y::L1.
  338 X ne Y, X-Y::L <=> A=R-S,setof(A,(member(A,L),S\==R),L1) | X-Y::L1.
  339% propagation between paired domains (path-consistency)
  340% X-Y::L1, Y-Z::L2 ==> intersect(A-B::L1,B-C::L2,A-C::L), X-Z::L.
  341% X-Y::L1, Z-Y::L2 ==> intersect(A-B::L1,C-B::L2,A-C::L), X-Z::L.
  342% X-Y::L1, X-Z::L2 ==> intersect(I-V::L1,I-W::L2,V-W::L), Y-Z::L.
  343% propagation to usual unary domains
  344 X-Y::L ==> A=R-S,setof(R,A^member(A,L),L1), X::L1,
  345	          setof(S,A^member(A,L),L2), Y::L2.
  346
  347
  348
  349% ATMOST/3 ===================================================================
  350
  351atmost(N,List,V):-length(List,K),atmost(N,List,V,K).
  352
  353constraints atmost/4.
  354
  355atmost(N,List,V,K) <=> K=<N | true.
  356atmost(0,List,V,K) <=> (ground(V);ground(List)) | outof(V,List).
  357atmost(N,List,V,K) <=> K>N,ground(V),delete_ground(X,List,L1) |
  358		(X==V -> N1 is N-1 ; N1=N),K1 is K-1, atmost(N1,L1,V,K1).
  359
  360	delete_ground(X,List,L1):- delete(X,List,L1),ground(X),!.
  361
  362delete( X, [X|Xs], Xs).
  363delete( Y, [X|Xs], [X|Xt]) :-
  364	delete( Y, Xs, Xt).
  365
  366
  367% ALLDISTINCT/1 ===============================================================
  368% uses ne/2 constraint
  369
  370constraints alldistinct/1.
  371
  372alldistinct([]) <=> true.
  373alldistinct([X]) <=> true.
  374alldistinct([X,Y]) <=> X ne Y.
  375alldistinct([A|L]) <=> delete_ground(X,[A|L],L1) | outof(X,L1),alldistinct(L1).
  376
  377alldistinct([]).
  378alldistinct([X|L]):-
  379	outof(X,L),
  380	alldistinct(L).
  381
  382outof(_,[]).
  383outof(X,[Y|L]):-
  384	X ne Y,
  385	outof(X,L).
  386
  387constraints alldistinct1/2.
  388		
  389alldistinct1(R,[]) <=> true.
  390alldistinct1(R,[X]), X::[A|L] <=> ground(R) | 
  391			remove_list(R,[A|L],T), X::T.
  392alldistinct1(R,[X]) <=> (ground(R);ground(X)) | outof(X,R).	
  393alldistinct1(R,[A|L]) <=> ground(R),delete_ground(X,[A|L],L1) | 
  394			(member(X,R) -> fail ; alldistinct1([X|R],L1)).
  395
  396
  397
  398% CIRCUIT/1 =================================================================
  399
  400% constraints circuit1/1, circuit/1.
  401% uses list domains and ne/2
  402
  403
  404% lazy version
  405
  406circuit1(L):-length(L,N),N>1,circuit1(N,L).
  407
  408circuit1(2,[2,1]).
  409circuit1(N,L):- N>2,
  410		interval(1,N,D),
  411		T=..[f|L],
  412		domains1(1,D,L),
  413		alldistinct1([],L),
  414		no_subtours(N,1,T,[]).	
  415
  416domains1(_,_,[]).
  417domains1(N,D,[X|L]):- 
  418		remove(N,D,DX),
  419		X::DX,
  420		N1 is N+1,
  421		domains1(N1,D,L).
  422
  423no_subtours(0,_,_,_):- !.
  424no_subtours(K,N,L,R):- 
  425	outof(N,R),
  426	(var(N) -> freeze(N,no_subtours1(K,N,L,R)) ; no_subtours1(K,N,L,R)).
  427% no_subtours(K,N,T,R) \ no_subtours(K1,N,T,_) <=> K<K1 | true.
  428
  429	no_subtours1(K,N,L,R):- 
  430		K>0,K1 is K-1,arg(N,L,A),no_subtours(K1,A,L,[N|R]).
  431
  432
  433% eager version
  434
  435circuit(L):- length(L,N),N>1,circuit(N,L).
  436
  437circuit(2,[2,1]).
  438%circuit(3,[2,3,1]).
  439%circuit(3,[3,1,2]).
  440circuit(N,L):- 	N>2,
  441		interval(1,N,D),
  442		T=..[f|L],
  443		N1 is N-1,
  444		domains(1,D,L,T,N1),
  445		alldistinct(L).		
  446
  447domains(_,_,[],_,_).
  448domains(N,D,[X|L],T,K):- 
  449		remove(N,D,DX),
  450		X::DX,
  451		N1 is N+1,
  452		no_subtours(K,N,T,[]),		% unfolded
  453		%no_subtours1(K,X,T,[N]),		
  454		domains(N1,D,L,T,K).
  455
  456
  457
  458
  459% remove*/3 auxiliary predicates =============================================
  460
  461remove(A,B,C):- 
  462	delete(A,B,C) -> true ; B=C.
  463
  464remove_list(_,[],T):- !, T=[].
  465remove_list([],S,T):- S=T.
  466remove_list([X|R],[Y|S],T):- remove(X,[Y|S],S1),remove_list(R,S1,T).
  467
  468remove_lower(_,[],L1):- !, L1=[].
  469remove_lower(Min,[X|L],L1):-
  470	X@<Min,
  471	!,
  472	remove_lower(Min,L,L1).
  473remove_lower(Min,[X|L],[X|L1]):-
  474	remove_lower(Min,L,L1).
  475
  476remove_higher(_,[],L1):- !, L1=[].
  477remove_higher(Max,[X|L],L1):-
  478	X@>Max,
  479	!,
  480	remove_higher(Max,L,L1).
  481remove_higher(Max,[X|L],[X|L1]):-
  482	remove_higher(Max,L,L1).
  483
  484
  485/* I wrote this! And it works!! MarcoA */
  486
  487constraints plus/3.
  488
  489plus_1 @
  490	plus(X,Y,Z)
  491	<=>
  492	number(X),
  493	number(Y)
  494	|
  495	Z is X+Y.
  496
  497plus_2 @
  498	plus(X,Y,Z)
  499	<=>
  500	number(X),
  501	number(Z)
  502	|
  503	Y is Z-X.
  504
  505plus_1 @
  506	plus(X,Y,Z)
  507	<=>
  508	number(Y),
  509	number(Z)
  510	|
  511	X is Z-Y.
  512
  513
  514% end of handler domain.chr =================================================
  515% ===========================================================================