1/*
    2
    3    Part of CLP(Q,R) (Constraint Logic Programming over Rationals and Reals)
    4
    5    Author:        Leslie De Koninck
    6    E-mail:        Leslie.DeKoninck@cs.kuleuven.be
    7    WWW:           http://www.swi-prolog.org
    8		   http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09
    9    Copyright (C): 2006, K.U. Leuven and
   10		   1992-1995, Austrian Research Institute for
   11		              Artificial Intelligence (OFAI),
   12			      Vienna, Austria
   13
   14    This software is based on CLP(Q,R) by Christian Holzbaur for SICStus
   15    Prolog and distributed under the license details below with permission from
   16    all mentioned authors.
   17
   18    This program is free software; you can redistribute it and/or
   19    modify it under the terms of the GNU General Public License
   20    as published by the Free Software Foundation; either version 2
   21    of the License, or (at your option) any later version.
   22
   23    This program is distributed in the hope that it will be useful,
   24    but WITHOUT ANY WARRANTY; without even the implied warranty of
   25    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   26    GNU General Public License for more details.
   27
   28    You should have received a copy of the GNU Lesser General Public
   29    License along with this library; if not, write to the Free Software
   30    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
   31
   32    As a special exception, if you link this library with other files,
   33    compiled with a Free Software compiler, to produce an executable, this
   34    library does not by itself cause the resulting executable to be covered
   35    by the GNU General Public License. This exception does not however
   36    invalidate any other reasons why the executable file might be covered by
   37    the GNU General Public License.
   38*/
   39
   40
   41% attribute = t(CLP,type(_),strictness(_),lin(_),order(_),class(_),forward(_),
   42%		nonzero,target,keep_indep,keep)
   43
   44:- module(clpcd_itf,
   45	[
   46	    dump_linear/3,
   47	    dump_nonzero/3
   48	]).   49
   50:- use_module(library(neck)).   51:- use_module(library(clpcd/class)).   52:- use_module(library(clpcd/domain_ops)).   53:- use_module(library(clpcd/indep)).   54:- use_module(library(clpcd/bv)).   55:- use_module(library(clpcd/nf)).   56:- use_module(library(clpcd/store)).   57:- use_module(library(clpcd/solve)).   58:- init_expansors.   59
   60% ----------------------------------- dump ------------------------------------
   61
   62% dump_strict(FilteredStrictness,Nonstrict,Strict,Res)
   63%
   64% Unifies Res with either Nonstrict or Strict depending on FilteredStrictness.
   65% FilteredStrictness is the component of strictness related to the bound: 0
   66% means nonstrict, 1 means strict upperbound, 2 means strict lowerbound,
   67% 3 is filtered out to either 1 or 2.
   68
   69dump_strict(0,Result,_,Result).
   70dump_strict(1,_,Result,Result).
   71dump_strict(2,_,Result,Result).
   72
   73% dump_nz(CLP,V,H,I,Dump,DumpTail)
   74%
   75% Returns in Dump a representation of the nonzero constraint of variable V
   76% which has linear
   77% equation H + I.
   78
   79dump_nz(CLP,_,H,I) -->
   80	{
   81	    H = [l(_*K,_)|_],
   82	    div_d(CLP, 1, K, Kr),
   83	    eval_d(CLP, -Kr*I, I1),
   84	    mult_hom(H,CLP,Kr,H1),
   85	    nf2sum(H1,CLP,0,Sum)
   86	},
   87	[Sum =\= I1].
   88
   89% dump_var(Type,Var,I,H,Dump,DumpTail)
   90%
   91% Returns in Dump a representation of the linear constraint on variable
   92% Var which has linear equation H + I and has type Type.
   93
   94dump_v(t_none,CLP,V,I,H) -->
   95	!,
   96	(   {
   97		H = [l(W*K,_)],
   98		V == W,
   99                compare_d(CLP, =, 0, I),
  100                compare_d(CLP, =, K, 1)
  101	    }
  102	->  % indep var
  103	    []
  104	;   {nf2sum(H,CLP,I,Sum)},
  105	    [V = Sum]
  106	).
  107dump_v(t_L(L),CLP,V,I,H) -->
  108	!,
  109	dump_v(t_l(L),CLP,V,I,H).
  110% case lowerbound: V >= L or V > L
  111% say V >= L, and V = K*V1 + ... + I, then K*V1 + ... + I >= L
  112% and K*V1 + ... >= L-I and V1 + .../K = (L-I)/K
  113dump_v(t_l(L),CLP,V,I,H) -->
  114	!,
  115	{
  116	    H = [l(_*K,_)|_], % avoid 1 >= 0
  117	    get_attr(V,clpcd_itf,Att),
  118	    arg(3,Att,strictness(Strict)),
  119	    Sm is Strict /\ 2,
  120            div_d(CLP, 1, K, Kr),
  121	    eval_d(CLP, Kr*(L - I), Li),
  122	    mult_hom(H,CLP,Kr,H1),
  123	    nf2sum(H1,CLP,0,Sum),
  124	    (   compare_d(CLP, <, 0, K)
  125	    ->	dump_strict(Sm,Sum >= Li,Sum > Li,Result)
  126	    ;   dump_strict(Sm,Sum =< Li,Sum < Li,Result)
  127	    )
  128	},
  129	[Result].
  130dump_v(t_U(U),C,V,I,H) -->
  131	!,
  132	dump_v(t_u(U),C,V,I,H).
  133dump_v(t_u(U),CLP,V,I,H) -->
  134	!,
  135	{
  136	    H = [l(_*K,_)|_], % avoid 0 =< 1
  137	    get_attr(V,clpcd_itf,Att),
  138	    arg(3,Att,strictness(Strict)),
  139	    Sm is Strict /\ 1,
  140	    div_d(CLP, 1, K, Kr),
  141	    eval_d(CLP, Kr*(U-I), Ui),
  142	    mult_hom(H,CLP,Kr,H1),
  143	    nf2sum(H1,CLP,0,Sum),
  144	    (   compare_d(CLP, <, 0, K)
  145	    ->	dump_strict(Sm,Sum =< Ui,Sum < Ui,Result)
  146	    ;   dump_strict(Sm,Sum >= Ui,Sum > Ui,Result)
  147	    )
  148	},
  149	[Result].
  150dump_v(t_Lu(L,U),C,V,I,H) -->
  151	!,
  152	dump_v(t_l(L),C,V,I,H),
  153	dump_v(t_u(U),C,V,I,H).
  154dump_v(t_lU(L,U),C,V,I,H) -->
  155	!,
  156	dump_v(t_l(L),C,V,I,H),
  157	dump_v(t_u(U),C,V,I,H).
  158dump_v(t_lu(L,U),C,V,I,H) -->
  159	!,
  160	dump_v(t_l(L),C,V,I,H),
  161	dump_v(t_U(U),C,V,I,H).
  162dump_v(T,_,V,I,H) --> % should not happen
  163	[V:T:I+H].
  164
  165
  166% nf2sum(Lin,CLP,Sofar,Term)
  167%
  168% Transforms a linear expression into a sum
  169% (e.g. the expression [5,_,[l(X*2,OrdX),l(Y*-1,OrdY)]] gets transformed into 5 + 2*X - Y)
  170
  171nf2sum([],_,I,I).
  172nf2sum([X|Xs],CLP,I,Sum) :-
  173	(   compare_d(CLP, =, 0, I)
  174	->  X = l(Var*K,_),
  175 	    (   % K =:= 1.0
  176                compare_d(CLP, =, K, 1)
  177	    ->  hom2sum(Xs,CLP,Var,Sum)
  178	    ;   % K =:= -1.0
  179                compare_d(CLP, =, K, -1)
  180	    ->  hom2sum(Xs,CLP,-Var,Sum)
  181	    ;	hom2sum(Xs,CLP,K*Var,Sum)
  182	    )
  183	;   hom2sum([X|Xs],CLP,I,Sum)
  184 	).
  185
  186% hom2sum(Hom,_,Sofar,Term)
  187%
  188% Transforms a linear expression into a sum
  189% this predicate handles all but the first term
  190% (the first term does not need a concatenation symbol + or -)
  191% see also nf2sum/3
  192
  193hom2sum([],_,Term,Term).
  194hom2sum([l(Var*K,_)|Cs],CLP,Sofar,Term) :-
  195	(   % K =:= 1.0
  196            compare_d(CLP, =, K, 1)
  197	->  Next = Sofar + Var
  198	;   % K =:= -1.0
  199            compare_d(CLP, =, K, -1)
  200	->  Next = Sofar - Var
  201	;   % K < 0.0
  202	    compare_d(CLP, <, K, 0)
  203	->  eval_d(CLP, -K, Ka),
  204	    Next = Sofar - Ka*Var
  205	;   Next = Sofar + K*Var
  206	),
  207	hom2sum(Cs,CLP,Next,Term).
  208
  209dump_linear(V) -->
  210	{
  211	    get_attr(V,clpcd_itf,Att),
  212	    arg(1,Att,CLP),
  213	    arg(2,Att,type(Type)),
  214	    arg(4,Att,lin(Lin)),
  215	    !,
  216	    Lin = [I,_|H]
  217	},
  218	(   {
  219		Type=t_none
  220	    ;	arg(9,Att,n)
  221	    }
  222	->  []
  223	;   dump_v(t_none,CLP,V,I,H)
  224	),
  225	(   {
  226		Type=t_none,
  227		arg(9,Att,n) % attribute should not have changed by dump_v...
  228	    }
  229	->  % nonzero produces such
  230	    []
  231	;   dump_v(Type,CLP,V,I,H)
  232	).
  233dump_linear(_) --> [].
  234
  235dump_nonzero(V) -->
  236	{
  237	    get_attr(V,clpcd_itf,Att),
  238	    arg(1,Att,CLP),
  239	    arg(4,Att,lin(Lin)),
  240	    arg(8,Att,nonzero),
  241	    !,
  242	    Lin = [I,_|H]
  243	},
  244	dump_nz(CLP,V,H,I).
  245dump_nonzero(_) --> [].
  246
  247attr_unify_hook(t(CLP,n,n,n,n,n,n,n,_,_,_),Y) :-
  248	!,
  249	(   get_attr(Y,clpcd_itf,AttY),
  250	    \+ arg(1,AttY,CLP)
  251	->  throw(error(permission_error('mix CLP(CD) variables with',
  252		'CLP(CD) variables of other subdomain:',Y),context(_)))
  253	;   true
  254	).
  255attr_unify_hook(t(CLP,Ty,St,Li,Or,Cl,_,No,_,_,_),Y) :-
  256	(   get_attr(Y,clpcd_itf,AttY),
  257	    \+ arg(1,AttY,CLP)
  258	->  throw(error(permission_error('mix CLP(CD) variables with',
  259		'CLP(CD) variables of other subdomain:',Y),context(_)))
  260	;   true
  261	),
  262	do_checks(CLP,Y,Ty,St,Li,Or,Cl,No,Later),
  263	maplist(call,Later).
  264
  265do_checks(CLP,Y,Ty,St,Li,Or,Cl,No,Later) :-
  266    numbers_only(CLP,Y),
  267    verify_nonzero(No,CLP,Y),
  268    verify_type(Ty,CLP,St,Y,Later,[]),
  269    verify_lin(Or,CLP,Cl,Li,Y),
  270    maplist(call,Later).
  271
  272% verify_nonzero(Nonzero,CLP,Y)
  273%
  274% if Nonzero = nonzero, then verify that Y is not zero
  275% (if possible, otherwise set Y to be nonzero)
  276
  277verify_nonzero(nonzero,CLP,Y) :-
  278    (   var(Y)
  279    ->  (   get_attr(Y,clpcd_itf,Att)
  280	->  setarg(8,Att,nonzero)
  281	;   put_attr(Y,clpcd_itf,t(CLP,n,n,n,n,n,n,nonzero,n,n,n))
  282	)
  283    ;   compare_d(CLP, \=, 0, Y)
  284    ).
  285verify_nonzero(n,_,_). % X is not nonzero
  286
  287% verify_type(type(Type),strictness(Strict),Y,[OL|OLT],OLT)
  288%
  289% if possible verifies whether Y satisfies the type and strictness of X
  290% if not possible to verify, then returns the constraints that follow from
  291% the type and strictness
  292
  293verify_type(type(Type),CLP,strictness(Strict),Y) -->
  294    verify_type2(CLP,Y,Type,Strict).
  295verify_type(n,_,n,_) --> [].
  296
  297verify_type2(C,Y,TypeX,StrictX) -->
  298    {var(Y)},
  299    !,
  300    verify_type_var(TypeX,C,Y,StrictX).
  301verify_type2(C,Y,TypeX,StrictX) -->
  302    {verify_type_nonvar(TypeX,C,Y,StrictX)}.
  303
  304% verify_type_nonvar(Type,CLP,Nonvar,Strictness)
  305%
  306% verifies whether the type and strictness are satisfied with the Nonvar
  307
  308verify_type_nonvar(t_none,_,_,_).
  309verify_type_nonvar(t_l(L),C,Value,S) :- ilb(S,C,L,Value).
  310verify_type_nonvar(t_u(U),C,Value,S) :- iub(S,C,U,Value).
  311verify_type_nonvar(t_lu(L,U),C,Value,S) :-
  312    ilb(S,C,L,Value),
  313    iub(S,C,U,Value).
  314verify_type_nonvar(t_L(L),C,Value,S) :- ilb(S,C,L,Value).
  315verify_type_nonvar(t_U(U),C,Value,S) :- iub(S,C,U,Value).
  316verify_type_nonvar(t_Lu(L,U),C,Value,S) :-
  317    ilb(S,C,L,Value),
  318    iub(S,C,U,Value).
  319verify_type_nonvar(t_lU(L,U),C,Value,S) :-
  320    ilb(S,C,L,Value),
  321    iub(S,C,U,Value).
  322
  323% ilb(Strict,CLP,Lower,Value) & iub(Strict,CLP,Upper,Value)
  324%
  325% check whether Value is satisfiable with the given lower/upper bound and
  326% strictness.
  327% strictness is encoded as follows:
  328% 2 = strict lower bound
  329% 1 = strict upper bound
  330% 3 = strict lower and upper bound
  331% 0 = no strict bounds
  332
  333ilb(S,C,L,V) :-
  334    between(0, 3, S),
  335    ( S /\ 2 =:= 0
  336    ->Op = (=<) % non-strict
  337    ; Op = (<)  % strict
  338    ),
  339    neck,
  340    compare_d(C, Op, L, V).
  341
  342iub(S,C,U,V) :-
  343    between(0, 3, S),
  344    ( S /\ 1 =:= 0
  345    ->Op = (=<) % non-strict
  346    ; Op = (<)  % strict
  347    ),
  348    neck,
  349    compare_d(C, Op, V, U).
  350
  351%
  352% Running some goals after X=Y simplifies the coding. It should be possible
  353% to run the goals here and taking care not to put_atts/2 on X ...
  354%
  355
  356% verify_type_var(Type,Var,Strictness,[OutList|OutListTail],OutListTail)
  357%
  358% returns the inequalities following from a type and strictness satisfaction
  359% test with Var
  360
  361verify_type_var(t_none,_,_,_) --> [].
  362verify_type_var(t_l(L),C,Y,S) --> llb(S,C,L,Y).
  363verify_type_var(t_u(U),C,Y,S) --> lub(S,C,U,Y).
  364verify_type_var(t_lu(L,U),C,Y,S) -->
  365    llb(S,C,L,Y),
  366    lub(S,C,U,Y).
  367verify_type_var(t_L(L),C,Y,S) --> llb(S,C,L,Y).
  368verify_type_var(t_U(U),C,Y,S) --> lub(S,C,U,Y).
  369verify_type_var(t_Lu(L,U),C,Y,S) -->
  370    llb(S,C,L,Y),
  371    lub(S,C,U,Y).
  372verify_type_var(t_lU(L,U),C,Y,S) -->
  373    llb(S,C,L,Y),
  374    lub(S,C,U,Y).
  375
  376% llb(Strict,Lower,Value,[OL|OLT],OLT) and lub(Strict,Upper,Value,[OL|OLT],OLT)
  377%
  378% returns the inequalities following from the lower and upper bounds and the
  379% strictness see also lb and ub
  380llb(S,C,L,V) -->
  381    {S /\ 2 =:= 0 },
  382    !,
  383    [C:{L =< V}].
  384llb(_,C,L,V) -->
  385    [C:{L < V}].
  386
  387lub(S,C,U,V) -->
  388    {S /\ 1 =:= 0 },
  389     !,
  390     [C:{V =< U}].
  391lub(_,C,U,V) -->
  392    [C:{V < U}].
  393
  394%
  395% We used to drop X from the class/basis to avoid trouble with subsequent
  396% put_atts/2 on X. Now we could let these dead but harmless updates happen.
  397% In R however, exported bindings might conflict, e.g. 0 \== 0.0
  398%
  399% If X is indep and we do _not_ solve for it, we are in deep shit
  400% because the ordering is violated.
  401%
  402verify_lin(order(OrdX),CLP,class(Class),lin(LinX),Y) :-
  403    !,
  404    (   indep(CLP,LinX,OrdX)
  405    ->  detach_bounds_vlv(CLP,OrdX,LinX,Class,Y,NewLinX),
  406	% if there were bounds, they are requeued already
  407	class_drop(Class,Y),
  408	nf(-Y, CLP, NfY),
  409	deref(CLP,NfY,LinY),
  410	add_linear_11(CLP, NewLinX, LinY, Lind),
  411	(   nf_coeff_of(Lind,OrdX,_)
  412	->  % X is element of Lind
  413	    solve_ord_x(CLP, Lind, OrdX, Class)
  414	;   solve(CLP, Lind)	% X is gone, can safely solve Lind
  415	)
  416    ;   class_drop(Class,Y),
  417	nf(-Y, CLP, NfY),
  418	deref(CLP,NfY,LinY),
  419	add_linear_11(CLP, LinX, LinY, Lind),
  420	solve(CLP, Lind)
  421    ).
  422verify_lin(_,_,_,_,_)