1/*
    2
    3    Part of CLP(Q) (Constraint Logic Programming over Rationals)
    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:- module(clpcd_solve,
   41          [
   42              'solve_='/2,
   43              'solve_=\\='/2,
   44              backsubst/4,
   45              backsubst_delta/5,
   46              basis_add/2,
   47              bs/4,
   48              dec_step/3,
   49              deref/3,
   50              deref_var/3,
   51              export_binding/1,
   52              inc_step/3,
   53              intro_at/4,
   54              lb/4,
   55              lb/5,
   56              log_deref/5,
   57              pivot/6,
   58              pivot_a/5,
   59              pivot_b/5,
   60              rcbl/4,
   61              rcbl_status/7,
   62              reconsider/2,
   63              same_class/2,
   64              select_active_bound/2,
   65              solve/2,
   66              solve/6,
   67              solve_ord_x/4,
   68              ub/4,
   69              ub/5,
   70              unconstrained/5,
   71              var_with_def_intern/5
   72          ]).   73
   74:- use_module(library(clpcd/class)).   75:- use_module(library(clpcd/domain_ops)).   76:- use_module(library(clpcd/ordering)).   77:- use_module(library(clpcd/split)).   78:- use_module(library(clpcd/store)).   79:- use_module(library(ordsets)).   80
   81% ----------------------------------- deref -----------------------------------
   82%
   83
   84% deref(CLP,Lin,Lind)
   85%
   86% Makes a linear equation of the form [v(I,[])|H] into a solvable linear
   87% equation.
   88% If the variables are new, they are initialized with the linear equation X=X.
   89
   90deref(CLP,Lin,Lind) :-
   91	split(Lin,H,I),
   92	normalize_scalar(I,Nonvar),
   93	length(H,Len),
   94	log_deref(Len,CLP,H,[],Restd),
   95	add_linear_11(CLP, Nonvar, Restd, Lind).
   96
   97% log_deref(Len,CLP,[Vs|VsTail],VsTail,Res)
   98%
   99% Logarithmically converts a linear equation in normal form ([v(_,_)|_]) into a
  100% linear equation in solver form ([I,R,K*X|_]). Res contains the result, Len is
  101% the length of the part to convert and [Vs|VsTail] is a difference list
  102% containing the equation in normal form.
  103
  104log_deref(0,_,Vs,Vs,Lin) :-
  105	!,
  106	Lin = [0,0].
  107log_deref(1,C,[v(K,[X^1])|Vs],Vs,Lin) :-
  108	!,
  109	deref_var(C, X, Lx),
  110	mult_linear_factor(C, Lx, K, Lin).
  111log_deref(2,C,[v(Kx,[X^1]),v(Ky,[Y^1])|Vs],Vs,Lin) :-
  112	!,
  113	deref_var(C,X,Lx),
  114	deref_var(C,Y,Ly),
  115	add_linear_ff(C, Lx, Kx, Ly, Ky, Lin).
  116log_deref(N,C,V0,V2,Lin) :-
  117	P is N >> 1,
  118	Q is N - P,
  119	log_deref(P,C,V0,V1,Lp),
  120	log_deref(Q,C,V1,V2,Lq),
  121	add_linear_11(C, Lp, Lq, Lin).
  122
  123% deref_var(CLP,X,Lin)
  124%
  125% Returns the equation of variable X. If X is a new variable, a new equation
  126% X = X is made.
  127
  128deref_var(CLP,X,Lin) :-
  129	(   get_attr(X,clpcd_itf,Att)
  130	->  (   \+ arg(1,Att,CLP)
  131	    ->  throw(error(permission_error('mix CLPCD variables of',
  132		'different domains:',X),context(_)))
  133	    ;   arg(4,Att,lin(Lin))
  134	    ->  true
  135	    ;   setarg(2,Att,type(t_none)),
  136		setarg(3,Att,strictness(0)),
  137		Lin = [0,0,l(X*1,Ord)],
  138		setarg(4,Att,lin(Lin)),
  139		setarg(5,Att,order(Ord))
  140	    )
  141	;   Lin = [0,0,l(X*1,Ord)],
  142	    put_attr(X,clpcd_itf,t(CLP,type(t_none),strictness(0),
  143		lin(Lin),order(Ord),n,n,n,n,n,n))
  144	).
  145
  146% var_with_def_intern(Type,CLP,Var,Lin,Strictness)
  147%
  148% Makes Lin the linear equation of new variable Var, makes all variables of
  149% Lin, and Var of the same class and bounds Var by type(Type) and
  150% strictness(Strictness)
  151
  152var_with_def_intern(Type,CLP,Var,Lin,Strict) :-
  153	put_attr(Var,clpcd_itf,t(CLP,type(Type),strictness(Strict),lin(Lin),
  154	    order(_),n,n,n,n,n,n)),	% check uses
  155	Lin = [_,_|Hom],
  156	get_or_add_class(Var,Class),
  157	same_class(Hom,Class).
  158
  159% -----------------------------------------------------------------------------
  160
  161% export_binding(Lst)
  162%
  163% Binds variables X to Y where Lst contains elements of the form [X-Y].
  164
  165export_binding([]).
  166export_binding([X-Y|Gs]) :-
  167	Y = X,
  168	export_binding(Gs).
  169
  170% 'solve_='(CLP,Nf)
  171%
  172% Solves linear equation Nf = 0 where Nf is in normal form.
  173
  174'solve_='(CLP,Nf) :-
  175	deref(CLP,Nf,Nfd),	% dereferences and turns Nf into solvable form Nfd
  176	solve(CLP,Nfd).
  177
  178% 'solve_=\\='(CLP,Nf)
  179%
  180% Solves linear inequality Nf =\= 0 where Nf is in normal form.
  181
  182'solve_=\\='(CLP,Nf) :-
  183	deref(CLP,Nf,Lind),	% dereferences and turns Nf into solvable form Lind
  184	Lind = [Inhom,_|Hom],
  185	(   Hom = []
  186	->  compare_d(CLP, \=, 0, Inhom)
  187	;   % make new variable Nz = Lind
  188	    var_with_def_intern(t_none,CLP,Nz,Lind,0),
  189	    % make Nz nonzero
  190	    get_attr(Nz,clpcd_itf,Att),
  191	    setarg(8,Att,nonzero)
  192	).
  193
  194%
  195% Status = {optimum,unlimited(Indep,DepT),applied}
  196% If Status = optimum, the tables have not been changed at all.
  197% Searches left to right, does not try to find the 'best' pivot
  198% Therefore we might discover unboundedness only after a few pivots
  199%
  200
  201
  202dec_step_cont([],_,optimum,Cont,Cont).
  203dec_step_cont([l(V*K,OrdV)|Vs],CLP,Status,ContIn,ContOut) :-
  204	get_attr(V,clpcd_itf,Att),
  205	arg(2,Att,type(W)),
  206	arg(6,Att,class(Class)),
  207	(   dec_step_2_cont(W,CLP,l(V*K,OrdV),Class,Status,ContIn,ContOut)
  208	->  true
  209	;   dec_step_cont(Vs,CLP,Status,ContIn,ContOut)
  210	).
  211
  212inc_step_cont([], _, optimum, Cont, Cont).
  213inc_step_cont([l(V*K,OrdV)|Vs], CLP, Status, ContIn, ContOut) :-
  214	get_attr(V,clpcd_itf,Att),
  215	arg(2,Att,type(W)),
  216	arg(6,Att,class(Class)),
  217	(   inc_step_2_cont(W, CLP, l(V*K,OrdV), Class, Status, ContIn, ContOut)
  218	->  true
  219	;   inc_step_cont(Vs, CLP, Status, ContIn, ContOut)
  220	).
  221
  222dec_step_2_cont(t_U(U),CLP,l(V*K,OrdV),Class,Status,ContIn,ContOut) :-
  223	compare_d(CLP, <, 0, K),
  224	(   lb(Class,CLP,OrdV,Vub-Vb-_)
  225	->  % found a lower bound
  226	    Status = applied,
  227	    pivot_a(CLP,Vub,V,Vb,t_u(U)),
  228	    replace_in_cont(ContIn,Vub,V,ContOut)
  229	;   Status = unlimited(V,t_u(U)),
  230	    ContIn = ContOut
  231	).
  232dec_step_2_cont(t_lU(L,U),CLP,l(V*K,OrdV),Class,applied,ContIn,ContOut) :-
  233	compare_d(CLP, <, 0, K),
  234	eval_d(CLP, L - U, Init),
  235	class_basis(Class,Deps),
  236	lb(Deps,CLP,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_),
  237	pivot_b(CLP,Vub,V,Vb,t_lu(L,U)),
  238	replace_in_cont(ContIn,Vub,V,ContOut).
  239dec_step_2_cont(t_L(L),CLP,l(V*K,OrdV),Class,Status,ContIn,ContOut) :-
  240	compare_d(CLP, >, 0, K),
  241	(   ub(Class,CLP,OrdV,Vub-Vb-_)
  242	->  Status = applied,
  243	    pivot_a(CLP,Vub,V,Vb,t_l(L)),
  244	    replace_in_cont(ContIn,Vub,V,ContOut)
  245	;   Status = unlimited(V,t_l(L)),
  246	    ContIn = ContOut
  247	).
  248dec_step_2_cont(t_Lu(L,U),CLP,l(V*K,OrdV),Class,applied,ContIn,ContOut) :-
  249	compare_d(CLP, >, 0, K),
  250	eval_d(CLP, U - L, Init),
  251	class_basis(Class,Deps),
  252	ub(Deps,CLP,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_),
  253	pivot_b(CLP,Vub,V,Vb,t_lu(L,U)),
  254	replace_in_cont(ContIn,Vub,V,ContOut).
  255dec_step_2_cont(t_none,_,l(V*_,_),_,unlimited(V,t_none),Cont,Cont).
  256
  257
  258
  259inc_step_2_cont(t_U(U),CLP,l(V*K,OrdV),Class,Status,ContIn,ContOut) :-
  260	compare_d(CLP, >, 0, K),
  261	(   lb(Class,CLP,OrdV,Vub-Vb-_)
  262	->  Status = applied,
  263	    pivot_a(CLP,Vub,V,Vb,t_u(U)),
  264	    replace_in_cont(ContIn,Vub,V,ContOut)
  265	;   Status = unlimited(V,t_u(U)),
  266	    ContIn = ContOut
  267	).
  268inc_step_2_cont(t_lU(L,U),CLP,l(V*K,OrdV),Class,applied,ContIn,ContOut) :-
  269	compare_d(CLP, >, 0, K),
  270	eval_d(CLP, L - U, Init),
  271	class_basis(Class,Deps),
  272	lb(Deps,CLP,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_),
  273	pivot_b(CLP,Vub,V,Vb,t_lu(L,U)),
  274	replace_in_cont(ContIn,Vub,V,ContOut).
  275inc_step_2_cont(t_L(L),CLP,l(V*K,OrdV),Class,Status,ContIn,ContOut) :-
  276	compare_d(CLP, <, 0, K),
  277	(   ub(Class,CLP,OrdV,Vub-Vb-_)
  278	->  Status = applied,
  279	    pivot_a(CLP,Vub,V,Vb,t_l(L)),
  280	    replace_in_cont(ContIn,Vub,V,ContOut)
  281	;   Status = unlimited(V,t_l(L)),
  282	    ContIn = ContOut
  283	).
  284inc_step_2_cont(t_Lu(L,U),CLP,l(V*K,OrdV),Class,applied,ContIn,ContOut) :-
  285	compare_d(CLP, <, 0, K),
  286	eval_d(CLP, U - L, Init),
  287	class_basis(Class,Deps),
  288	ub(Deps,CLP,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_),
  289	pivot_b(CLP,Vub,V,Vb,t_lu(L,U)),
  290	replace_in_cont(ContIn,Vub,V,ContOut).
  291inc_step_2_cont(t_none,_,l(V*_,_),_,unlimited(V,t_none),Cont,Cont).
  292
  293replace_in_cont([],_,_,[]).
  294replace_in_cont([H1|T1],X,Y,[H2|T2]) :-
  295	(   H1 == X
  296	->  H2 = Y,
  297	    T1 = T2
  298	;   H2 = H1,
  299	    replace_in_cont(T1,X,Y,T2)
  300	).
  301
  302dec_step([], _, optimum).
  303dec_step([l(V*K,OrdV)|Vs], CLP, Status) :-
  304	get_attr(V,clpcd_itf,Att),
  305	arg(2,Att,type(W)),
  306	arg(6,Att,class(Class)),
  307	(   dec_step_2(W, CLP, l(V*K,OrdV), Class, Status)
  308	->  true
  309	;   dec_step(Vs, CLP, Status)
  310	).
  311
  312dec_step_2(t_U(U),CLP,l(V*K,OrdV),Class,Status) :-
  313	compare_d(CLP, <, 0, K),
  314	(   lb(Class,CLP,OrdV,Vub-Vb-_)
  315	->  % found a lower bound
  316	    Status = applied,
  317	    pivot_a(CLP,Vub,V,Vb,t_u(U))
  318	;   Status = unlimited(V,t_u(U))
  319	).
  320dec_step_2(t_lU(L,U),CLP,l(V*K,OrdV),Class,applied) :-
  321	compare_d(CLP, <, 0, K),
  322	eval_d(CLP, L - U, Init),
  323	class_basis(Class,Deps),
  324	lb(Deps,CLP,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_),
  325	pivot_b(CLP,Vub,V,Vb,t_lu(L,U)).
  326dec_step_2(t_L(L),CLP,l(V*K,OrdV),Class,Status) :-
  327	compare_d(CLP, >, 0, K),
  328	(   ub(Class,CLP,OrdV,Vub-Vb-_)
  329	->  Status = applied,
  330	    pivot_a(CLP,Vub,V,Vb,t_l(L))
  331	;   Status = unlimited(V,t_l(L))
  332	).
  333dec_step_2(t_Lu(L,U),CLP,l(V*K,OrdV),Class,applied) :-
  334	compare_d(CLP, >, 0, K),
  335	eval_d(CLP, U - L, Init),
  336	class_basis(Class,Deps),
  337	ub(Deps,CLP,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_),
  338	pivot_b(CLP,Vub,V,Vb,t_lu(L,U)).
  339dec_step_2(t_none,_,l(V*_,_),_,unlimited(V,t_none)).
  340
  341inc_step([], _, optimum).	% if status has not been set yet: no changes
  342inc_step([l(V*K,OrdV)|Vs], CLP, Status) :-
  343	get_attr(V,clpcd_itf,Att),
  344	arg(2,Att,type(W)),
  345	arg(6,Att,class(Class)),
  346	(   inc_step_2(W, CLP, l(V*K,OrdV), Class, Status)
  347	->  true
  348	;   inc_step(Vs, CLP, Status)
  349	).
  350
  351inc_step_2(t_U(U),CLP,l(V*K,OrdV),Class,Status) :-
  352	compare_d(CLP, >, 0, K),
  353	(   lb(Class,CLP,OrdV,Vub-Vb-_)
  354	->  Status = applied,
  355	    pivot_a(CLP,Vub,V,Vb,t_u(U))
  356	;   Status = unlimited(V,t_u(U))
  357	).
  358inc_step_2(t_lU(L,U),CLP,l(V*K,OrdV),Class,applied) :-
  359	compare_d(CLP, >, 0, K),
  360	eval_d(CLP, L - U, Init),
  361	class_basis(Class,Deps),
  362	lb(Deps,CLP,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_),
  363	pivot_b(CLP,Vub,V,Vb,t_lu(L,U)).
  364inc_step_2(t_L(L),CLP,l(V*K,OrdV),Class,Status) :-
  365	compare_d(CLP, <, 0, K),
  366	(   ub(Class,CLP,OrdV,Vub-Vb-_)
  367	->  Status = applied,
  368	    pivot_a(CLP,Vub,V,Vb,t_l(L))
  369	;   Status = unlimited(V,t_l(L))
  370	).
  371inc_step_2(t_Lu(L,U),CLP,l(V*K,OrdV),Class,applied) :-
  372	compare_d(CLP, <, 0, K),
  373	eval_d(CLP, U - L, Init),
  374	class_basis(Class,Deps),
  375	ub(Deps,CLP,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_),
  376	pivot_b(CLP,Vub,V,Vb,t_lu(L,U)).
  377inc_step_2(t_none,_,l(V*_,_),_,unlimited(V,t_none)).
  378
  379% ------------------------- find the most constraining row --------------------
  380%
  381% The code for the lower and the upper bound are dual versions of each other.
  382% The only difference is in the orientation of the comparisons.
  383% Indeps are ruled out by their types.
  384% If there is no bound, this fails.
  385%
  386% *** The actual lb and ub on an indep variable X are [lu]b + b(X), where b(X)
  387% is the value of the active bound.
  388%
  389% Nota bene: We must NOT consider infeasible rows as candidates to
  390%	     leave the basis!
  391%
  392% ub(Class,CLP,OrdX,Ub)
  393%
  394% See lb/4: this is similar
  395
  396ub(Class,CLP,OrdX,Ub) :-
  397	class_basis(Class,Deps),
  398	ub_first(Deps,CLP,OrdX,Ub).
  399
  400% ub_first(Deps,X,Dep-W-Ub)
  401%
  402% Finds the tightest upperbound for variable X from the linear equations of
  403% basis variables Deps, and puts the resulting bound in Ub. Dep is the basis
  404% variable that generates the bound, and W is bound of that variable that has
  405% to be activated to achieve this.
  406
  407ub_first([Dep|Deps],CLP,OrdX,Tightest) :-
  408	(   get_attr(Dep,clpcd_itf,Att),
  409	    arg(2,Att,type(Type)),
  410	    arg(4,Att,lin(Lin)),
  411	    ub_inner(Type,CLP,OrdX,Lin,W,Ub)
  412	->  ub(Deps,CLP,OrdX,Dep-W-Ub,Tightest)
  413	;   ub_first(Deps,CLP,OrdX,Tightest)
  414	).
  415
  416% ub(Deps,OrdX,TightestIn,TightestOut)
  417%
  418% See lb/5: this is similar
  419
  420ub([],_,_,T0,T0).
  421ub([Dep|Deps],CLP,OrdX,T0,T1) :-
  422	(   get_attr(Dep,clpcd_itf,Att),
  423	    arg(2,Att,type(Type)),
  424	    arg(4,Att,lin(Lin)),
  425	    ub_inner(Type,CLP,OrdX,Lin,W,Ub),
  426            T0 = _-Ubb,
  427            compare_d(CLP, <, Ub, Ubb)
  428	->  ub(Deps,CLP,OrdX,Dep-W-Ub,T1)	% tighter bound, use new bound
  429	;   ub(Deps,CLP,OrdX,T0,T1)	% no tighter bound, keep current one
  430	).
  431
  432% ub_inner(Type,CLP,OrdX,Lin,W,Ub)
  433%
  434% See lb_inner/6: this is similar
  435
  436ub_inner(t_l(L),CLP,OrdX,Lin,t_L(L),Ub) :-
  437	nf_rhs_x(CLP,Lin,OrdX,Rhs,K),
  438	% Rhs is right hand side of lin. eq. Lin containing term X*K
  439	compare_d(CLP, >, 0, K),
  440        compare_d(CLP, =<, L, Rhs),
  441        div_d(CLP, L - Rhs, K, Ub).
  442ub_inner(t_u(U),CLP,OrdX,Lin,t_U(U),Ub) :-
  443	nf_rhs_x(CLP,Lin,OrdX,Rhs,K),
  444	compare_d(CLP, <, 0, K),
  445        compare_d(CLP, =<, Rhs, U),
  446        div_d(CLP, U - Rhs, K, Ub).
  447ub_inner(t_lu(L,U),CLP,OrdX,Lin,W,Ub) :-
  448	nf_rhs_x(CLP,Lin,OrdX,Rhs,K),
  449	(   compare_d(CLP, >, 0, K) % use lowerbound
  450	->  W = t_Lu(L,U),
  451            compare_d(CLP, =<, L, Rhs),
  452	    div_d(CLP, L - Rhs, K, Ub)
  453	;   compare_d(CLP, <, 0, K) % use upperbound
  454	->  W = t_lU(L,U),
  455            compare_d(CLP, =<, Rhs, U),
  456	    div_d(CLP, U - Rhs, K, Ub)
  457	).
  458
  459
  460% lb(Class,OrdX,Lb)
  461%
  462% Returns in Lb how much we can lower the upperbound of X without violating
  463% a bound of the basisvariables.
  464% Lb has the form Dep-W-Lb with Dep the variable whose bound is violated when
  465% lowering the bound for X more, W the actual bound that has to be activated
  466% and Lb the amount that the upperbound can be lowered.
  467% X has ordering OrdX and class Class.
  468
  469lb(Class,CLP,OrdX,Lb) :-
  470	class_basis(Class,Deps),
  471	lb_first(Deps,CLP,OrdX,Lb).
  472
  473% lb_first(Deps,CLP,OrdX,Tightest)
  474%
  475% Returns in Tightest how much we can lower the upperbound of X without
  476% violating a bound of Deps.
  477% Tightest has the form Dep-W-Lb with Dep the variable whose bound is violated
  478% when lowering the bound for X more, W the actual bound that has to be
  479% activated and Lb the amount that the upperbound can be lowered. X has
  480% ordering attribute OrdX.
  481
  482lb_first([Dep|Deps],CLP,OrdX,Tightest) :-
  483	(   get_attr(Dep,clpcd_itf,Att),
  484	    arg(2,Att,type(Type)),
  485	    arg(4,Att,lin(Lin)),
  486	    lb_inner(Type, CLP, OrdX, Lin, W, Lb)
  487	->  lb(Deps,CLP,OrdX,Dep-W-Lb,Tightest)
  488	;   lb_first(Deps, CLP, OrdX, Tightest)
  489	).
  490
  491% lb(Deps,CLP,OrdX,TightestIn,TightestOut)
  492%
  493% See lb_first/3: this one does the same thing, but is used for the steps after
  494% the first one and remembers the tightest bound so far.
  495
  496lb([],_,_,T0,T0).
  497lb([Dep|Deps],CLP,OrdX,T0,T1) :-
  498	(   get_attr(Dep,clpcd_itf,Att),
  499	    arg(2,Att,type(Type)),
  500	    arg(4,Att,lin(Lin)),
  501	    lb_inner(Type, CLP, OrdX, Lin, W, Lb),
  502	    T0 = _-Lbb,
  503	    compare_d(CLP, >, Lb, Lbb)  % Lb > Lbb: choose the least lowering, others
  504					% might violate bounds
  505	->  lb(Deps,CLP,OrdX,Dep-W-Lb,T1)
  506	;   lb(Deps,CLP,OrdX,T0,T1)
  507	).
  508
  509% lb_inner(Type,CLP,X,Lin,W,Lb)
  510%
  511% Returns in Lb how much lower we can make X without violating a bound
  512% by using the linear equation Lin of basis variable B which has type
  513% Type and which has to activate a bound (type W) to do so.
  514%
  515% E.g. when B has a lowerbound L, then L should always be smaller than I + R.
  516% So a lowerbound of X (which has scalar K in Lin), could be at most
  517% (L-(I+R))/K lower than its upperbound (if K is positive).
  518% Also note that Lb should always be smaller than 0, otherwise the row is
  519% not feasible.
  520% X has ordering attribute OrdX.
  521
  522lb_inner(t_l(L),CLP,OrdX,Lin,t_L(L),Lb) :-
  523	nf_rhs_x(CLP,Lin,OrdX,Rhs,K), % if linear equation Lin contains the term
  524				  % X*K, Rhs is the right hand side of that
  525				  % equation
  526        compare_d(CLP, <, 0, K),
  527        compare_d(CLP, =<, L, Rhs),
  528	div_d(CLP, L - Rhs, K, Lb).
  529lb_inner(t_u(U),CLP,OrdX,Lin,t_U(U),Lb) :-
  530	nf_rhs_x(CLP,Lin,OrdX,Rhs,K),
  531        compare_d(CLP, >, 0, K),
  532        compare_d(CLP, =<, Rhs, U),
  533	div_d(CLP, U - Rhs, K, Lb).
  534lb_inner(t_lu(L,U),CLP,OrdX,Lin,W,Lb) :-
  535	nf_rhs_x(CLP,Lin,OrdX,Rhs,K),
  536	(   compare_d(CLP, >, 0, K)
  537	->  compare_d(CLP, =<, Rhs, U),
  538            W = t_lU(L,U),
  539            div_d(CLP, U - Rhs, K, Lb)
  540	;   compare_d(CLP, <, 0, K)
  541	->  compare_d(CLP, =<, L, Rhs),
  542	    W = t_Lu(L,U),
  543            div_d(CLP, L - Rhs, K, Lb)
  544	).
  545
  546% ---------------------------------- equations --------------------------------
  547%
  548% backsubstitution will not make the system infeasible, if the bounds on the
  549% indep vars are obeyed, but some implied values might pop up in rows where X
  550% occurs
  551%	-) special case X=Y during bs -> get rid of dependend var(s), alias
  552%
  553
  554solve(CLP,Lin) :-
  555	Lin = [I,_|H],
  556	solve(H,CLP,Lin,I,Bindings,[]),
  557	export_binding(Bindings).
  558
  559% solve(Hom,Lin,I,Bind,BindT)
  560%
  561% Solves a linear equation Lin = [I,_|H] = 0 and exports the generated bindings
  562
  563solve([],CLP,_,I,Bind0,Bind0) :-
  564	!,
  565	compare_d(CLP, =, 0, I). % redundant or trivially unsat
  566solve(H,CLP,Lin,_,Bind0,BindT) :-
  567	sd(H,CLP,[],ClassesUniq,9-9-0,Category-Selected-_,NV,NVT),
  568	get_attr(Selected,clpcd_itf,Att),
  569	arg(5,Att,order(Ord)),
  570	isolate(CLP, Ord, Lin, Lin1),	% Lin = 0 => Selected = Lin1
  571	(   Category = 1 % classless variable, no bounds
  572	->  setarg(4,Att,lin(Lin1)),
  573	    Lin1 = [Inhom,_|Hom],
  574	    bs_collect_binding(Hom,Selected,Inhom,Bind0,BindT),
  575	    eq_classes(CLP, NV, NVT, ClassesUniq)
  576	;   Category = 2 % class variable, no bounds
  577	->  arg(6,Att,class(NewC)),
  578	    class_allvars(NewC,Deps),
  579	    (   ClassesUniq = [_] % rank increasing
  580	    ->	bs_collect_bindings(Deps, CLP, Ord, Lin1, Bind0, BindT)
  581	    ;   Bind0 = BindT,
  582		bs(Deps, CLP, Ord, Lin1)
  583	    ),
  584	    eq_classes(CLP, NV, NVT, ClassesUniq)
  585	;   Category = 3 % classless variable, all variables in Lin and
  586			 % Selected are bounded
  587	->  arg(2,Att,type(Type)),
  588	    setarg(4,Att,lin(Lin1)),
  589	    deactivate_bound(Type,Selected),
  590	    eq_classes(CLP, NV, NVT, ClassesUniq),
  591	    basis_add(Selected,Basis),
  592	    undet_active(CLP,Lin1),	% we can't tell which bound will likely be a
  593				% problem at this point
  594	    Lin1 = [Inhom,_|Hom],
  595	    bs_collect_binding(Hom,Selected,Inhom,Bind0,Bind1),	% only if
  596								% Hom = []
  597	    rcbl(Basis,CLP,Bind1,BindT) % reconsider entire basis
  598	;   Category = 4 % class variable, all variables in Lin and Selected
  599			 % are bounded
  600	->  arg(2,Att,type(Type)),
  601	    arg(6,Att,class(NewC)),
  602	    class_allvars(NewC,Deps),
  603	    (   ClassesUniq = [_] % rank increasing
  604	    ->	bs_collect_bindings(Deps, CLP, Ord, Lin1, Bind0, Bind1)
  605	    ;   Bind0 = Bind1,
  606		bs(Deps, CLP, Ord, Lin1)
  607	    ),
  608	    deactivate_bound(Type,Selected),
  609	    basis_add(Selected,Basis),
  610	    % eq_classes( NV, NVT, ClassesUniq),
  611	    %  4 -> var(NV)
  612	    equate(ClassesUniq,_),
  613	    undet_active(CLP,Lin1),
  614	    rcbl(Basis,CLP,Bind1,BindT)
  615	).
  616
  617% solve_ord_x(CLP,Lin,OrdX,ClassX)
  618%
  619% Much like solve, but we solve for a particular variable of type t_none, it has
  620% the ordering of X and its class as input. This also means that X has a class
  621% which is not sure in solve_x/2.
  622
  623solve_ord_x(CLP,Lin,OrdX,ClassX) :-
  624	Lin = [I,_|H],
  625	solve_ord_x(H,CLP,Lin,I,OrdX,ClassX,Bindings,[]),
  626	export_binding(Bindings).
  627
  628solve_ord_x([],_,_,I,_,_,Bind0,Bind0) :-
  629	I =:= 0.
  630solve_ord_x([_|_],CLP,Lin,_,OrdX,ClassX,Bind0,BindT) :-
  631	isolate(CLP, OrdX, Lin, Lin1),
  632	Lin1 = [_,_|H1],
  633	sd(H1,CLP,[],ClassesUniq1,9-9-0,_,NV,NVT), % do sd on Lin without X, then
  634					       % add class of X
  635	ord_add_element(ClassesUniq1,ClassX,ClassesUniq),
  636	class_allvars(ClassX,Deps),
  637	(   ClassesUniq = [_] % rank increasing
  638	->  bs_collect_bindings(Deps, CLP, OrdX, Lin1, Bind0, BindT)
  639	;   Bind0 = BindT,
  640	    bs(Deps, CLP, OrdX, Lin1)
  641	),
  642	eq_classes(CLP,NV,NVT,ClassesUniq).
  643
  644% sd(H,[],ClassesUniq,9-9-0,Category-Selected-_,NV,NVT)
  645
  646% sd(Hom,ClassesIn,ClassesOut,PreferenceIn,PreferenceOut,[NV|NVTail],NVTail)
  647%
  648% ClassesOut is a sorted list of the different classes that are either in
  649% ClassesIn or that are the classes of the variables in Hom. Variables that do
  650% not belong to a class yet, are put in the difference list NV.
  651
  652sd([],_,Class0,Class0,Preference0,Preference0,NV0,NV0).
  653sd([l(X*K,_)|Xs],CLP,Class0,ClassN,Preference0,PreferenceN,NV0,NVt) :-
  654	get_attr(X,clpcd_itf,Att),
  655	(   arg(6,Att,class(Xc)) % old: has class
  656	->  NV0 = NV1,
  657	    ord_add_element(Class0,Xc,Class1),
  658	    (   arg(2,Att,type(t_none))
  659	    ->  preference(CLP,Preference0,2-X-K,Preference1)
  660		    % has class, no bounds => category 2
  661	    ;   preference(CLP,Preference0,4-X-K,Preference1)
  662		    % has class, is bounded => category 4
  663	    )
  664	;   % new: has no class
  665	    Class1 = Class0,
  666	    NV0 = [X|NV1], % X has no class yet, add to list of new variables
  667	    (   arg(2,Att,type(t_none))
  668	    ->  preference(CLP,Preference0,1-X-K,Preference1)
  669		    % no class, no bounds => category 1
  670	    ;   preference(CLP,Preference0,3-X-K,Preference1)
  671		    % no class, is bounded => category 3
  672	    )
  673	),
  674	sd(Xs,CLP,Class1,ClassN,Preference1,PreferenceN,NV1,NVt).
  675
  676%
  677% A is best sofar, B is current
  678% smallest preferred
  679preference(CLP,A,B,Pref) :-
  680	A = Px-_-Ka,
  681	B = Py-_-Kb,
  682	(   Px < Py
  683	->  Pref = A
  684	;   Px > Py
  685        ->  Pref = B
  686        ;   compare_d(CLP, >, abs(Ka), abs(Kb))
  687        ->  Pref = A
  688        ;   Pref = B
  689	).
  690
  691% eq_classes(CLP,NV,NVTail,Cs)
  692%
  693% Attaches all classless variables NV to a new class and equates all other
  694% classes with this class. The equate operation only happens after attach_class
  695% because the unification of classes can bind the tail of the AllVars attribute
  696% to a nonvar and then the attach_class operation wouldn't work.
  697
  698eq_classes(_,NV,_,Cs) :-
  699	var(NV),
  700	!,
  701	equate(Cs,_).
  702eq_classes(CLP,NV,NVT,Cs) :-
  703	class_new(Su,CLP,NV,NVT,[]), % make a new class Su with NV as the variables
  704	attach_class(NV,Su), % attach the variables NV to Su
  705	equate(Cs,Su).
  706
  707equate([],_).
  708equate([X|Xs],X) :- equate(Xs,X).
  709
  710%
  711% assert: none of the Vars has a class attribute yet
  712%
  713attach_class(Xs,_) :-
  714	var(Xs), % Tail
  715	!.
  716attach_class([X|Xs],Class) :-
  717	get_attr(X,clpcd_itf,Att),
  718	setarg(6,Att,class(Class)),
  719	attach_class(Xs,Class).
  720
  721% unconstrained(CLP,Lin,Uc,Kuc,Rest)
  722%
  723% Finds an unconstrained variable Uc (type(t_none)) in Lin with scalar Kuc and
  724% removes it from Lin to return Rest.
  725
  726unconstrained(CLP,Lin,Uc,Kuc,Rest) :-
  727	Lin = [_,_|H],
  728	sd(H,CLP,[],_,9-9-0,Category-Uc-_,_,_),
  729	Category =< 2,
  730	get_attr(Uc,clpcd_itf,Att),
  731	arg(5,Att,order(OrdUc)),
  732	delete_factor(OrdUc,Lin,Rest,Kuc).
  733
  734%
  735% point the vars in Lin into the same equivalence class
  736% maybe join some global data
  737%
  738same_class([],_).
  739same_class([l(X*_,_)|Xs],Class) :-
  740	get_or_add_class(X,Class),
  741	same_class(Xs,Class).
  742
  743% deactivate_bound(Type,Variable)
  744%
  745% The Type of the variable is changed to reflect the deactivation of its
  746% bounds.
  747% t_L(_) becomes t_l(_), t_lU(_,_) becomes t_lu(_,_) and so on.
  748
  749deactivate_bound(t_l(_),_).
  750deactivate_bound(t_u(_),_).
  751deactivate_bound(t_lu(_,_),_).
  752deactivate_bound(t_L(L),X) :-
  753	get_attr(X,clpcd_itf,Att),
  754	setarg(2,Att,type(t_l(L))).
  755deactivate_bound(t_Lu(L,U),X) :-
  756	get_attr(X,clpcd_itf,Att),
  757	setarg(2,Att,type(t_lu(L,U))).
  758deactivate_bound(t_U(U),X) :-
  759	get_attr(X,clpcd_itf,Att),
  760	setarg(2,Att,type(t_u(U))).
  761deactivate_bound(t_lU(L,U),X) :-
  762	get_attr(X,clpcd_itf,Att),
  763	setarg(2,Att,type(t_lu(L,U))).
  764
  765% intro_at(CLP,X,Value,Type)
  766%
  767% Variable X gets new type Type which reflects the activation of a bound with
  768% value Value. In the linear equations of all the variables belonging to the
  769% same class as X, X is substituted by [0,Value,X] to reflect the new active
  770% bound.
  771
  772intro_at(CLP,X,Value,Type) :-
  773	get_attr(X,clpcd_itf,Att),
  774	arg(5,Att,order(Ord)),
  775	arg(6,Att,class(Class)),
  776	setarg(2,Att,type(Type)),
  777	(   compare_d(CLP, =, 0, Value)
  778	->  true
  779	;   backsubst_delta(CLP, Class, Ord, X, Value)
  780	).
  781
  782% undet_active(CLP, Lin)
  783%
  784% For each variable in the homogene part of Lin, a bound is activated
  785% if an inactive bound exists. (t_l(L) becomes t_L(L) and so on)
  786
  787undet_active(CLP, [_,_|H]) :-
  788	undet_active_h(H, CLP).
  789
  790% undet_active_h(Hom)
  791%
  792% For each variable in homogene part Hom, a bound is activated if an
  793% inactive bound exists (t_l(L) becomes t_L(L) and so on)
  794
  795undet_active_h([],_).
  796undet_active_h([l(X*_,_)|Xs],CLP) :-
  797	get_attr(X,clpcd_itf,Att),
  798	arg(2,Att,type(Type)),
  799	undet_active(Type,CLP,X),
  800	undet_active_h(Xs,CLP).
  801
  802% undet_active(Type,CLP,Var)
  803%
  804% An inactive bound of Var is activated if such exists
  805% t_lu(L,U) is arbitrarily chosen to become t_Lu(L,U)
  806
  807undet_active(t_none,_,_).	% type_activity
  808undet_active(t_L(_),_,_).
  809undet_active(t_Lu(_,_),_,_).
  810undet_active(t_U(_),_,_).
  811undet_active(t_lU(_,_),_,_).
  812undet_active(t_l(L),CLP,X) :- intro_at(CLP,X,L,t_L(L)).
  813undet_active(t_u(U),CLP,X) :- intro_at(CLP,X,U,t_U(U)).
  814undet_active(t_lu(L,U),CLP,X) :- intro_at(CLP,X,L,t_Lu(L,U)).
  815
  816% ----------------------------- manipulate the basis --------------------------
  817
  818% basis_add(X,NewBasis)
  819%
  820% NewBasis is the result of adding X to the basis of the class to which X
  821% belongs.
  822
  823basis_add(X,NewBasis) :-
  824	get_attr(X,clpcd_itf,Att),
  825	arg(6,Att,class(Cv)),
  826	class_basis_add(Cv,X,NewBasis).
  827
  828% basis_pivot(Leave,Enter)
  829%
  830% Removes Leave from the basis of the class to which it belongs, and adds
  831% Enter to that basis.
  832
  833basis_pivot(Leave,Enter) :-
  834	get_attr(Leave,clpcd_itf,Att),
  835	arg(6,Att,class(Cv)),
  836	class_basis_pivot(Cv,Enter,Leave).
  837
  838% ----------------------------------- pivot -----------------------------------
  839
  840% pivot_a(Dep,Indep,IndepT,DepT)
  841%
  842% Removes Dep from the basis, puts Indep in, and pivots the equation of
  843% Dep to become one of Indep. The type of Dep becomes DepT (which means
  844% it gets deactivated), the type of Indep becomes IndepT (which means it
  845% gets activated)
  846
  847
  848pivot_a(CLP,Dep,Indep,Vb,Wd) :-
  849	basis_pivot(Dep,Indep),
  850	get_attr(Indep,clpcd_itf,Att),
  851	arg(2,Att,type(Type)),
  852	arg(5,Att,order(Ord)),
  853	arg(6,Att,class(Class)),
  854	pivot(CLP,Dep,Class,Ord,Vb,Type),
  855	get_attr(Indep,clpcd_itf,Att2), %changed?
  856	setarg(2,Att2,type(Wd)).
  857
  858pivot_b(CLP,Vub,V,Vb,Wd) :-
  859	(   Vub == V
  860	->  get_attr(V,clpcd_itf,Att),
  861	    arg(5,Att,order(Ord)),
  862	    arg(6,Att,class(Class)),
  863	    setarg(2,Att,type(Vb)),
  864	    pivot_b_delta(Vb, CLP, Delta), % nonzero(Delta)
  865	    backsubst_delta(CLP, Class, Ord, V, Delta)
  866	;   pivot_a(CLP,Vub,V,Vb,Wd)
  867	).
  868
  869pivot_b_delta(t_Lu(L,U), CLP, Delta) :- eval_d(CLP, L-U, Delta).
  870pivot_b_delta(t_lU(L,U), CLP, Delta) :- eval_d(CLP, U-L, Delta).
  871
  872% select_active_bound(Type,Bound)
  873%
  874% Returns the bound that is active in Type (if such exists, 0 otherwise)
  875
  876select_active_bound(t_L(L),L).
  877select_active_bound(t_Lu(L,_),L).
  878select_active_bound(t_U(U),U).
  879select_active_bound(t_lU(_,U),U).
  880select_active_bound(t_none,0).
  881%
  882% for project.pl
  883%
  884select_active_bound(t_l(_),0).
  885select_active_bound(t_u(_),0).
  886select_active_bound(t_lu(_,_),0).
  887
  888% pivot(Dep,Class,IndepOrd,DepAct,IndAct)
  889%
  890% See pivot/2.
  891% In addition, variable Indep with ordering IndepOrd has an active bound IndAct
  892
  893%
  894%
  895% Pivot taking care of rhs and active states
  896%
  897pivot(CLP,Dep,Class,IndepOrd,DepAct,IndAct) :-
  898	get_attr(Dep,clpcd_itf,Att),
  899	arg(4,Att,lin(H)),
  900	arg(5,Att,order(DepOrd)),
  901	setarg(2,Att,type(DepAct)),
  902	select_active_bound(DepAct,AbvD), % New current value for Dep
  903	select_active_bound(IndAct,AbvI), % Old current value of Indep
  904	delete_factor(IndepOrd,H,H0,Coeff), % Dep = ... + Coeff*Indep + ...
  905	eval_d(CLP, -AbvD, AbvDm),
  906	eval_d(CLP, -AbvI, AbvIm),
  907	add_linear_f1(CLP, [0,AbvIm], Coeff, H0, H1),
  908	div_d(CLP, -1, Coeff, K),
  909	add_linear_ff(CLP, H1, K, [0,AbvDm,l(Dep* -1,DepOrd)], K, H2),
  910	    % Indep = -1/Coeff*... + 1/Coeff*Dep
  911	add_linear_11(CLP, H2, [0,AbvIm], Lin),
  912	backsubst(CLP, Class, IndepOrd, Lin).
  913
  914% backsubst_delta(CLP,Class,OrdX,X,Delta)
  915%
  916% X with ordering attribute OrdX, is substituted in all linear equations of
  917% variables in the class Class, by linear equation [0,Delta,l(X*1,OrdX)]. This
  918% reflects the activation of a bound.
  919
  920backsubst_delta(CLP, Class, OrdX, X, Delta) :-
  921	backsubst(CLP, Class, OrdX, [0,Delta,l(X*1,OrdX)]).
  922
  923% backsubst(Class,OrdX,Lin)
  924%
  925% X with ordering OrdX is substituted in all linear equations of variables in
  926% the class Class, by linear equation Lin
  927
  928backsubst(CLP, Class, OrdX, Lin) :-
  929	class_allvars(Class,Allvars),
  930	bs(Allvars, CLP, OrdX, Lin).
  931
  932% bs(Vars,OrdV,Lin)
  933%
  934% In all linear equations of the variables Vars, variable V with ordering
  935% attribute OrdV is substituted by linear equation Lin.
  936%
  937% valid if nothing will go ground
  938%
  939
  940bs(Xs, _, _, _) :-
  941	var(Xs),
  942	!.
  943bs([X|Xs], CLP, OrdV, Lin) :-
  944	(   get_attr(X,clpcd_itf,Att),
  945	    arg(4,Att,lin(LinX)),
  946	    nf_substitute(CLP, OrdV, Lin, LinX, LinX1) % does not change attributes
  947	->  setarg(4,Att,lin(LinX1)),
  948	    bs(Xs, CLP, OrdV, Lin)
  949	;   bs(Xs, CLP, OrdV, Lin)
  950	).
  951
  952%
  953% rank increasing backsubstitution
  954%
  955
  956% bs_collect_bindings(Deps,SelectedOrd,Lin,Bind,BindT)
  957%
  958% Collects bindings (of the form [X-I] where X = I is the binding) by
  959% substituting Selected in all linear equations of the variables Deps (which
  960% are of the same class), by Lin. Selected has ordering attribute SelectedOrd.
  961%
  962% E.g. when V = 2X + 3Y + 4, X = 3V + 2Z and Y = 4X + 3
  963% we can substitute V in the linear equation of X: X = 6X + 9Y + 2Z + 12
  964% we can't substitute V in the linear equation of Y of course.
  965
  966bs_collect_bindings(Xs, _, _, _, Bind0, BindT) :-
  967	var(Xs),
  968	!,
  969	Bind0 = BindT.
  970bs_collect_bindings([X|Xs], CLP, OrdV, Lin, Bind0, BindT) :-
  971	(   get_attr(X,clpcd_itf,Att),
  972	    arg(4,Att,lin(LinX)),
  973	    nf_substitute(CLP, OrdV, Lin, LinX, LinX1) % does not change attributes
  974	->  setarg(4,Att,lin(LinX1)),
  975	    LinX1 = [Inhom,_|Hom],
  976	    bs_collect_binding(Hom,X,Inhom,Bind0,Bind1),
  977	    bs_collect_bindings(Xs, CLP, OrdV, Lin, Bind1, BindT)
  978	;   bs_collect_bindings(Xs, CLP, OrdV, Lin, Bind0, BindT)
  979	).
  980
  981% bs_collect_binding(Hom,Selected,Inhom,Bind,BindT)
  982%
  983% Collects binding following from Selected = Hom + Inhom.
  984% If Hom = [], returns the binding Selected-Inhom (=0)
  985%
  986bs_collect_binding([],X,Inhom) --> [X-Inhom].
  987bs_collect_binding([_|_],_,_) --> [].
  988
  989%
  990% reconsider the basis
  991%
  992
  993% rcbl(Basis,Bind,BindT)
  994%
  995%
  996
  997rcbl([],_,Bind0,Bind0).
  998rcbl([X|Continuation],CLP,Bind0,BindT) :-
  999	(   rcb_cont(CLP,X,Status,Violated,Continuation,NewContinuation) % have a culprit
 1000	->  rcbl_status(Status,CLP,X,NewContinuation,Bind0,BindT,Violated)
 1001	;   rcbl(Continuation,CLP,Bind0,BindT)
 1002	).
 1003
 1004rcb_cont(CLP,X,Status,Violated,ContIn,ContOut) :-
 1005	get_attr(X,clpcd_itf,Att),
 1006	arg(2,Att,type(Type)),
 1007	arg(4,Att,lin([I,R|H])),
 1008	(   Type = t_l(L) % case 1: lowerbound: R + I should always be larger
 1009			  % than the lowerbound
 1010	->  compare_d(CLP, =<, R + I, L),
 1011	    Violated = l(L),
 1012	    inc_step_cont(H, CLP, Status, ContIn, ContOut)
 1013	;   Type = t_u(U) % case 2: upperbound: R + I should always be smaller
 1014			  % than the upperbound
 1015	->  compare_d(CLP, >, R + I, U),
 1016	    Violated = u(U),
 1017	    dec_step_cont(H,CLP,Status,ContIn,ContOut)
 1018	;   Type = t_lu(L,U) % case 3: check both
 1019	->  eval_d(CLP, R + I, At),
 1020	    (   compare_d(CLP, =<, At, L)
 1021	    ->	Violated = l(L),
 1022		inc_step_cont(H, CLP, Status, ContIn, ContOut)
 1023	    ;   compare_d(CLP, >=, At, U)
 1024	    ->	Violated = u(U),
 1025		dec_step_cont(H,CLP,Status,ContIn,ContOut)
 1026	    )
 1027	). % other types imply nonbasic variable or unbounded variable
 1028
 1029%
 1030% reconsider one element of the basis
 1031% later: lift the binds
 1032%
 1033reconsider(CLP,X) :-
 1034	rcb(CLP,X,Status,Violated),
 1035	!,
 1036	rcbl_status(Status,CLP,X,[],Binds,[],Violated),
 1037	export_binding(Binds).
 1038reconsider(_,_).
 1039
 1040%
 1041% Find a basis variable out of its bound or at its bound
 1042% Try to move it into whithin its bound
 1043%   a) impossible -> fail
 1044%   b) optimum at the bound -> implied value
 1045%   c) else look at the remaining basis variables
 1046%
 1047%
 1048% Idea: consider a variable V with linear equation Lin.
 1049% When a bound on a variable X of Lin gets activated, its value, multiplied
 1050% with the scalar of X, is added to the R component of Lin.
 1051% When we consider the lowerbound of V, it must be smaller than R + I, since R
 1052% contains at best the lowerbounds of the variables in Lin (but could contain
 1053% upperbounds, which are of course larger). So checking this can show the
 1054% violation of a bound of V. A similar case works for the upperbound.
 1055
 1056rcb(CLP,X,Status,Violated) :-
 1057	get_attr(X,clpcd_itf,Att),
 1058	arg(2,Att,type(Type)),
 1059	arg(4,Att,lin([I,R|H])),
 1060	(   Type = t_l(L) % case 1: lowerbound: R + I should always be larger
 1061			  % than the lowerbound
 1062	->  compare_d(CLP, =<, R + I, L),
 1063	    Violated = l(L),
 1064	    inc_step(H, CLP, Status)
 1065	;   Type = t_u(U) % case 2: upperbound: R + I should always be smaller
 1066			  % than the upperbound
 1067	->  compare_d(CLP, >=, R + I, U),
 1068	    Violated = u(U),
 1069	    dec_step(H, CLP, Status)
 1070	;   Type = t_lu(L,U) % case 3: check both
 1071	->  eval_d(CLP, R + I, At),
 1072	    (   compare_d(CLP, =<, At, L)
 1073	    ->	Violated = l(L),
 1074		inc_step(H, CLP, Status)
 1075	    ;   compare_d(CLP, >=, At, U)
 1076	    ->	Violated = u(U),
 1077		dec_step(H, CLP, Status)
 1078	    )
 1079	). % other types imply nonbasic variable or unbounded variable
 1080
 1081% rcbl_status(Status,CLP,X,Continuation,[Bind|BindT],BindT,Violated)
 1082%
 1083%
 1084
 1085rcbl_status(optimum,CLP,X,Cont,B0,Bt,Violated) :- rcbl_opt(Violated,CLP,X,Cont,B0,Bt).
 1086rcbl_status(applied,CLP,X,Cont,B0,Bt,Violated) :- rcbl_app(Violated,CLP,X,Cont,B0,Bt).
 1087rcbl_status(unlimited(Indep,DepT),CLP,X,Cont,B0,Bt,Violated) :-
 1088	rcbl_unl(Violated,CLP,X,Cont,B0,Bt,Indep,DepT).
 1089
 1090%
 1091% Might reach optimum immediately without changing the basis,
 1092% but in general we must assume that there were pivots.
 1093% If the optimum meets the bound, we backsubstitute the implied
 1094% value, solve will call us again to check for further implied
 1095% values or unsatisfiability in the rank increased system.
 1096%
 1097rcbl_opt(l(L),CLP,X,Continuation,B0,B1) :-
 1098	get_attr(X,clpcd_itf,Att),
 1099	arg(2,Att,type(Type)),
 1100	arg(3,Att,strictness(Strict)),
 1101	arg(4,Att,lin(Lin)),
 1102	Lin = [I,R|_],
 1103	eval_d(CLP, R + I, Opt),
 1104	(   compare_d(CLP, <, L, Opt)
 1105	->  narrow_u(Type,X,Opt), % { X =< Opt }
 1106	    rcbl(Continuation,CLP,B0,B1)
 1107	;   compare_d(CLP, =, L, Opt),
 1108	    Strict /\ 2 =:= 0, % meets lower
 1109	    eval_d(CLP, -Opt, Mop),
 1110	    normalize_scalar(Mop,MopN),
 1111	    add_linear_11(CLP, MopN, Lin, Lin1),
 1112	    Lin1 = [Inhom,_|Hom],
 1113	    (   Hom = []
 1114	    ->  rcbl(Continuation,CLP,B0,B1) % would not callback
 1115	    ;   solve(Hom,CLP,Lin1,Inhom,B0,B1)
 1116	    )
 1117	).
 1118rcbl_opt(u(U),CLP,X,Continuation,B0,B1) :-
 1119	get_attr(X,clpcd_itf,Att),
 1120	arg(2,Att,type(Type)),
 1121	arg(3,Att,strictness(Strict)),
 1122	arg(4,Att,lin(Lin)),
 1123	Lin = [I,R|_],
 1124	eval_d(CLP, R + I, Opt),
 1125	(   compare_d(CLP, >, U, Opt)
 1126	->  narrow_l(Type,X,Opt), % { X >= Opt }
 1127	    rcbl(Continuation,CLP,B0,B1)
 1128	;   compare_d(CLP, =, U, Opt),
 1129	    Strict /\ 1 =:= 0, % meets upper
 1130	    eval_d(CLP, -Opt, Mop),
 1131	    normalize_scalar(Mop,MopN),
 1132	    add_linear_11(CLP, MopN, Lin, Lin1),
 1133	    Lin1 = [Inhom,_|Hom],
 1134	    (   Hom = []
 1135	    ->  rcbl(Continuation,CLP,B0,B1) % would not callback
 1136	    ;   solve(Hom,CLP,Lin1,Inhom,B0,B1)
 1137	    )
 1138	).
 1139
 1140%
 1141% Basis has already changed when this is called
 1142%
 1143rcbl_app(l(L),CLP,X,Continuation,B0,B1) :-
 1144	get_attr(X,clpcd_itf,Att),
 1145	arg(4,Att,lin([I,R|H])),
 1146	(   compare_d(CLP, >, R + I, L) % within bound now
 1147	->  rcbl(Continuation,CLP,B0,B1)
 1148	;   inc_step(H, CLP, Status),
 1149	    rcbl_status(Status,CLP,X,Continuation,B0,B1,l(L))
 1150	).
 1151rcbl_app(u(U),CLP,X,Continuation,B0,B1) :-
 1152	get_attr(X,clpcd_itf,Att),
 1153	arg(4,Att,lin([I,R|H])),
 1154	(   compare_d(CLP, <, R + I, U) % within bound now
 1155	->  rcbl(Continuation,CLP,B0,B1)
 1156	;   dec_step(H, CLP, Status),
 1157	    rcbl_status(Status,CLP,X,Continuation,B0,B1,u(U))
 1158	).
 1159%
 1160% This is never called for a t_lu culprit
 1161%
 1162rcbl_unl(l(L),CLP,X,Continuation,B0,B1,Indep,DepT) :-
 1163	pivot_a(CLP,X,Indep,t_L(L),DepT), % changes the basis
 1164	rcbl(Continuation,CLP,B0,B1).
 1165rcbl_unl(u(U),CLP,X,Continuation,B0,B1,Indep,DepT) :-
 1166	pivot_a(CLP,X,Indep,t_U(U),DepT), % changes the basis
 1167	rcbl(Continuation,CLP,B0,B1).
 1168
 1169% narrow_u(Type,X,U)
 1170%
 1171% Narrows down the upperbound of X (type Type) to U.
 1172% Fails if Type is not t_u(_) or t_lu(_)
 1173
 1174narrow_u(t_u(_),X,U) :-
 1175	get_attr(X,clpcd_itf,Att),
 1176	setarg(2,Att,type(t_u(U))).
 1177narrow_u(t_lu(L,_),X,U) :-
 1178	get_attr(X,clpcd_itf,Att),
 1179	setarg(2,Att,type(t_lu(L,U))).
 1180
 1181% narrow_l(Type,X,L)
 1182%
 1183% Narrows down the lowerbound of X (type Type) to L.
 1184% Fails if Type is not t_l(_) or t_lu(_)
 1185
 1186narrow_l( t_l(_),    X, L) :-
 1187	get_attr(X,clpcd_itf,Att),
 1188	setarg(2,Att,type(t_l(L))).
 1189
 1190narrow_l( t_lu(_,U), X, L) :-
 1191	get_attr(X,clpcd_itf,Att),
 1192	setarg(2,Att,type(t_lu(L,U)))