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