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