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(bv_q,
   41	[
   42	    allvars/2,
   43	    backsubst/3,
   44	    backsubst_delta/4,
   45	    basis_add/2,
   46	    dec_step/2,
   47	    deref/2,
   48	    deref_var/2,
   49	    detach_bounds/1,
   50	    detach_bounds_vlv/5,
   51	    determine_active_dec/1,
   52	    determine_active_inc/1,
   53	    dump_var/6,
   54	    dump_nz/5,
   55	    export_binding/1,
   56	    get_or_add_class/2,
   57	    inc_step/2,
   58	    intro_at/3,
   59	    iterate_dec/2,
   60	    lb/3,
   61	    pivot_a/4,
   62	    pivot/5,
   63	    rcbl_status/6,
   64	    reconsider/1,
   65	    same_class/2,
   66	    solve/1,
   67	    solve_ord_x/3,
   68	    ub/3,
   69	    unconstrained/4,
   70	    var_intern/2,
   71	    var_intern/3,
   72	    var_intern/4,
   73	    var_with_def_assign/2,
   74	    var_with_def_intern/4,
   75	    maximize/1,
   76	    minimize/1,
   77	    sup/2,
   78	    sup/4,
   79	    inf/2,
   80	    inf/4,
   81	    'solve_<'/1,
   82	    'solve_=<'/1,
   83	    'solve_=\\='/1,
   84	    log_deref/4
   85	]).   86:- use_module(store_q,
   87	[
   88	    add_linear_11/3,
   89	    add_linear_f1/4,
   90	    add_linear_ff/5,
   91	    delete_factor/4,
   92	    indep/2,
   93	    isolate/3,
   94	    nf2sum/3,
   95	    nf_rhs_x/4,
   96	    nf_substitute/4,
   97	    normalize_scalar/2,
   98	    mult_hom/3,
   99	    mult_linear_factor/3
  100	]).  101:- use_module('../clpqr/class',
  102	[
  103	    class_allvars/2,
  104	    class_basis/2,
  105	    class_basis_add/3,
  106	    class_basis_drop/2,
  107	    class_basis_pivot/3,
  108	    class_new/5
  109	]).  110:- use_module(ineq_q,
  111	[
  112	    ineq/4
  113	]).  114:- use_module(nf_q,
  115	[
  116	    {}/1,
  117	    split/3,
  118	    wait_linear/3
  119	]).  120:- use_module(bb_q,
  121	[
  122	    vertex_value/2
  123	]).  124:- use_module(library(ordsets),
  125	[
  126	    ord_add_element/3
  127	]).  128
  129% For the rhs maint. the following events are important:
  130%
  131%	-) introduction of an indep var at active bound B
  132%	-) narrowing of active bound
  133%	-) swap active bound
  134%	-) pivot
  135%
  136
  137% a variables bound (L/U) can have the states:
  138%
  139%	-) t_none	no bounds
  140%	-) t_l		inactive lower bound
  141%	-) t_u		inactive upper bound
  142%	-) t_L		active lower bound
  143%	-) t_U		active upper bound
  144%	-) t_lu		inactive lower and upper bound
  145%	-) t_Lu		active lower bound and inactive upper bound
  146%	-) t_lU		inactive lower bound and active upper bound
  147
  148% ----------------------------------- deref -----------------------------------
  149%
  150
  151% deref(Lin,Lind)
  152%
  153% Makes a linear equation of the form [v(I,[])|H] into a solvable linear
  154% equation.
  155% If the variables are new, they are initialized with the linear equation X=X.
  156
  157deref(Lin,Lind) :-
  158	split(Lin,H,I),
  159	normalize_scalar(I,Nonvar),
  160	length(H,Len),
  161	log_deref(Len,H,[],Restd),
  162	add_linear_11(Nonvar,Restd,Lind).
  163
  164% log_deref(Len,[Vs|VsTail],VsTail,Res)
  165%
  166% Logarithmically converts a linear equation in normal form ([v(_,_)|_]) into a
  167% linear equation in solver form ([I,R,K*X|_]). Res contains the result, Len is
  168% the length of the part to convert and [Vs|VsTail] is a difference list
  169% containing the equation in normal form.
  170
  171log_deref(0,Vs,Vs,Lin) :-
  172	!,
  173	Lin = [0,0].
  174log_deref(1,[v(K,[X^1])|Vs],Vs,Lin) :-
  175	!,
  176	deref_var(X,Lx),
  177	mult_linear_factor(Lx,K,Lin).
  178log_deref(2,[v(Kx,[X^1]),v(Ky,[Y^1])|Vs],Vs,Lin) :-
  179	!,
  180	deref_var(X,Lx),
  181	deref_var(Y,Ly),
  182	add_linear_ff(Lx,Kx,Ly,Ky,Lin).
  183log_deref(N,V0,V2,Lin) :-
  184	P is N >> 1,
  185	Q is N - P,
  186	log_deref(P,V0,V1,Lp),
  187	log_deref(Q,V1,V2,Lq),
  188	add_linear_11(Lp,Lq,Lin).
  189
  190% deref_var(X,Lin)
  191%
  192% Returns the equation of variable X. If X is a new variable, a new equation
  193% X = X is made.
  194
  195deref_var(X,Lin) :-
  196	(   get_attr(X,clpqr_itf,Att)
  197	->  (   \+ arg(1,Att,clpq)
  198	    ->  throw(error(permission_error('mix CLP(Q) variables with',
  199		'CLP(R) variables:',X),context(_,_)))
  200	    ;   arg(4,Att,lin(Lin))
  201	    ->  true
  202	    ;   setarg(2,Att,type(t_none)),
  203		setarg(3,Att,strictness(0)),
  204		Lin = [0,0,l(X*1,Ord)],
  205		setarg(4,Att,lin(Lin)),
  206		setarg(5,Att,order(Ord))
  207	    )
  208	;   Lin = [0,0,l(X*1,Ord)],
  209	    put_attr(X,clpqr_itf,t(clpq,type(t_none),strictness(0),
  210		lin(Lin),order(Ord),n,n,n,n,n,n))
  211	).
  212
  213% TODO
  214%
  215%
  216
  217var_with_def_assign(Var,Lin) :-
  218	Lin = [I,_|Hom],
  219	(   Hom = []
  220	->  % X=k
  221	    Var = I
  222	;   Hom = [l(V*K,_)|Cs]
  223	->  (   Cs = [],
  224		K =:= 1,
  225		I =:= 0
  226	    ->	% X=Y
  227		Var = V
  228	    ;	% general case
  229		var_with_def_intern(t_none,Var,Lin,0)
  230	    )
  231	).
  232
  233% var_with_def_intern(Type,Var,Lin,Strictness)
  234%
  235% Makes Lin the linear equation of new variable Var, makes all variables of
  236% Lin, and Var of the same class and bounds Var by type(Type) and
  237% strictness(Strictness)
  238%
  239% Var is always a variable the solver invents for itself: a slack variable,
  240% the target of an optimisation or the witness of a disequation.  It is
  241% therefore marked auxiliary (argument 7), so that projection can eliminate
  242% it again rather than report it as part of an answer.
  243
  244var_with_def_intern(Type,Var,Lin,Strict) :-
  245	put_attr(Var,clpqr_itf,t(clpq,type(Type),strictness(Strict),lin(Lin),
  246	    order(_),n,aux,n,n,n,n)),	% check uses
  247	Lin = [_,_|Hom],
  248	get_or_add_class(Var,Class),
  249	same_class(Hom,Class).
  250
  251% TODO
  252%
  253%
  254
  255var_intern(Type,Var,Strict) :-
  256	var_intern(Type,Var,Strict,n).
  257
  258% var_intern(Type,Var,Strictness,Aux)
  259%
  260% As var_intern/3.  Aux is aux if Var is a slack variable introduced by the
  261% solver and n if it is a variable the user wrote down; see
  262% var_with_def_intern/4.
  263
  264var_intern(Type,Var,Strict,Aux) :-
  265	put_attr(Var,clpqr_itf,t(clpq,type(Type),strictness(Strict),
  266	    lin([0,0,l(Var*1,Ord)]),order(Ord),n,Aux,n,n,n,n)),
  267	get_or_add_class(Var,_Class).
  268
  269% TODO
  270%
  271%
  272
  273var_intern(Var,Class) :-	% for ordered/1 but otherwise free vars
  274	get_attr(Var,clpqr_itf,Att),
  275	arg(2,Att,type(_)),
  276	arg(4,Att,lin(_)),
  277	!,
  278	get_or_add_class(Var,Class).
  279var_intern(Var,Class) :-
  280	put_attr(Var,clpqr_itf,t(clpq,type(t_none),strictness(0),
  281	    lin([0,0,l(Var*1,Ord)]),order(Ord),n,n,n,n,n,n)),
  282	get_or_add_class(Var,Class).
  283
  284% -----------------------------------------------------------------------------
  285
  286% export_binding(Lst)
  287%
  288% Binds variables X to Y where Lst contains elements of the form [X-Y].
  289
  290export_binding([]).
  291export_binding([X-Y|Gs]) :-
  292	Y = X,
  293	export_binding(Gs).
  294
  295% 'solve_=\\='(Nf)
  296%
  297% Solves linear inequality Nf =\= 0 where Nf is in normal form.
  298
  299'solve_=\\='(Nf) :-
  300	deref(Nf,Lind),	% dereferences and turns Nf into solvable form Lind
  301	Lind = [Inhom,_|Hom],
  302	(   Hom = []
  303	->  Inhom =\= 0
  304	;   % make new variable Nz = Lind
  305	    var_with_def_intern(t_none,Nz,Lind,0),
  306	    % make Nz nonzero
  307	    get_attr(Nz,clpqr_itf,Att),
  308	    setarg(8,Att,nonzero)
  309	).
  310
  311% 'solve_<'(Nf)
  312%
  313% Solves linear inequality Nf < 0 where Nf is in normal form.
  314
  315'solve_<'(Nf) :-
  316	split(Nf,H,I),
  317	ineq(H,I,Nf,strict).
  318
  319% 'solve_=<'(Nf)
  320%
  321% Solves linear inequality Nf =< 0 where Nf is in normal form.
  322
  323'solve_=<'(Nf) :-
  324	split(Nf,H,I),
  325	ineq(H,I,Nf,nonstrict).
  326
  327maximize(Term) :-
  328	minimize(-Term).
  329
  330%
  331% This is NOT coded as minimize(Expr) :- inf(Expr,Expr).
  332%
  333% because the new version of inf/2 only visits
  334% the vertex where the infimum is assumed and returns
  335% to the 'current' vertex via backtracking.
  336% The rationale behind this construction is to eliminate
  337% all garbage in the solver data structures produced by
  338% the pivots on the way to the extremal point caused by
  339% {inf,sup}/{2,4}.
  340%
  341% If we are after the infimum/supremum for minimizing/maximizing,
  342% this strategy may have adverse effects on performance because
  343% the simplex algorithm is forced to re-discover the
  344% extremal vertex through the equation {Inf =:= Expr}.
  345%
  346% Thus the extra code for {minimize,maximize}/1.
  347%
  348% In case someone comes up with an example where
  349%
  350%   inf(Expr,Expr)
  351%
  352% outperforms the provided formulation for minimize - so be it.
  353% Both forms are available to the user.
  354%
  355minimize(Term) :-
  356	wait_linear(Term,Nf,minimize_lin(Nf)).
  357
  358% minimize_lin(Lin)
  359%
  360% Minimizes the linear expression Lin. It does so by making a new
  361% variable Dep and minimizes its value.
  362
  363minimize_lin(Lin) :-
  364	deref(Lin,Lind),
  365	var_with_def_intern(t_none,Dep,Lind,0),
  366	determine_active_dec(Lind),
  367	iterate_dec(Dep,Inf),
  368	{ Dep =:= Inf }.
  369
  370sup(Expression,Sup) :-
  371	sup(Expression,Sup,[],[]).
  372
  373sup(Expression,Sup,Vector,Vertex) :-
  374	inf(-Expression,-Sup,Vector,Vertex).
  375
  376inf(Expression,Inf) :-
  377	inf(Expression,Inf,[],[]).
  378
  379inf(Expression,Inf,Vector,Vertex) :-
  380	% wait until Expression becomes linear, Nf contains linear Expression
  381	% in normal form
  382	wait_linear(Expression,Nf,inf_lin(Nf,Inf,Vector,Vertex)).
  383
  384% The optimum is found by pivoting and then thrown away again by the
  385% failure driven loop, which undoes those pivots.  It is carried across in
  386% a mutable term local to this call rather than in a global variable, which
  387% would clobber a global of the same name in the calling program.
  388
  389inf_lin(Lin,Infimum,Vector,Vertex) :-
  390	State = state(none),
  391	(   deref(Lin,Lind),
  392	    var_with_def_intern(t_none,Dep,Lind,0),	% make new variable Dep = Lind
  393	    determine_active_dec(Lind),	% minimizes Lind
  394	    iterate_dec(Dep,Inf),
  395	    vertex_value(Vector,Values),
  396	    nb_setarg(1,State,[Inf|Values]),
  397	    fail
  398	;   arg(1,State,L),
  399	    L = [_|_],
  400	    assign([Infimum|Vertex],L)
  401	).
  402
  403% assign(L1,L2)
  404%
  405% The elements of L1 are pairwise assigned to the elements of L2
  406% by means of asserting {X =:= Y} where X is an element of L1 and Y
  407% is the corresponding element of L2.
  408
  409assign([],[]).
  410assign([X|Xs],[Y|Ys]) :-
  411	{X =:= Y},		  % more defensive/expressive than X=Y
  412	assign(Xs,Ys).
  413
  414% --------------------------------- optimization ------------------------------
  415%
  416% The _sn(S) =< 0 row might be temporarily infeasible.
  417% We use reconsider/1 to fix this.
  418%
  419%   s(S) e [_,0] = d +xi ... -xj, Rhs > 0 so we want to decrease s(S)
  420%
  421%   positive xi would have to be moved towards their lower bound,
  422%   negative xj would have to be moved towards their upper bound,
  423%
  424%   the row s(S) does not limit the lower bound of xi
  425%   the row s(S) does not limit the upper bound of xj
  426%
  427%   a) if some other row R is limiting xk, we pivot(R,xk),
  428%      s(S) will decrease and get more feasible until (b)
  429%   b) if there is no limiting row for some xi: we pivot(s(S),xi)
  430%					    xj: we pivot(s(S),xj)
  431%      which cures the infeasibility in one step
  432%
  433
  434
  435% iterate_dec(OptVar,Opt)
  436%
  437% Decreases the bound on the variables of the linear equation of OptVar as much
  438% as possible and returns the resulting optimal bound in Opt. Fails if for some
  439% variable, a status of unlimited is found.
  440
  441iterate_dec(OptVar,Opt) :-
  442	get_attr(OptVar,clpqr_itf,Att),
  443	arg(4,Att,lin([I,R|H])),
  444	dec_step(H,Status),
  445	(   Status = applied
  446	->  iterate_dec(OptVar,Opt)
  447	;   Status = optimum,
  448	    Opt is R + I
  449	).
  450
  451%
  452% Status = {optimum,unlimited(Indep,DepT),applied}
  453% If Status = optimum, the tables have not been changed at all.
  454% Searches left to right, does not try to find the 'best' pivot
  455% Therefore we might discover unboundedness only after a few pivots
  456%
  457
  458
  459dec_step_cont([],optimum,Cont,Cont).
  460dec_step_cont([l(V*K,OrdV)|Vs],Status,ContIn,ContOut) :-
  461	get_attr(V,clpqr_itf,Att),
  462	arg(2,Att,type(W)),
  463	arg(6,Att,class(Class)),
  464	(   dec_step_2_cont(W,l(V*K,OrdV),Class,Status,ContIn,ContOut)
  465	->  true
  466	;   dec_step_cont(Vs,Status,ContIn,ContOut)
  467	).
  468
  469inc_step_cont([],optimum,Cont,Cont).
  470inc_step_cont([l(V*K,OrdV)|Vs],Status,ContIn,ContOut) :-
  471	get_attr(V,clpqr_itf,Att),
  472	arg(2,Att,type(W)),
  473	arg(6,Att,class(Class)),
  474	(   inc_step_2_cont(W,l(V*K,OrdV),Class,Status,ContIn,ContOut)
  475	->  true
  476	;   inc_step_cont(Vs,Status,ContIn,ContOut)
  477	).
  478
  479dec_step_2_cont(t_U(U),l(V*K,OrdV),Class,Status,ContIn,ContOut) :-
  480	K > 0,
  481	(   lb(Class,OrdV,Vub-Vb-_)
  482	->  % found a lower bound
  483	    Status = applied,
  484	    pivot_a(Vub,V,Vb,t_u(U)),
  485	    replace_in_cont(ContIn,Vub,V,ContOut)
  486	;   Status = unlimited(V,t_u(U)),
  487	    ContIn = ContOut
  488	).
  489dec_step_2_cont(t_lU(L,U),l(V*K,OrdV),Class,applied,ContIn,ContOut) :-
  490	K > 0,
  491	Init is L - U,
  492	class_basis(Class,Deps),
  493	lb(Deps,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_),
  494	pivot_b(Vub,V,Vb,t_lu(L,U)),
  495	replace_in_cont(ContIn,Vub,V,ContOut).
  496dec_step_2_cont(t_L(L),l(V*K,OrdV),Class,Status,ContIn,ContOut) :-
  497	K < 0,
  498	(   ub(Class,OrdV,Vub-Vb-_)
  499	->  Status = applied,
  500	    pivot_a(Vub,V,Vb,t_l(L)),
  501	    replace_in_cont(ContIn,Vub,V,ContOut)
  502	;   Status = unlimited(V,t_l(L)),
  503	    ContIn = ContOut
  504	).
  505dec_step_2_cont(t_Lu(L,U),l(V*K,OrdV),Class,applied,ContIn,ContOut) :-
  506	K < 0,
  507	Init is U - L,
  508	class_basis(Class,Deps),
  509	ub(Deps,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_),
  510	pivot_b(Vub,V,Vb,t_lu(L,U)),
  511	replace_in_cont(ContIn,Vub,V,ContOut).
  512dec_step_2_cont(t_none,l(V*_,_),_,unlimited(V,t_none),Cont,Cont).
  513
  514
  515
  516inc_step_2_cont(t_U(U),l(V*K,OrdV),Class,Status,ContIn,ContOut) :-
  517	K < 0,
  518	(   lb(Class,OrdV,Vub-Vb-_)
  519	->  Status = applied,
  520	    pivot_a(Vub,V,Vb,t_u(U)),
  521	    replace_in_cont(ContIn,Vub,V,ContOut)
  522	;   Status = unlimited(V,t_u(U)),
  523	    ContIn = ContOut
  524	).
  525inc_step_2_cont(t_lU(L,U),l(V*K,OrdV),Class,applied,ContIn,ContOut) :-
  526	K < 0,
  527	Init is L - U,
  528	class_basis(Class,Deps),
  529	lb(Deps,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_),
  530	pivot_b(Vub,V,Vb,t_lu(L,U)),
  531	replace_in_cont(ContIn,Vub,V,ContOut).
  532inc_step_2_cont(t_L(L),l(V*K,OrdV),Class,Status,ContIn,ContOut) :-
  533	K > 0,
  534	(   ub(Class,OrdV,Vub-Vb-_)
  535	->  Status = applied,
  536	    pivot_a(Vub,V,Vb,t_l(L)),
  537	    replace_in_cont(ContIn,Vub,V,ContOut)
  538	;   Status = unlimited(V,t_l(L)),
  539	    ContIn = ContOut
  540	).
  541inc_step_2_cont(t_Lu(L,U),l(V*K,OrdV),Class,applied,ContIn,ContOut) :-
  542	K > 0,
  543	Init is U - L,
  544	class_basis(Class,Deps),
  545	ub(Deps,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_),
  546	pivot_b(Vub,V,Vb,t_lu(L,U)),
  547	replace_in_cont(ContIn,Vub,V,ContOut).
  548inc_step_2_cont(t_none,l(V*_,_),_,unlimited(V,t_none),Cont,Cont).
  549
  550replace_in_cont([],_,_,[]).
  551replace_in_cont([H1|T1],X,Y,[H2|T2]) :-
  552	(   H1 == X
  553	->  H2 = Y,
  554	    T1 = T2
  555	;   H2 = H1,
  556	    replace_in_cont(T1,X,Y,T2)
  557	).
  558
  559dec_step([],optimum).
  560dec_step([l(V*K,OrdV)|Vs],Status) :-
  561	get_attr(V,clpqr_itf,Att),
  562	arg(2,Att,type(W)),
  563	arg(6,Att,class(Class)),
  564	(   dec_step_2(W,l(V*K,OrdV),Class,Status)
  565	->  true
  566	;   dec_step(Vs,Status)
  567	).
  568
  569dec_step_2(t_U(U),l(V*K,OrdV),Class,Status) :-
  570	K > 0,
  571	(   lb(Class,OrdV,Vub-Vb-_)
  572	->  % found a lower bound
  573	    Status = applied,
  574	    pivot_a(Vub,V,Vb,t_u(U))
  575	;   Status = unlimited(V,t_u(U))
  576	).
  577dec_step_2(t_lU(L,U),l(V*K,OrdV),Class,applied) :-
  578	K > 0,
  579	Init is L - U,
  580	class_basis(Class,Deps),
  581	lb(Deps,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_),
  582	pivot_b(Vub,V,Vb,t_lu(L,U)).
  583dec_step_2(t_L(L),l(V*K,OrdV),Class,Status) :-
  584	K < 0,
  585	(   ub(Class,OrdV,Vub-Vb-_)
  586	->  Status = applied,
  587	    pivot_a(Vub,V,Vb,t_l(L))
  588	;   Status = unlimited(V,t_l(L))
  589	).
  590dec_step_2(t_Lu(L,U),l(V*K,OrdV),Class,applied) :-
  591	K < 0,
  592	Init is U - L,
  593	class_basis(Class,Deps),
  594	ub(Deps,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_),
  595	pivot_b(Vub,V,Vb,t_lu(L,U)).
  596dec_step_2(t_none,l(V*_,_),_,unlimited(V,t_none)).
  597
  598inc_step([],optimum).	% if status has not been set yet: no changes
  599inc_step([l(V*K,OrdV)|Vs],Status) :-
  600	get_attr(V,clpqr_itf,Att),
  601	arg(2,Att,type(W)),
  602	arg(6,Att,class(Class)),
  603	(   inc_step_2(W,l(V*K,OrdV),Class,Status)
  604	->  true
  605	;   inc_step(Vs,Status)
  606	).
  607
  608inc_step_2(t_U(U),l(V*K,OrdV),Class,Status) :-
  609	K < 0,
  610	(   lb(Class,OrdV,Vub-Vb-_)
  611	->  Status = applied,
  612	    pivot_a(Vub,V,Vb,t_u(U))
  613	;   Status = unlimited(V,t_u(U))
  614	).
  615inc_step_2(t_lU(L,U),l(V*K,OrdV),Class,applied) :-
  616	K < 0,
  617	Init is L - U,
  618	class_basis(Class,Deps),
  619	lb(Deps,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_),
  620	pivot_b(Vub,V,Vb,t_lu(L,U)).
  621inc_step_2(t_L(L),l(V*K,OrdV),Class,Status) :-
  622	K > 0,
  623	(   ub(Class,OrdV,Vub-Vb-_)
  624	->  Status = applied,
  625	    pivot_a(Vub,V,Vb,t_l(L))
  626	;   Status = unlimited(V,t_l(L))
  627	).
  628inc_step_2(t_Lu(L,U),l(V*K,OrdV),Class,applied) :-
  629	K > 0,
  630	Init is U - L,
  631	class_basis(Class,Deps),
  632	ub(Deps,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_),
  633	pivot_b(Vub,V,Vb,t_lu(L,U)).
  634inc_step_2(t_none,l(V*_,_),_,unlimited(V,t_none)).
  635
  636% ------------------------- find the most constraining row --------------------
  637%
  638% The code for the lower and the upper bound are dual versions of each other.
  639% The only difference is in the orientation of the comparisons.
  640% Indeps are ruled out by their types.
  641% If there is no bound, this fails.
  642%
  643% *** The actual lb and ub on an indep variable X are [lu]b + b(X), where b(X)
  644% is the value of the active bound.
  645%
  646% Nota bene: We must NOT consider infeasible rows as candidates to
  647%	     leave the basis!
  648%
  649% ub(Class,OrdX,Ub)
  650%
  651% See lb/3: this is similar
  652
  653ub(Class,OrdX,Ub) :-
  654	class_basis(Class,Deps),
  655	ub_first(Deps,OrdX,Ub).
  656
  657% ub_first(Deps,X,Dep-W-Ub)
  658%
  659% Finds the tightest upperbound for variable X from the linear equations of
  660% basis variables Deps, and puts the resulting bound in Ub. Dep is the basis
  661% variable that generates the bound, and W is bound of that variable that has
  662% to be activated to achieve this.
  663
  664ub_first([Dep|Deps],OrdX,Tightest) :-
  665	(   get_attr(Dep,clpqr_itf,Att),
  666	    arg(2,Att,type(Type)),
  667	    arg(4,Att,lin(Lin)),
  668	    ub_inner(Type,OrdX,Lin,W,Ub),
  669	    Ub >= 0
  670	->  ub(Deps,OrdX,Dep-W-Ub,Tightest)
  671	;   ub_first(Deps,OrdX,Tightest)
  672	).
  673
  674% ub(Deps,OrdX,TightestIn,TightestOut)
  675%
  676% See lb/4: this is similar
  677
  678ub([],_,T0,T0).
  679ub([Dep|Deps],OrdX,T0,T1) :-
  680	(   get_attr(Dep,clpqr_itf,Att),
  681	    arg(2,Att,type(Type)),
  682	    arg(4,Att,lin(Lin)),
  683	    ub_inner(Type,OrdX,Lin,W,Ub),
  684	    T0 = _-Ubb,
  685	    Ub < Ubb,
  686	    Ub >= 0
  687	->  ub(Deps,OrdX,Dep-W-Ub,T1)	% tighter bound, use new bound
  688	;   ub(Deps,OrdX,T0,T1)	% no tighter bound, keep current one
  689	).
  690
  691% ub_inner(Type,OrdX,Lin,W,Ub)
  692%
  693% See lb_inner/5: this is similar
  694
  695ub_inner(t_l(L),OrdX,Lin,t_L(L),Ub) :-
  696	nf_rhs_x(Lin,OrdX,Rhs,K),
  697	K < 0,
  698	Ub is (L - Rhs) rdiv K.
  699ub_inner(t_u(U),OrdX,Lin,t_U(U),Ub) :-
  700	nf_rhs_x(Lin,OrdX,Rhs,K),
  701	K > 0,
  702	Ub is (U - Rhs) rdiv K.
  703ub_inner(t_lu(L,U),OrdX,Lin,W,Ub) :-
  704	nf_rhs_x(Lin,OrdX,Rhs,K),
  705	(   K < 0 % use lowerbound
  706	->  W = t_Lu(L,U),
  707	    Ub = (L - Rhs) rdiv K
  708	;   K > 0 % use upperbound
  709	->  W = t_lU(L,U),
  710	    Ub = (U - Rhs) rdiv K
  711	).
  712
  713% lb(Class,OrdX,Lb)
  714%
  715% Returns in Lb how much we can lower the upperbound of X without violating
  716% a bound of the basisvariables.
  717% Lb has the form Dep-W-Lb with Dep the variable whose bound is violated when
  718% lowering the bound for X more, W the actual bound that has to be activated
  719% and Lb the amount that the upperbound can be lowered.
  720% X has ordering OrdX and class Class.
  721
  722lb(Class,OrdX,Lb) :-
  723	class_basis(Class,Deps),
  724	lb_first(Deps,OrdX,Lb).
  725
  726% lb_first(Deps,OrdX,Tightest)
  727%
  728% Returns in Tightest how much we can lower the upperbound of X without
  729% violating a bound of Deps.
  730% Tightest has the form Dep-W-Lb with Dep the variable whose bound is violated
  731% when lowering the bound for X more, W the actual bound that has to be
  732% activated and Lb the amount that the upperbound can be lowered. X has
  733% ordering attribute OrdX.
  734
  735lb_first([Dep|Deps],OrdX,Tightest) :-
  736	(   get_attr(Dep,clpqr_itf,Att),
  737	    arg(2,Att,type(Type)),
  738	    arg(4,Att,lin(Lin)),
  739	    lb_inner(Type,OrdX,Lin,W,Lb),
  740	    Lb =< 0 % Lb > 0 means a violated bound
  741	->  lb(Deps,OrdX,Dep-W-Lb,Tightest)
  742	;   lb_first(Deps,OrdX,Tightest)
  743	).
  744
  745% lb(Deps,OrdX,TightestIn,TightestOut)
  746%
  747% See lb_first/3: this one does the same thing, but is used for the steps after
  748% the first one and remembers the tightest bound so far.
  749
  750lb([],_,T0,T0).
  751lb([Dep|Deps],OrdX,T0,T1) :-
  752	(   get_attr(Dep,clpqr_itf,Att),
  753	    arg(2,Att,type(Type)),
  754	    arg(4,Att,lin(Lin)),
  755	    lb_inner(Type,OrdX,Lin,W,Lb),
  756	    T0 = _-Lbb,
  757	    Lb > Lbb,	% choose the least lowering, others might violate
  758			% bounds
  759	    Lb =< 0	% violation of a bound (without lowering)
  760	->  lb(Deps,OrdX,Dep-W-Lb,T1)
  761	;   lb(Deps,OrdX,T0,T1)
  762	).
  763
  764% lb_inner(Type,X,Lin,W,Lb)
  765%
  766% Returns in Lb how much lower we can make X without violating a bound
  767% by using the linear equation Lin of basis variable B which has type
  768% Type and which has to activate a bound (type W) to do so.
  769%
  770% E.g. when B has a lowerbound L, then L should always be smaller than I + R.
  771% So a lowerbound of X (which has scalar K in Lin), could be at most
  772% (L-(I+R))/K lower than its upperbound (if K is positive).
  773% Also note that Lb should always be smaller than 0, otherwise the row is
  774% not feasible.
  775% X has ordering attribute OrdX.
  776
  777lb_inner(t_l(L),OrdX,Lin,t_L(L),Lb) :-
  778	nf_rhs_x(Lin,OrdX,Rhs,K), % if linear equation Lin contains the term
  779				  % X*K, Rhs is the right hand side of that
  780				  % equation
  781	K > 0,
  782	Lb is (L - Rhs) rdiv K.
  783lb_inner(t_u(U),OrdX,Lin,t_U(U),Lb) :-
  784	nf_rhs_x(Lin,OrdX,Rhs,K),
  785	K < 0, % K < 0
  786	Lb is (U - Rhs) rdiv K.
  787lb_inner(t_lu(L,U),OrdX,Lin,W,Lb) :-
  788	nf_rhs_x(Lin,OrdX,Rhs,K),
  789	(   K < 0
  790	->  W = t_lU(L,U),
  791	    Lb is (U - Rhs) rdiv K
  792	;   K > 0
  793	->  W = t_Lu(L,U),
  794	    Lb is (L - Rhs) rdiv K
  795	).
  796
  797% ---------------------------------- equations --------------------------------
  798%
  799% backsubstitution will not make the system infeasible, if the bounds on the
  800% indep vars are obeyed, but some implied values might pop up in rows where X
  801% occurs
  802%	-) special case X=Y during bs -> get rid of dependend var(s), alias
  803%
  804
  805solve(Lin) :-
  806	Lin = [I,_|H],
  807	solve(H,Lin,I,Bindings,[]),
  808	export_binding(Bindings).
  809
  810% solve(Hom,Lin,I,Bind,BindT)
  811%
  812% Solves a linear equation Lin = [I,_|H] = 0 and exports the generated bindings
  813
  814solve([],_,I,Bind0,Bind0) :-
  815	!,
  816	I =:= 0.
  817solve(H,Lin,_,Bind0,BindT) :-
  818	sd(H,[],ClassesUniq,9-9-0,Category-Selected-_,NV,NVT),
  819	get_attr(Selected,clpqr_itf,Att),
  820	arg(5,Att,order(Ord)),
  821	isolate(Ord,Lin,Lin1),	% Lin = 0 => Selected = Lin1
  822	(   Category = 1 % classless variable, no bounds
  823	->  setarg(4,Att,lin(Lin1)),
  824	    Lin1 = [Inhom,_|Hom],
  825	    bs_collect_binding(Hom,Selected,Inhom,Bind0,BindT),
  826	    eq_classes(NV,NVT,ClassesUniq)
  827	;   Category = 2 % class variable, no bounds
  828	->  arg(6,Att,class(NewC)),
  829	    class_allvars(NewC,Deps),
  830	    (   ClassesUniq = [_] % rank increasing
  831	    ->	bs_collect_bindings(Deps,Ord,Lin1,Bind0,BindT)
  832	    ;   Bind0 = BindT,
  833		bs(Deps,Ord,Lin1)
  834	    ),
  835	    eq_classes(NV,NVT,ClassesUniq)
  836	;   Category = 3 % classless variable, all variables in Lin and
  837			 % Selected are bounded
  838	->  arg(2,Att,type(Type)),
  839	    setarg(4,Att,lin(Lin1)),
  840	    deactivate_bound(Type,Selected),
  841	    eq_classes(NV,NVT,ClassesUniq),
  842	    basis_add(Selected,Basis),
  843	    undet_active(Lin1),	% we can't tell which bound will likely be a
  844				% problem at this point
  845	    Lin1 = [Inhom,_|Hom],
  846	    bs_collect_binding(Hom,Selected,Inhom,Bind0,Bind1),	% only if
  847								% Hom = []
  848	    rcbl(Basis,Bind1,BindT) % reconsider entire basis
  849	;   Category = 4 % class variable, all variables in Lin and Selected
  850			 % are bounded
  851	->  arg(2,Att,type(Type)),
  852	    arg(6,Att,class(NewC)),
  853	    class_allvars(NewC,Deps),
  854	    (   ClassesUniq = [_] % rank increasing
  855	    ->	bs_collect_bindings(Deps,Ord,Lin1,Bind0,Bind1)
  856	    ;   Bind0 = Bind1,
  857		bs(Deps,Ord,Lin1)
  858	    ),
  859	    deactivate_bound(Type,Selected),
  860	    basis_add(Selected,Basis),
  861	    % eq_classes( NV, NVT, ClassesUniq),
  862	    %  4 -> var(NV)
  863	    equate(ClassesUniq,_),
  864	    undet_active(Lin1),
  865	    rcbl(Basis,Bind1,BindT)
  866	).
  867
  868% solve_ord_x(Lin,OrdX,ClassX)
  869%
  870% Like solve/1, but solves for the particular variable with ordering OrdX,
  871% whose class ClassX is known.
  872
  873solve_ord_x(Lin,OrdX,ClassX) :-
  874	Lin = [I,_|H],
  875	solve_ord_x(H,Lin,I,OrdX,ClassX,Bindings,[]),
  876	export_binding(Bindings).
  877
  878solve_ord_x([],_,I,_,_,Bind0,Bind0) :-
  879	I =:= 0.
  880solve_ord_x([_|_],Lin,_,OrdX,ClassX,Bind0,BindT) :-
  881	isolate(OrdX,Lin,Lin1),
  882	Lin1 = [_,_|H1],
  883	sd(H1,[],ClassesUniq1,9-9-0,_,NV,NVT), % do sd on Lin without X, then
  884					       % add class of X
  885	ord_add_element(ClassesUniq1,ClassX,ClassesUniq),
  886	class_allvars(ClassX,Deps),
  887	(   ClassesUniq = [_] % rank increasing
  888	->  bs_collect_bindings(Deps,OrdX,Lin1,Bind0,BindT)
  889	;   Bind0 = BindT,
  890	    bs(Deps,OrdX,Lin1)
  891	),
  892	eq_classes(NV,NVT,ClassesUniq).
  893
  894% sd(H,[],ClassesUniq,9-9-0,Category-Selected-_,NV,NVT)
  895
  896% sd(Hom,ClassesIn,ClassesOut,PreferenceIn,PreferenceOut,[NV|NVTail],NVTail)
  897%
  898% ClassesOut is a sorted list of the different classes that are either in
  899% ClassesIn or that are the classes of the variables in Hom. Variables that do
  900% not belong to a class yet, are put in the difference list NV.
  901
  902sd([],Class0,Class0,Preference0,Preference0,NV0,NV0).
  903sd([l(X*K,_)|Xs],Class0,ClassN,Preference0,PreferenceN,NV0,NVt) :-
  904	get_attr(X,clpqr_itf,Att),
  905	(   arg(6,Att,class(Xc)) % old: has class
  906	->  NV0 = NV1,
  907	    ord_add_element(Class0,Xc,Class1),
  908	    (   arg(2,Att,type(t_none))
  909	    ->  preference(Preference0,2-X-K,Preference1)
  910		    % has class, no bounds => category 2
  911	    ;   preference(Preference0,4-X-K,Preference1)
  912		    % has class, is bounded => category 4
  913	    )
  914	;   % new: has no class
  915	    Class1 = Class0,
  916	    NV0 = [X|NV1], % X has no class yet, add to list of new variables
  917	    (   arg(2,Att,type(t_none))
  918	    ->  preference(Preference0,1-X-K,Preference1)
  919		    % no class, no bounds => category 1
  920	    ;   preference(Preference0,3-X-K,Preference1)
  921		    % no class, is bounded => category 3
  922	    )
  923	),
  924	sd(Xs,Class1,ClassN,Preference1,PreferenceN,NV1,NVt).
  925
  926%
  927% A is best sofar, B is current
  928% smallest prefered
  929preference(A,B,Pref) :-
  930	A = Px-_-_,
  931	B = Py-_-_,
  932	(   Px < Py
  933	->  Pref = A
  934	;   Pref = B
  935	).
  936
  937% eq_classes(NV,NVTail,Cs)
  938%
  939% Attaches all classless variables NV to a new class and equates all other
  940% classes with this class. The equate operation only happens after attach_class
  941% because the unification of classes can bind the tail of the AllVars attribute
  942% to a nonvar and then the attach_class operation wouldn't work.
  943
  944eq_classes(NV,_,Cs) :-
  945	var(NV),
  946	!,
  947	equate(Cs,_).
  948eq_classes(NV,NVT,Cs) :-
  949	class_new(Su,clpq,NV,NVT,[]), % make a new class Su with NV as the variables
  950	attach_class(NV,Su), % attach the variables NV to Su
  951	equate(Cs,Su).
  952
  953equate([],_).
  954equate([X|Xs],X) :- equate(Xs,X).
  955
  956%
  957% assert: none of the Vars has a class attribute yet
  958%
  959attach_class(Xs,_) :-
  960	var(Xs), % Tail
  961	!.
  962attach_class([X|Xs],Class) :-
  963	get_attr(X,clpqr_itf,Att),
  964	setarg(6,Att,class(Class)),
  965	attach_class(Xs,Class).
  966
  967% unconstrained(Lin,Uc,Kuc,Rest)
  968%
  969% Finds an unconstrained variable Uc (type(t_none)) in Lin with scalar Kuc and
  970% removes it from Lin to return Rest.
  971
  972unconstrained(Lin,Uc,Kuc,Rest) :-
  973	Lin = [_,_|H],
  974	sd(H,[],_,9-9-0,Category-Uc-_,_,_),
  975	Category =< 2,
  976	get_attr(Uc,clpqr_itf,Att),
  977	arg(5,Att,order(OrdUc)),
  978	delete_factor(OrdUc,Lin,Rest,Kuc).
  979
  980%
  981% point the vars in Lin into the same equivalence class
  982% maybe join some global data
  983%
  984same_class([],_).
  985same_class([l(X*_,_)|Xs],Class) :-
  986	get_or_add_class(X,Class),
  987	same_class(Xs,Class).
  988
  989% get_or_add_class(X,Class)
  990%
  991% Returns in Class the class of X if X has one, or a new class where X now
  992% belongs to if X didn't have one.
  993
  994get_or_add_class(X,Class) :-
  995	get_attr(X,clpqr_itf,Att),
  996	arg(1,Att,CLP),
  997	(   arg(6,Att,class(ClassX))
  998	->  ClassX = Class
  999	;   setarg(6,Att,class(Class)),
 1000	    class_new(Class,CLP,[X|Tail],Tail,[])
 1001	).
 1002
 1003% allvars(X,Allvars)
 1004%
 1005% Allvars is a list of all variables in the class to which X belongs.
 1006
 1007allvars(X,Allvars) :-
 1008	get_attr(X,clpqr_itf,Att),
 1009	arg(6,Att,class(C)),
 1010	class_allvars(C,Allvars).
 1011
 1012% deactivate_bound(Type,Variable)
 1013%
 1014% The Type of the variable is changed to reflect the deactivation of its
 1015% bounds.
 1016% t_L(_) becomes t_l(_), t_lU(_,_) becomes t_lu(_,_) and so on.
 1017
 1018deactivate_bound(t_l(_),_).
 1019deactivate_bound(t_u(_),_).
 1020deactivate_bound(t_lu(_,_),_).
 1021deactivate_bound(t_L(L),X) :-
 1022	get_attr(X,clpqr_itf,Att),
 1023	setarg(2,Att,type(t_l(L))).
 1024deactivate_bound(t_Lu(L,U),X) :-
 1025	get_attr(X,clpqr_itf,Att),
 1026	setarg(2,Att,type(t_lu(L,U))).
 1027deactivate_bound(t_U(U),X) :-
 1028	get_attr(X,clpqr_itf,Att),
 1029	setarg(2,Att,type(t_u(U))).
 1030deactivate_bound(t_lU(L,U),X) :-
 1031	get_attr(X,clpqr_itf,Att),
 1032	setarg(2,Att,type(t_lu(L,U))).
 1033
 1034% intro_at(X,Value,Type)
 1035%
 1036% Variable X gets new type Type which reflects the activation of a bound with
 1037% value Value. In the linear equations of all the variables belonging to the
 1038% same class as X, X is substituted by [0,Value,X] to reflect the new active
 1039% bound.
 1040
 1041intro_at(X,Value,Type) :-
 1042	get_attr(X,clpqr_itf,Att),
 1043	arg(5,Att,order(Ord)),
 1044	arg(6,Att,class(Class)),
 1045	setarg(2,Att,type(Type)),
 1046	(   Value =:= 0
 1047	->  true
 1048	;   backsubst_delta(Class,Ord,X,Value)
 1049	).
 1050
 1051% undet_active(Lin)
 1052%
 1053% For each variable in the homogene part of Lin, a bound is activated
 1054% if an inactive bound exists. (t_l(L) becomes t_L(L) and so on)
 1055
 1056undet_active([_,_|H]) :-
 1057	undet_active_h(H).
 1058
 1059% undet_active_h(Hom)
 1060%
 1061% For each variable in homogene part Hom, a bound is activated if an
 1062% inactive bound exists (t_l(L) becomes t_L(L) and so on)
 1063
 1064undet_active_h([]).
 1065undet_active_h([l(X*_,_)|Xs]) :-
 1066	get_attr(X,clpqr_itf,Att),
 1067	arg(2,Att,type(Type)),
 1068	undet_active(Type,X),
 1069	undet_active_h(Xs).
 1070
 1071% undet_active(Type,Var)
 1072%
 1073% An inactive bound of Var is activated if such exists
 1074% t_lu(L,U) is arbitrarily chosen to become t_Lu(L,U)
 1075
 1076undet_active(t_none,_).	% type_activity
 1077undet_active(t_L(_),_).
 1078undet_active(t_Lu(_,_),_).
 1079undet_active(t_U(_),_).
 1080undet_active(t_lU(_,_),_).
 1081undet_active(t_l(L),X) :- intro_at(X,L,t_L(L)).
 1082undet_active(t_u(U),X) :- intro_at(X,U,t_U(U)).
 1083undet_active(t_lu(L,U),X) :- intro_at(X,L,t_Lu(L,U)).
 1084
 1085% determine_active_dec(Lin)
 1086%
 1087% Activates inactive bounds on the variables of Lin if such bounds exist.
 1088% If the type of a variable is t_none, this fails. This version is aimed
 1089% to make the R component of Lin as small as possible in order not to violate
 1090% an upperbound (see reconsider/1)
 1091
 1092determine_active_dec([_,_|H]) :-
 1093	determine_active(H,-1).
 1094
 1095% determine_active_inc(Lin)
 1096%
 1097% Activates inactive bounds on the variables of Lin if such bounds exist.
 1098% If the type of a variable is t_none, this fails. This version is aimed
 1099% to make the R component of Lin as large as possible in order not to violate
 1100% a lowerbound (see reconsider/1)
 1101
 1102determine_active_inc([_,_|H]) :-
 1103	determine_active(H,1).
 1104
 1105% determine_active(Hom,S)
 1106%
 1107% For each variable in Hom, activates its bound if it is not yet activated.
 1108% For the case of t_lu(_,_) the lower or upper bound is activated depending on
 1109% K and S:
 1110% If sign of K*S is negative, then lowerbound, otherwise upperbound.
 1111
 1112determine_active([],_).
 1113determine_active([l(X*K,_)|Xs],S) :-
 1114	get_attr(X,clpqr_itf,Att),
 1115	arg(2,Att,type(Type)),
 1116	determine_active(Type,X,K,S),
 1117	determine_active(Xs,S).
 1118
 1119determine_active(t_L(_),_,_,_).
 1120determine_active(t_Lu(_,_),_,_,_).
 1121determine_active(t_U(_),_,_,_).
 1122determine_active(t_lU(_,_),_,_,_).
 1123determine_active(t_l(L),X,_,_) :- intro_at(X,L,t_L(L)).
 1124determine_active(t_u(U),X,_,_) :- intro_at(X,U,t_U(U)).
 1125determine_active(t_lu(L,U),X,K,S) :-
 1126	KS is K*S,
 1127	(   KS < 0
 1128	->  intro_at(X,L,t_Lu(L,U))
 1129	;   KS > 0
 1130	->  intro_at(X,U,t_lU(L,U))
 1131	).
 1132
 1133%
 1134% Careful when an indep turns into t_none !!!
 1135%
 1136
 1137detach_bounds(V) :-
 1138	get_attr(V,clpqr_itf,Att),
 1139	arg(2,Att,type(Type)),
 1140	arg(4,Att,lin(Lin)),
 1141	arg(5,Att,order(OrdV)),
 1142	arg(6,Att,class(Class)),
 1143	setarg(2,Att,type(t_none)),
 1144	setarg(3,Att,strictness(0)),
 1145	(   indep(Lin,OrdV)
 1146	->  (   ub(Class,OrdV,Vub-Vb-_)
 1147	    ->	% exchange against thightest
 1148		class_basis_drop(Class,Vub),
 1149		pivot(Vub,Class,OrdV,Vb,Type)
 1150	    ;   lb(Class,OrdV,Vlb-Vb-_)
 1151	    ->  class_basis_drop(Class,Vlb),
 1152		pivot(Vlb,Class,OrdV,Vb,Type)
 1153	    ;   true
 1154	    )
 1155	;   class_basis_drop(Class,V)
 1156	).
 1157
 1158detach_bounds_vlv(OrdV,Lin,Class,Var,NewLin) :-
 1159	(   indep(Lin,OrdV)
 1160	->  Lin = [_,R|_],
 1161	    (   ub(Class,OrdV,Vub-Vb-_)
 1162	    ->  % in verify_lin, class might contain two occurrences of Var,
 1163		% but it doesn't matter which one we delete
 1164		class_basis_drop(Class,Var),
 1165		pivot_vlv(Vub,Class,OrdV,Vb,R,NewLin)
 1166	    ;   lb(Class,OrdV,Vlb-Vb-_)
 1167	    ->  class_basis_drop(Class,Var),
 1168		pivot_vlv(Vlb,Class,OrdV,Vb,R,NewLin)
 1169	    ;   NewLin = Lin
 1170	    )
 1171	;   NewLin = Lin,
 1172	    class_basis_drop(Class,Var)
 1173	).
 1174
 1175% ----------------------------- manipulate the basis --------------------------
 1176
 1177% basis_add(X,NewBasis)
 1178%
 1179% NewBasis is the result of adding X to the basis of the class to which X
 1180% belongs.
 1181
 1182basis_add(X,NewBasis) :-
 1183	get_attr(X,clpqr_itf,Att),
 1184	arg(6,Att,class(Cv)),
 1185	class_basis_add(Cv,X,NewBasis).
 1186
 1187% basis_pivot(Leave,Enter)
 1188%
 1189% Removes Leave from the basis of the class to which it belongs, and adds
 1190% Enter to that basis.
 1191
 1192basis_pivot(Leave,Enter) :-
 1193	get_attr(Leave,clpqr_itf,Att),
 1194	arg(6,Att,class(Cv)),
 1195	class_basis_pivot(Cv,Enter,Leave).
 1196
 1197% ----------------------------------- pivot -----------------------------------
 1198
 1199% pivot_a(Dep,Indep,IndepT,DepT)
 1200%
 1201% Removes Dep from the basis, puts Indep in, and pivots the equation of
 1202% Dep to become one of Indep. The type of Dep becomes DepT (which means
 1203% it gets deactivated), the type of Indep becomes IndepT (which means it
 1204% gets activated)
 1205
 1206
 1207pivot_a(Dep,Indep,Vb,Wd) :-
 1208	basis_pivot(Dep,Indep),
 1209	get_attr(Indep,clpqr_itf,Att),
 1210	arg(2,Att,type(Type)),
 1211	arg(5,Att,order(Ord)),
 1212	arg(6,Att,class(Class)),
 1213	pivot(Dep,Class,Ord,Vb,Type),
 1214	get_attr(Indep,clpqr_itf,Att2), %changed?
 1215	setarg(2,Att2,type(Wd)).
 1216
 1217pivot_b(Vub,V,Vb,Wd) :-
 1218	(   Vub == V
 1219	->  get_attr(V,clpqr_itf,Att),
 1220	    arg(5,Att,order(Ord)),
 1221	    arg(6,Att,class(Class)),
 1222	    setarg(2,Att,type(Vb)),
 1223	    pivot_b_delta(Vb,Delta), % nonzero(Delta)
 1224	    backsubst_delta(Class,Ord,V,Delta)
 1225	;   pivot_a(Vub,V,Vb,Wd)
 1226	).
 1227
 1228pivot_b_delta(t_Lu(L,U),Delta) :- Delta is L-U.
 1229pivot_b_delta(t_lU(L,U),Delta) :- Delta is U-L.
 1230
 1231% select_active_bound(Type,Bound)
 1232%
 1233% Returns the bound that is active in Type (if such exists, 0 otherwise)
 1234
 1235select_active_bound(t_L(L),L).
 1236select_active_bound(t_Lu(L,_),L).
 1237select_active_bound(t_U(U),U).
 1238select_active_bound(t_lU(_,U),U).
 1239select_active_bound(t_none,0).
 1240%
 1241% for project.pl
 1242%
 1243select_active_bound(t_l(_),0).
 1244select_active_bound(t_u(_),0).
 1245select_active_bound(t_lu(_,_),0).
 1246
 1247
 1248% pivot(Dep,Class,IndepOrd,DepAct,IndAct)
 1249%
 1250% See pivot/2.
 1251% In addition, variable Indep with ordering IndepOrd has an active bound IndAct
 1252
 1253%
 1254%
 1255% Pivot taking care of rhs and active states
 1256%
 1257pivot(Dep,Class,IndepOrd,DepAct,IndAct) :-
 1258	get_attr(Dep,clpqr_itf,Att),
 1259	arg(4,Att,lin(H)),
 1260	arg(5,Att,order(DepOrd)),
 1261	setarg(2,Att,type(DepAct)),
 1262	select_active_bound(DepAct,AbvD), % New current value for Dep
 1263	select_active_bound(IndAct,AbvI), % Old current value of Indep
 1264	delete_factor(IndepOrd,H,H0,Coeff), % Dep = ... + Coeff*Indep + ...
 1265	AbvDm is -AbvD,
 1266	AbvIm is -AbvI,
 1267	add_linear_f1([0,AbvIm],Coeff,H0,H1),
 1268	K is -1 rdiv Coeff,
 1269	add_linear_ff(H1,K,[0,AbvDm,l(Dep* -1,DepOrd)],K,H2),
 1270	    % Indep = -1/Coeff*... + 1/Coeff*Dep
 1271	add_linear_11(H2,[0,AbvIm],Lin),
 1272	backsubst(Class,IndepOrd,Lin).
 1273
 1274% Rewrite Dep = ... + Coeff*Indep + ...
 1275% into Indep = ... + -1/Coeff*Dep + ...
 1276%
 1277% For backsubstitution, old current value of Indep must be removed from RHS
 1278% New current value of Dep must be added to RHS
 1279% For solving: old current value of Indep should be out of RHS
 1280
 1281pivot_vlv(Dep,Class,IndepOrd,DepAct,AbvI,Lin) :-
 1282	get_attr(Dep,clpqr_itf,Att),
 1283	arg(4,Att,lin(H)),
 1284	arg(5,Att,order(DepOrd)),
 1285	setarg(2,Att,type(DepAct)),
 1286	select_active_bound(DepAct,AbvD), % New current value for Dep
 1287	delete_factor(IndepOrd,H,H0,Coeff), % Dep = ... + Coeff*Indep + ...
 1288	AbvDm is -AbvD,
 1289	AbvIm is -AbvI,
 1290	add_linear_f1([0,AbvIm],Coeff,H0,H1),
 1291	K is -1 rdiv Coeff,
 1292	add_linear_ff(H1,K,[0,AbvDm,l(Dep* -1,DepOrd)],K,Lin),
 1293	    % Indep = -1/Coeff*... + 1/Coeff*Dep
 1294	add_linear_11(Lin,[0,AbvIm],SubstLin),
 1295	backsubst(Class,IndepOrd,SubstLin).
 1296
 1297% backsubst_delta(Class,OrdX,X,Delta)
 1298%
 1299% X with ordering attribute OrdX, is substituted in all linear equations of
 1300% variables in the class Class, by linear equation [0,Delta,l(X*1,OrdX)]. This
 1301% reflects the activation of a bound.
 1302
 1303backsubst_delta(Class,OrdX,X,Delta) :-
 1304	backsubst(Class,OrdX,[0,Delta,l(X*1,OrdX)]).
 1305
 1306% backsubst(Class,OrdX,Lin)
 1307%
 1308% X with ordering OrdX is substituted in all linear equations of variables in
 1309% the class Class, by linear equation Lin
 1310
 1311backsubst(Class,OrdX,Lin) :-
 1312	class_allvars(Class,Allvars),
 1313	bs(Allvars,OrdX,Lin).
 1314
 1315% bs(Vars,OrdV,Lin)
 1316%
 1317% In all linear equations of the variables Vars, variable V with ordering
 1318% attribute OrdV is substituted by linear equation Lin.
 1319%
 1320% valid if nothing will go ground
 1321%
 1322
 1323bs(Xs,_,_) :-
 1324	var(Xs),
 1325	!.
 1326bs([X|Xs],OrdV,Lin) :-
 1327	(   get_attr(X,clpqr_itf,Att),
 1328	    arg(4,Att,lin(LinX)),
 1329	    nf_substitute(OrdV,Lin,LinX,LinX1) % does not change attributes
 1330	->  setarg(4,Att,lin(LinX1)),
 1331	    bs(Xs,OrdV,Lin)
 1332	;   bs(Xs,OrdV,Lin)
 1333	).
 1334
 1335%
 1336% rank increasing backsubstitution
 1337%
 1338
 1339% bs_collect_bindings(Deps,SelectedOrd,Lin,Bind,BindT)
 1340%
 1341% Collects bindings (of the form [X-I] where X = I is the binding) by
 1342% substituting Selected in all linear equations of the variables Deps (which
 1343% are of the same class), by Lin. Selected has ordering attribute SelectedOrd.
 1344%
 1345% E.g. when V = 2X + 3Y + 4, X = 3V + 2Z and Y = 4X + 3
 1346% we can substitute V in the linear equation of X: X = 6X + 9Y + 2Z + 12
 1347% we can't substitute V in the linear equation of Y of course.
 1348
 1349bs_collect_bindings(Xs,_,_,Bind0,BindT) :-
 1350	var(Xs),
 1351	!,
 1352	Bind0 = BindT.
 1353bs_collect_bindings([X|Xs],OrdV,Lin,Bind0,BindT) :-
 1354	(   get_attr(X,clpqr_itf,Att),
 1355	    arg(4,Att,lin(LinX)),
 1356	    nf_substitute(OrdV,Lin,LinX,LinX1) % does not change attributes
 1357	->  setarg(4,Att,lin(LinX1)),
 1358	    LinX1 = [Inhom,_|Hom],
 1359	    bs_collect_binding(Hom,X,Inhom,Bind0,Bind1),
 1360	    bs_collect_bindings(Xs,OrdV,Lin,Bind1,BindT)
 1361	;   bs_collect_bindings(Xs,OrdV,Lin,Bind0,BindT)
 1362	).
 1363
 1364% bs_collect_binding(Hom,Selected,Inhom,Bind,BindT)
 1365%
 1366% Collects binding following from Selected = Hom + Inhom.
 1367% If Hom = [], returns the binding Selected-Inhom (=0)
 1368%
 1369bs_collect_binding([],X,Inhom) --> [X-Inhom].
 1370bs_collect_binding([_|_],_,_) --> [].
 1371
 1372%
 1373% reconsider the basis
 1374%
 1375
 1376% rcbl(Basis,Bind,BindT)
 1377%
 1378%
 1379
 1380rcbl([],Bind0,Bind0).
 1381rcbl([X|Continuation],Bind0,BindT) :-
 1382	(   rcb_cont(X,Status,Violated,Continuation,NewContinuation) % have a culprit
 1383	->  rcbl_status(Status,X,NewContinuation,Bind0,BindT,Violated)
 1384	;   rcbl(Continuation,Bind0,BindT)
 1385	).
 1386
 1387rcb_cont(X,Status,Violated,ContIn,ContOut) :-
 1388	get_attr(X,clpqr_itf,Att),
 1389	arg(2,Att,type(Type)),
 1390	arg(4,Att,lin([I,R|H])),
 1391	(   Type = t_l(L) % case 1: lowerbound: R + I should always be larger
 1392			  % than the lowerbound
 1393	->  R + I =< L,
 1394	    Violated = l(L),
 1395	    inc_step_cont(H,Status,ContIn,ContOut)
 1396	;   Type = t_u(U) % case 2: upperbound: R + I should always be smaller
 1397			  % than the upperbound
 1398	->  R + I >= U,
 1399	    Violated = u(U),
 1400	    dec_step_cont(H,Status,ContIn,ContOut)
 1401	;   Type = t_lu(L,U) % case 3: check both
 1402	->  At is R + I,
 1403	    (   At =< L
 1404	    ->	Violated = l(L),
 1405		inc_step_cont(H,Status,ContIn,ContOut)
 1406	    ;   At >= U
 1407	    ->	Violated = u(U),
 1408		dec_step_cont(H,Status,ContIn,ContOut)
 1409	    )
 1410	). % other types imply nonbasic variable or unbounded variable
 1411
 1412
 1413
 1414%
 1415% reconsider one element of the basis
 1416% later: lift the binds
 1417%
 1418reconsider(X) :-
 1419	rcb(X,Status,Violated),
 1420	!,
 1421	rcbl_status(Status,X,[],Binds,[],Violated),
 1422	export_binding(Binds).
 1423reconsider(_).
 1424
 1425%
 1426% Find a basis variable out of its bound or at its bound
 1427% Try to move it into whithin its bound
 1428%   a) impossible -> fail
 1429%   b) optimum at the bound -> implied value
 1430%   c) else look at the remaining basis variables
 1431%
 1432%
 1433% Idea: consider a variable V with linear equation Lin.
 1434% When a bound on a variable X of Lin gets activated, its value, multiplied
 1435% with the scalar of X, is added to the R component of Lin.
 1436% When we consider the lowerbound of V, it must be smaller than R + I, since R
 1437% contains at best the lowerbounds of the variables in Lin (but could contain
 1438% upperbounds, which are of course larger). So checking this can show the
 1439% violation of a bound of V. A similar case works for the upperbound.
 1440
 1441rcb(X,Status,Violated) :-
 1442	get_attr(X,clpqr_itf,Att),
 1443	arg(2,Att,type(Type)),
 1444	arg(4,Att,lin([I,R|H])),
 1445	(   Type = t_l(L) % case 1: lowerbound: R + I should always be larger
 1446			  % than the lowerbound
 1447	->  R + I =< L,
 1448	    Violated = l(L),
 1449	    inc_step(H,Status)
 1450	;   Type = t_u(U) % case 2: upperbound: R + I should always be smaller
 1451			  % than the upperbound
 1452	->  R + I >= U,
 1453	    Violated = u(U),
 1454	    dec_step(H,Status)
 1455	;   Type = t_lu(L,U) % case 3: check both
 1456	->  At is R + I,
 1457	    (   At =< L
 1458	    ->	Violated = l(L),
 1459		inc_step(H,Status)
 1460	    ;   At >= U
 1461	    ->	Violated = u(U),
 1462		dec_step(H,Status)
 1463	    )
 1464	). % other types imply nonbasic variable or unbounded variable
 1465
 1466% rcbl_status(Status,X,Continuation,[Bind|BindT],BindT,Violated)
 1467%
 1468%
 1469
 1470rcbl_status(optimum,X,Cont,B0,Bt,Violated) :- rcbl_opt(Violated,X,Cont,B0,Bt).
 1471rcbl_status(applied,X,Cont,B0,Bt,Violated) :- rcbl_app(Violated,X,Cont,B0,Bt).
 1472rcbl_status(unlimited(Indep,DepT),X,Cont,B0,Bt,Violated) :-
 1473	rcbl_unl(Violated,X,Cont,B0,Bt,Indep,DepT).
 1474
 1475%
 1476% Might reach optimum immediately without changing the basis,
 1477% but in general we must assume that there were pivots.
 1478% If the optimum meets the bound, we backsubstitute the implied
 1479% value, solve will call us again to check for further implied
 1480% values or unsatisfiability in the rank increased system.
 1481%
 1482rcbl_opt(l(L),X,Continuation,B0,B1) :-
 1483	get_attr(X,clpqr_itf,Att),
 1484	arg(2,Att,type(Type)),
 1485	arg(3,Att,strictness(Strict)),
 1486	arg(4,Att,lin(Lin)),
 1487	Lin = [I,R|_],
 1488	Opt is R + I,
 1489	(   L < Opt
 1490	->  narrow_u(Type,X,Opt), % { X =< Opt }
 1491	    rcbl(Continuation,B0,B1)
 1492	;   L =:= Opt,
 1493	    Strict /\ 2 =:= 0, % meets lower
 1494	    Mop is -Opt,
 1495	    normalize_scalar(Mop,MopN),
 1496	    add_linear_11(MopN,Lin,Lin1),
 1497	    Lin1 = [Inhom,_|Hom],
 1498	    (   Hom = []
 1499	    ->  rcbl(Continuation,B0,B1) % would not callback
 1500	    ;   solve(Hom,Lin1,Inhom,B0,B1)
 1501	    )
 1502	).
 1503rcbl_opt(u(U),X,Continuation,B0,B1) :-
 1504	get_attr(X,clpqr_itf,Att),
 1505	arg(2,Att,type(Type)),
 1506	arg(3,Att,strictness(Strict)),
 1507	arg(4,Att,lin(Lin)),
 1508	Lin = [I,R|_],
 1509	Opt is R + I,
 1510	(   U > Opt
 1511	->  narrow_l(Type,X,Opt), % { X >= Opt }
 1512	    rcbl(Continuation,B0,B1)
 1513	;   U =:= Opt,
 1514	    Strict /\ 1 =:= 0, % meets upper
 1515	    Mop is -Opt,
 1516	    normalize_scalar(Mop,MopN),
 1517	    add_linear_11(MopN,Lin,Lin1),
 1518	    Lin1 = [Inhom,_|Hom],
 1519	    (   Hom = []
 1520	    ->  rcbl(Continuation,B0,B1) % would not callback
 1521	    ;   solve(Hom,Lin1,Inhom,B0,B1)
 1522	    )
 1523	).
 1524
 1525%
 1526% Basis has already changed when this is called
 1527%
 1528rcbl_app(l(L),X,Continuation,B0,B1) :-
 1529	get_attr(X,clpqr_itf,Att),
 1530	arg(4,Att,lin([I,R|H])),
 1531	(   R + I > L % within bound now
 1532	->  rcbl(Continuation,B0,B1)
 1533	;   inc_step(H,Status),
 1534	    rcbl_status(Status,X,Continuation,B0,B1,l(L))
 1535	).
 1536rcbl_app(u(U),X,Continuation,B0,B1) :-
 1537	get_attr(X,clpqr_itf,Att),
 1538	arg(4,Att,lin([I,R|H])),
 1539	(   R + I < U % within bound now
 1540	->  rcbl(Continuation,B0,B1)
 1541	;   dec_step(H,Status),
 1542	    rcbl_status(Status,X,Continuation,B0,B1,u(U))
 1543	).
 1544%
 1545% This is never called for a t_lu culprit
 1546%
 1547rcbl_unl(l(L),X,Continuation,B0,B1,Indep,DepT) :-
 1548	pivot_a(X,Indep,t_L(L),DepT), % changes the basis
 1549	rcbl(Continuation,B0,B1).
 1550rcbl_unl(u(U),X,Continuation,B0,B1,Indep,DepT) :-
 1551	pivot_a(X,Indep,t_U(U),DepT), % changes the basis
 1552	rcbl(Continuation,B0,B1).
 1553
 1554% narrow_u(Type,X,U)
 1555%
 1556% Narrows down the upperbound of X (type Type) to U.
 1557% Fails if Type is not t_u(_) or t_lu(_)
 1558
 1559narrow_u(t_u(_),X,U) :-
 1560	get_attr(X,clpqr_itf,Att),
 1561	setarg(2,Att,type(t_u(U))).
 1562narrow_u(t_lu(L,_),X,U) :-
 1563	get_attr(X,clpqr_itf,Att),
 1564	setarg(2,Att,type(t_lu(L,U))).
 1565
 1566% narrow_l(Type,X,L)
 1567%
 1568% Narrows down the lowerbound of X (type Type) to L.
 1569% Fails if Type is not t_l(_) or t_lu(_)
 1570
 1571narrow_l( t_l(_),    X, L) :-
 1572	get_attr(X,clpqr_itf,Att),
 1573	setarg(2,Att,type(t_l(L))).
 1574
 1575narrow_l( t_lu(_,U), X, L) :-
 1576	get_attr(X,clpqr_itf,Att),
 1577	setarg(2,Att,type(t_lu(L,U))).
 1578
 1579% ----------------------------------- dump ------------------------------------
 1580
 1581% dump_var(Type,Var,I,H,Dump,DumpTail)
 1582%
 1583% Returns in Dump a representation of the linear constraint on variable
 1584% Var which has linear equation H + I and has type Type.
 1585
 1586dump_var(t_none,V,I,H) -->
 1587	!,
 1588	(   {
 1589		H = [l(W*K,_)],
 1590		V == W,
 1591		I =:= 0,
 1592		K =:= 1
 1593	    }
 1594	->  % indep var
 1595	    []
 1596	;   {nf2sum(H,I,Sum)},
 1597	    [V = Sum]
 1598	).
 1599dump_var(t_L(L),V,I,H) -->
 1600	!,
 1601	dump_var(t_l(L),V,I,H).
 1602% case lowerbound: V >= L or V > L
 1603% say V >= L, and V = K*V1 + ... + I, then K*V1 + ... + I >= L
 1604% and K*V1 + ... >= L-I and V1 + .../K = (L-I)/K
 1605dump_var(t_l(L),V,I,H) -->
 1606	!,
 1607	{
 1608	    H = [l(_*K,_)|_], % avoid 1 >= 0
 1609	    get_attr(V,clpqr_itf,Att),
 1610	    arg(3,Att,strictness(Strict)),
 1611	    Sm is Strict /\ 2,
 1612	    Kr is 1 rdiv K,
 1613	    Li is Kr*(L - I),
 1614	    mult_hom(H,Kr,H1),
 1615	    nf2sum(H1,0,Sum),
 1616	    (   K > 0 % K > 0
 1617	    ->	dump_strict(Sm,Sum >= Li,Sum > Li,Result)
 1618	    ;   dump_strict(Sm,Sum =< Li,Sum < Li,Result)
 1619	    )
 1620	},
 1621	[Result].
 1622dump_var(t_U(U),V,I,H) -->
 1623	!,
 1624	dump_var(t_u(U),V,I,H).
 1625dump_var(t_u(U),V,I,H) -->
 1626	!,
 1627	{
 1628	    H = [l(_*K,_)|_], % avoid 0 =< 1
 1629	    get_attr(V,clpqr_itf,Att),
 1630	    arg(3,Att,strictness(Strict)),
 1631	    Sm is Strict /\ 1,
 1632	    Kr is 1 rdiv K,
 1633	    Ui is Kr*(U-I),
 1634	    mult_hom(H,Kr,H1),
 1635	    nf2sum(H1,0.0,Sum),
 1636	    (   K > 0
 1637	    ->	dump_strict(Sm,Sum =< Ui,Sum < Ui,Result)
 1638	    ;   dump_strict(Sm,Sum >= Ui,Sum > Ui,Result)
 1639	    )
 1640	},
 1641	[Result].
 1642dump_var(t_Lu(L,U),V,I,H) -->
 1643	!,
 1644	dump_var(t_l(L),V,I,H),
 1645	dump_var(t_u(U),V,I,H).
 1646dump_var(t_lU(L,U),V,I,H) -->
 1647	!,
 1648	dump_var(t_l(L),V,I,H),
 1649	dump_var(t_u(U),V,I,H).
 1650dump_var(t_lu(L,U),V,I,H) -->
 1651	!,
 1652	dump_var(t_l(L),V,I,H),
 1653	dump_var(t_U(U),V,I,H).
 1654dump_var(T,V,I,H) --> % should not happen
 1655	[V:T:I+H].
 1656
 1657% dump_strict(FilteredStrictness,Nonstrict,Strict,Res)
 1658%
 1659% Unifies Res with either Nonstrict or Strict depending on FilteredStrictness.
 1660% FilteredStrictness is the component of strictness related to the bound: 0
 1661% means nonstrict, 1 means strict upperbound, 2 means strict lowerbound,
 1662% 3 is filtered out to either 1 or 2.
 1663
 1664dump_strict(0,Result,_,Result).
 1665dump_strict(1,_,Result,Result).
 1666dump_strict(2,_,Result,Result).
 1667
 1668% dump_nz(V,H,I,Dump,DumpTail)
 1669%
 1670% Returns in Dump a representation of the nonzero constraint of variable V
 1671% which has linear
 1672% equation H + I.
 1673
 1674dump_nz(_,H,I) -->
 1675	{
 1676	    H = [l(_*K,_)|_],
 1677	    Kr is 1 rdiv K,
 1678	    I1 is -Kr*I,
 1679	    mult_hom(H,Kr,H1),
 1680	    nf2sum(H1,0,Sum)
 1681	},
 1682	[Sum =\= I1].
 1683
 1684		 /*******************************
 1685		 *	       SANDBOX		*
 1686		 *******************************/
 1687:- multifile
 1688	sandbox:safe_primitive/1. 1689
 1690sandbox:safe_primitive(bv_q:inf(_,_)).
 1691sandbox:safe_primitive(bv_q:inf(_,_,_,_)).
 1692sandbox:safe_primitive(bv_q:sup(_,_)).
 1693sandbox:safe_primitive(bv_q:sup(_,_,_,_)).
 1694sandbox:safe_primitive(bv_q:maximize(_)).
 1695sandbox:safe_primitive(bv_q:minimize(_))