yap-6.3/CLPQR/clpqr/bv.yap

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  clp(q,r)                                         version 1.3.3 %
%                                                                 %
%  (c) Copyright 1992,1993,1994,1995                              %
%  Austrian Research Institute for Artificial Intelligence (OFAI) %
%  Schottengasse 3                                                %
%  A-1010 Vienna, Austria                                         %
%                                                                 %
%  File:   bv.pl                                                  %
%  Author: Christian Holzbaur           christian@ai.univie.ac.at %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%
%  simplex with bounded variables, ch, 93/12
%

%
% TODO: +) var/bound/state classification and maintainance
%	+) inc/dec_step: take the best?, at least find unconstrained var first
%	+) trivially implied values
%	+) avoid eval_rhs through an extra column (Coeff=Rhs)
%	+) if an optimum is encountered, record the value as bound !!!
%	+) generalized (transparent) attribute handling
%	+) coordinate reconsideration cascades
%	+) =\=
%	+) strict inequalities via =\=
%	-) decompose via nonvar test -> no symbolic constants any more ?
%	   constants complicate the nonlin solver anyway ...
%	+) join t_l,l(L), .... into t_l(L), ...
%	+) shortcuts for strict ineqs
%	-) extra types for vars with l/u bound zero
%	-) occurrence lists for indep vars (with coeffs) ???
%		each solve produces one dep var -> push
%		only complication: pivots
%	-) *incremental* REVISED simplex ?!!
%
% sicstus2.1.9.clp conversion:
%
%	-) stable ordering through extra attribute ...
%	   interpreted vs compiled yields different var order
%	   -> nasty in R (need different eps)
%
%	-) check determinism again
%
%

:- public {}/1, maximize/1, minimize/1, sup/2, inf/2, imin/2. % xref.pl

:- use_module( library(ordsets), [ord_add_element/3]).

% :- use_module( library(deterministic)).

%
% For the rhs maint. the following events are important:
%
%	-) introduction of an indep var at active bound B
%	-) narrowing of active bound
%	-) swap active bound
%	-) pivot
%

%
% a variables bound (L/U) can have the states:
%
%	-) t_none
%	-) t_l			has a lower bound (not active yet)
%	-) t_u
%	-) t_L			has an active lower bound
%	-) t_U
%	-) t_lu
%	-) t_Lu
%	-) t_lU
%

% ----------------------------------- deref ------------------------------------ %

%:- mode deref( +, -).
%
deref( Lin, Lind) :-
  split( Lin, H, I),
  normalize_scalar( I, Nonvar),
  length( H, Len),
  log_deref( Len, H, [], Restd),
  add_linear_11( Nonvar, Restd, Lind).

%:- mode log_deref( +, +, -, -).
%
log_deref( 0, Vs,	      Vs, Lin) :- !,
  arith_eval( 0, Z),
  Lin = [Z,Z].
log_deref( 1, [v(K,[X^1])|Vs],	     Vs, Lin) :- !,
  deref_var( X, Lx),
  mult_linear_factor( Lx, K, Lin).
log_deref( 2, [v(Kx,[X^1]),v(Ky,[Y^1])|Vs], Vs, Lin) :- !,
  deref_var( X, Lx),
  deref_var( Y, Ly),
  add_linear_ff( Lx, Kx, Ly, Ky, Lin).
log_deref( N, V0,	      V2, Lin) :-
  P is N >> 1,
  Q is N - P,
  log_deref( P, V0,V1, Lp),
  log_deref( Q, V1,V2, Lq),
  add_linear_11( Lp, Lq, Lin).

/*
%
% tail recursive version
%
deref( Lin, Lind) :-
  split( Lin, H, I),
  normalize_scalar( I, Nonvar),
  lin_deref( H, Nonvar, Lind).

log_deref( _, Lin, [], Res) :-			% called from nf.pl
  arith_eval( 0, Z),
  lin_deref( Lin, [Z,Z], Res).

lin_deref( [],		    Ld, Ld).
lin_deref( [v(K,[X^1])|Vs], Li, Lo) :-
  deref_var( X, Lx),
  add_linear_f1( Lx, K, Li, Lii),
  lin_deref( Vs, Lii, Lo).
*/

%
% If we see a nonvar here, this is a fault
%

deref_var( X, Lin) :-
  get_atts( X, lin(Lin)), !.

deref_var( X, Lin) :-		% create a linear var
  arith_eval( 0, Z),
  arith_eval( 1, One),
  Lin = [Z,Z,X*One],
  put_atts( X, [order(_),lin(Lin),type(t_none),strictness(2'00)]).

var_with_def_assign( Var, Lin) :-
  decompose( Lin, Hom, _, I),
  ( Hom = [],					% X=k
      Var = I
  ; Hom = [V*K|Cs],
      ( Cs = [],
	arith_eval(K=:=1),
	arith_eval(I=:=0) ->			% X=Y
	  Var = V
      ; 					% general case
	  var_with_def_intern( t_none, Var, Lin, 2'00)
      )
  ).

var_with_def_intern( Type, Var, Lin, Strict) :-
  put_atts( Var, [order(_),lin(Lin),type(Type),strictness(Strict)]),
  decompose( Lin, Hom, _, _),
  get_or_add_class( Var, Class),
  same_class( Hom, Class).

var_intern( Type, Var, Strict) :-
  arith_eval( 0, Z),
  arith_eval( 1, One),
  Lin = [Z,Z,Var*One],
  put_atts( Var, [order(_),lin(Lin),type(Type),strictness(Strict)]),
  get_or_add_class( Var, _Class).

% ------------------------------------------------------------------------------

%
% [V-Binding]*
% Only place where the linear solver binds variables
%
export_binding( []).
export_binding( [X-Y|Gs]) :-
  export_binding( Y, X),
  export_binding( Gs).

%
% numerical stabilizer, clp(r) only
%
export_binding( Y, X) :- var(Y), !, Y=X. %vsc: added cut here (01/06/06)
export_binding( Y, X) :- nonvar(Y),
  ( arith_eval( Y=:=0) ->
      arith_eval( 0, X)
  ;
      Y = X
  ).

'solve_='( Nf) :-
  deref( Nf, Nfd),
  solve( Nfd).

'solve_=\\='( Nf) :- % vsc
  deref( Nf, Lind),
  decompose( Lind, Hom, _, Inhom),
  ( Hom = [],	 arith_eval( Inhom =\= 0)
  ; Hom = [_|_], var_with_def_intern( t_none, Nz, Lind, 2'00),
		 put_atts( Nz, nonzero)
  ).

'solve_<'( Nf) :-
  split( Nf, H, I),
  ineq( H, I, Nf, strict).

'solve_=<'( Nf) :-
  split( Nf, H, I),
  ineq( H, I, Nf, nonstrict).

maximize( Term) :-
	minimize( -Term).

%
% This is NOT coded as minimize(Expr) :- inf(Expr,Expr).
%
% because the new version of inf/2 only visits
% the vertex where the infimum is assumed and returns
% to the 'current' vertex via backtracking.
% The rationale behind this construction is to eliminate
% all garbage in the solver data structures produced by
% the pivots on the way to the extremal point caused by
% {inf,sup}/{2,4}.
%
% If we are after the infimum/supremum for minimizing/maximizing,
% this strategy may have adverse effects on performance because
% the simplex algorithm is forced to re-discover the
% extremal vertex through the equation {Inf =:= Expr}.
%
% Thus the extra code for {minimize,maximize}/1.
%
% In case someone comes up with an example where
%
%   inf(Expr,Expr)
%
% outperforms the provided formulation for minimize - so be it.
% Both forms are available to the user.
%
minimize( Term) :-
	wait_linear( Term, Nf, minimize_lin(Nf)).

minimize_lin( Lin) :-
	deref( Lin, Lind),
	var_with_def_intern( t_none, Dep, Lind, 2'00),
	determine_active_dec( Lind),
	iterate_dec( Dep, Inf),
	{ Dep =:= Inf }.

sup( Expression, Sup) :-
	sup( Expression, Sup, [], []).

sup( Expression, Sup, Vector, Vertex) :-
	inf( -Expression, -Sup, Vector, Vertex).

inf( Expression, Inf) :-
	inf( Expression, Inf, [], []).

inf( Expression, Inf, Vector, Vertex) :-
	wait_linear( Expression, Nf, inf_lin(Nf,Inf,Vector,Vertex)).

inf_lin( Lin, _, Vector, _) :-
	deref( Lin, Lind),
	var_with_def_intern( t_none, Dep, Lind, 2'00),
	determine_active_dec( Lind),
	iterate_dec( Dep, Inf),
	vertex_value( Vector, Values),
	bb_put( inf, [Inf|Values]),
	fail.
inf_lin( _, Infimum, _, Vertex) :-
	bb_delete( inf, L),
	assign( [Infimum|Vertex], L).

assign( [],     []).
assign( [X|Xs], [Y|Ys]) :-
	{X =:= Y},				  % more defensive/expressive than X=Y
	assign( Xs, Ys).

% --------------------------------- optimization ------------------------------- %
%
% The _sn(S) =< 0 row might be temporarily infeasible.
% We use reconsider/1 to fix this.
%
%   s(S) e [_,0] = d +xi ... -xj, Rhs > 0 so we want to decrease s(S)
%
%   positive xi would have to be moved towards their lower bound,
%   negative xj would have to be moved towards their upper bound,
%
%   the row s(S) does not limit the lower bound of xi
%   the row s(S) does not limit the upper bound of xj
%
%   a) if some other row R is limiting xk, we pivot(R,xk),
%      s(S) will decrease and get more feasible until (b)
%   b) if there is no limiting row for some xi: we pivot(s(S),xi)
%					    xj: we pivot(s(S),xj)
%      which cures the infeasibility in one step
%


%
% fails if Status = unlimited/2
%
iterate_dec( OptVar, Opt) :-
  get_atts( OptVar, lin(Lin)),
  decompose( Lin, H, R, I),

  % arith_eval( R+I, Now), print(min(Now)), nl,

  % dec_step_best( H, Status),
  %vsc: added -> (01/06/06)
  dec_step( H, Status),
  ( Status = applied -> iterate_dec( OptVar, Opt)
  ; Status = optimum -> arith_eval( R+I, Opt)
  ).

iterate_inc( OptVar, Opt) :-
  get_atts( OptVar, lin(Lin)),
  decompose( Lin, H, R, I),
  inc_step( H, Status),
  %vsc: added -> (01/06/06)
  ( Status = applied -> iterate_inc( OptVar, Opt)
  ; Status = optimum -> arith_eval( R+I, Opt)
  ).

%
% Status = {optimum,unlimited(Indep,DepT),applied}
% If Status = optimum, the tables have not been changed at all.
% Searches left to right, does not try to find the 'best' pivot
% Therefore we might discover unboundedness only after a few pivots
%
dec_step( [],	    optimum).
dec_step( [V*K|Vs], Status) :-
  get_atts( V, type(W)),
  %vsc: added -> (01/06/06)
  ( W = t_U(U) ->
      ( arith_eval( K > 0) ->
	 ( lb( V, Vub-Vb-_) ->
	     Status = applied,
	     pivot_a(Vub,V,Vb,t_u(U))
	 ;
	     Status = unlimited(V,t_u(U))
	 )
      ;
	 dec_step( Vs, Status)
      )
  ; W = t_lU(L,U) ->
      ( arith_eval( K > 0) ->
	 Status = applied,
	 arith_eval( L-U, Init),
	 basis( V, Deps),
	 lb( Deps, V, V-t_Lu(L,U)-Init, Vub-Vb-_),
	 pivot_b(Vub,V,Vb,t_lu(L,U))
      ;
	 dec_step( Vs, Status)
      )
  ; W = t_L(L) ->
      ( arith_eval( K < 0) ->
	 ( ub( V, Vub-Vb-_) ->
	     Status = applied,
	     pivot_a(Vub,V,Vb,t_l(L))
	 ;
	     Status = unlimited(V,t_l(L))
	 )
      ;
	 dec_step( Vs, Status)
      )
  ; W = t_Lu(L,U) ->
      ( arith_eval( K < 0) ->
	 Status = applied,
	 arith_eval( U-L, Init),
	 basis( V, Deps),
	 ub( Deps, V, V-t_lU(L,U)-Init, Vub-Vb-_),
	 pivot_b(Vub,V,Vb,t_lu(L,U))
      ;
	 dec_step( Vs, Status)
      )
  ; W = t_none ->
      Status = unlimited(V,t_none)
  ).

inc_step( [],	    optimum).
inc_step( [V*K|Vs], Status) :-
  get_atts( V, type(W)),
  %vsc: added -> (01/06/06)
  ( W = t_U(U) ->
      ( arith_eval( K < 0) ->
	 ( lb( V, Vub-Vb-_) ->
	     Status = applied,
	     pivot_a(Vub,V,Vb,t_u(U))
	 ;
	     Status = unlimited(V,t_u(U))
	 )
      ;
	 inc_step( Vs, Status)
      )
  ; W = t_lU(L,U) ->
      ( arith_eval( K < 0) ->
	 Status = applied,
	 arith_eval( L-U, Init),
	 basis( V, Deps),
	 lb( Deps, V, V-t_Lu(L,U)-Init, Vub-Vb-_),
	 pivot_b(Vub,V,Vb,t_lu(L,U))
      ;
	 inc_step( Vs, Status)
      )
  ; W = t_L(L) ->
      ( arith_eval( K > 0) ->
	 ( ub( V, Vub-Vb-_) ->
	     Status = applied,
	     pivot_a(Vub,V,Vb,t_l(L))
	 ;
	     Status = unlimited(V,t_l(L))
	 )
      ;
	 inc_step( Vs, Status)
      )
  ; W = t_Lu(L,U) ->
      ( arith_eval( K > 0) ->
	 Status = applied,
	 arith_eval( U-L, Init),
	 basis( V, Deps),
	 ub( Deps, V, V-t_lU(L,U)-Init, Vub-Vb-_),
	 pivot_b(Vub,V,Vb,t_lu(L,U))
      ;
	 inc_step( Vs, Status)
      )
  ; W = t_none ->
      Status = unlimited(V,t_none)
  ).

% ------------------------------ best first heuristic -------------------------- %
%
% A replacement for dec_step/2 that uses a local best first heuristic.
%
%

dec_step_best( H, Status) :-
  dec_eval( H, E),
  ( E = unlimited(_,_),
     Status = E
  ; E = [],
     Status = optimum
  ; E = [_|_],
     Status = applied,
     keysort( E, [_-Best|_]),
     ( Best = pivot_a(Vub,V,Vb,Wd), pivot_a(Vub,V,Vb,Wd)
     ; Best = pivot_b(Vub,V,Vb,Wd), pivot_b(Vub,V,Vb,Wd)
     )
  ).

dec_eval( [],	    []).
dec_eval( [V*K|Vs], Res) :-
  get_atts( V, type(W)),
  ( W = t_U(U),
      ( arith_eval( K > 0) ->
	 ( lb( V, Vub-Vb-Limit) ->
	     arith_eval( float(Limit*K), Delta),
	     Res = [Delta-pivot_a(Vub,V,Vb,t_u(U)) | Tail],
	     dec_eval( Vs, Tail)
	 ;
	     Res = unlimited(V,t_u(U))
	 )
      ;
	 dec_eval( Vs, Res)
      )
  ; W = t_lU(L,U),
      ( arith_eval( K > 0) ->
	 arith_eval( L-U, Init),
	 basis( V, Deps),
	 lb( Deps, V, V-t_Lu(L,U)-Init, Vub-Vb-Limit),
	 arith_eval( float(Limit*K), Delta),
	 Res = [Delta-pivot_b(Vub,V,Vb,t_lu(L,U)) | Tail],
	 dec_eval( Vs, Tail)
      ;
	 dec_eval( Vs, Res)
      )
  ; W = t_L(L),
      ( arith_eval( K < 0) ->
	 ( ub( V, Vub-Vb-Limit) ->
	     arith_eval( float(Limit*K), Delta),
	     Res = [Delta-pivot_a(Vub,V,Vb,t_l(L)) | Tail],
	     dec_eval( Vs, Tail)
	 ;
	     Res = unlimited(V,t_l(L))
	 )
      ;
	 dec_eval( Vs, Res)
      )
  ; W = t_Lu(L,U),
      ( arith_eval( K < 0) ->
	 arith_eval( U-L, Init),
	 basis( V, Deps),
	 ub( Deps, V, V-t_lU(L,U)-Init, Vub-Vb-Limit),
	 arith_eval( float(Limit*K), Delta),
	 Res = [Delta-pivot_b(Vub,V,Vb,t_lu(L,U)) | Tail],
	 dec_eval( Vs, Tail)
      ;
	 dec_eval( Vs, Res)
      )
  ; W = t_none,
      Res = unlimited(V,t_none)
  ).

% ------------------------- find the most constraining row --------------------- %
%
% The code for the lower and the upper bound are dual versions of each other.
% The only difference is in the orientation of the comparisons.
% Indeps are ruled out by their types.
% If there is no bound, this fails.
%
% *** The actual lb and ub on an indep variable X are [lu]b + b(X), where b(X)
% is the value of the active bound.
%
% Nota bene: We must NOT consider infeasible rows as candidates to
%	     leave the basis!
%

ub( X, Ub) :-
  basis( X, Deps),
  ub_first( Deps, X, Ub).

%:- mode ub_first( +, ?, -).
%
ub_first( [Dep|Deps], X, Tightest) :-
  ( get_atts( Dep, [lin(Lin),type(Type)]),
    ub_inner( Type, X, Lin, W, Ub),
    arith_eval( Ub >= 0) ->
      ub( Deps, X, Dep-W-Ub, Tightest)
  ;
      ub_first( Deps, X, Tightest)
  ).

%
% Invariant: Ub >= 0 and decreasing
%
%:- mode ub( +, ?, +, -).
%
ub( [], 	_, T0,T0).
ub( [Dep|Deps], X, T0,T1) :-
  ( get_atts( Dep, [lin(Lin),type(Type)]),
    ub_inner( Type, X, Lin, W, Ub),
    T0 = _-Ubb,
    arith_eval( Ub < Ubb),
    arith_eval( Ub >= 0) ->			% rare failure
      ub( Deps, X, Dep-W-Ub,T1)
  ;
      ub( Deps, X, T0,T1)
  ).

lb( X, Lb) :-
  basis( X, Deps),
  lb_first( Deps, X, Lb).

%:- mode lb_first( +, ?, -).
%
lb_first( [Dep|Deps], X, Tightest) :-
  ( get_atts( Dep, [lin(Lin),type(Type)]),
    lb_inner( Type, X, Lin, W, Lb),
    arith_eval( Lb =< 0) ->
      lb( Deps, X, Dep-W-Lb, Tightest)
  ;
      lb_first( Deps, X, Tightest)
  ).

%
% Invariant: Lb =< 0 and increasing
%
%:- mode lb( +, ?, +, -).
%
lb( [], 	_, T0,T0).
lb( [Dep|Deps], X, T0,T1) :-
  ( get_atts( Dep, [lin(Lin),type(Type)]),
    lb_inner( Type, X, Lin, W, Lb),
    T0 = _-Lbb,
    arith_eval( Lb > Lbb),
    arith_eval( Lb =< 0) ->			% rare failure
      lb( Deps, X, Dep-W-Lb,T1)
  ;
      lb( Deps, X, T0,T1)
  ).

%
% Lb =< 0 for feasible rows
%
%:- mode lb_inner( +, ?, +, -, -).
%
lb_inner( t_l(L), X, Lin, t_L(L), Lb) :-
  nf_rhs_x( Lin, X, Rhs, K),
  arith_eval( K > 0),
  arith_eval( (L-Rhs)/K, Lb).
lb_inner( t_u(U), X, Lin, t_U(U), Lb) :-
  nf_rhs_x( Lin, X, Rhs, K),
  arith_eval( K < 0),
  arith_eval( (U-Rhs)/K, Lb).
lb_inner( t_lu(L,U), X, Lin, W, Lb) :-
  nf_rhs_x( Lin, X, Rhs, K),
  case_signum( K,
   (
    W = t_lU(L,U),
    arith_eval( (U-Rhs)/K, Lb)
   ),
    fail,
   (
    W = t_Lu(L,U),
    arith_eval( (L-Rhs)/K, Lb)
   )).

%
% Ub >= 0 for feasible rows
%
%:- mode ub_inner( +, ?, +, -, -).
%
ub_inner( t_l(L), X, Lin, t_L(L), Ub) :-
  nf_rhs_x( Lin, X, Rhs, K),
  arith_eval( K < 0),
  arith_eval( (L-Rhs)/K, Ub).
ub_inner( t_u(U), X, Lin, t_U(U), Ub) :-
  nf_rhs_x( Lin, X, Rhs, K),
  arith_eval( K > 0),
  arith_eval( (U-Rhs)/K, Ub).
ub_inner( t_lu(L,U), X, Lin, W, Ub) :-
  nf_rhs_x( Lin, X, Rhs, K),
  case_signum( K,
   (
    W = t_Lu(L,U),
    arith_eval( (L-Rhs)/K, Ub)
   ),
    fail,
   (
    W = t_lU(L,U),
    arith_eval( (U-Rhs)/K, Ub)
   )).

% ---------------------------------- equations --------------------------------- %
%
% backsubstitution will not make the system infeasible, if the bounds on the indep
% vars are obeyed, but some implied values might pop up in rows where X occurs
%	-) special case X=Y during bs -> get rid of dependend var(s), alias
%

solve( Lin) :-
  decompose( Lin, H, _, I),
  solve( H, Lin, I, Bindings, []),
  export_binding( Bindings).

solve( [], _,  I, Bind0,Bind0) :-
  arith_eval( I=:=0).					% redundant or trivially unsat
%vsc: changed to list in head  (01/06/06)
solve( [HHd|HTl], Lin, _, Bind0,BindT) :-
  %
  % [] is an empty ord_set, anything will be preferred
  % over 9-9
  %
  sd( [HHd|HTl], [],ClassesUniq, 9-9-0,Category-Selected-_, NV,NVT),

  isolate( Selected, Lin, Lin1),

  %vsc: added -> (01/06/06)
  ( Category = 1 ->
      put_atts( Selected, lin(Lin1)),
      decompose( Lin1, Hom, _, Inhom),
      bs_collect_binding( Hom, Selected, Inhom, Bind0,BindT),
      eq_classes( NV, NVT, ClassesUniq)
  ; Category = 2 ->
      get_atts( Selected, class(NewC)),
      class_allvars( NewC, Deps),
      ( ClassesUniq = [_] ->				% rank increasing
	  bs_collect_bindings( Deps, Selected, Lin1, Bind0,BindT)
      ;
	  Bind0 = BindT,
	  bs( Deps, Selected, Lin1)
      ),
      eq_classes( NV, NVT, ClassesUniq)
  ; Category = 3 ->
      put_atts( Selected, lin(Lin1)),
      get_atts( Selected, type(Type)),
      deactivate_bound( Type, Selected),
      eq_classes( NV, NVT, ClassesUniq),
      basis_add( Selected, Basis),
      undet_active( Lin1),
      decompose( Lin1, Hom, _, Inhom),
      bs_collect_binding( Hom, Selected, Inhom, Bind0,Bind1),
      rcbl( Basis, Bind1,BindT)
  ; Category = 4 ->
      get_atts( Selected, [type(Type),class(NewC)]),
      class_allvars( NewC, Deps),
      ( ClassesUniq = [_] ->				% rank increasing
	  bs_collect_bindings( Deps, Selected, Lin1, Bind0,Bind1)
      ;
	  Bind0 = Bind1,
	  bs( Deps, Selected, Lin1)
      ),
      deactivate_bound( Type, Selected),
      basis_add( Selected, Basis),
      % eq_classes( NV, NVT, ClassesUniq),		%  4 -> var(NV)
      equate( ClassesUniq, _),
      undet_active( Lin1),
      rcbl( Basis, Bind1,BindT)
  ).

%
% Much like solve, but we solve for a particular variable of type
% t_none
%
solve_x( Lin, X) :-
  decompose( Lin, H, _, I),
  solve_x( H, Lin, I, X, Bindings, []),
  export_binding( Bindings).

solve_x( [], _,  I, _,	      Bind0,Bind0) :-
  arith_eval( I=:=0).					% redundant or trivially unsat
solve_x( H, Lin, _, Selected, Bind0,BindT) :-
  H = [_|_],						% indexing
  sd( H, [],ClassesUniq, 9-9-0,_, NV,NVT),

  isolate( Selected, Lin, Lin1),

  ( get_atts( Selected, class(NewC)) ->
      class_allvars( NewC, Deps),
      ( ClassesUniq = [_] ->				% rank increasing
	  bs_collect_bindings( Deps, Selected, Lin1, Bind0,BindT)
      ;
	  Bind0 = BindT,
	  bs( Deps, Selected, Lin1)
      ),
      eq_classes( NV, NVT, ClassesUniq)
  ;
      put_atts( Selected, lin(Lin1)),
      decompose( Lin1, Hom, _, Inhom),
      bs_collect_binding( Hom, Selected, Inhom, Bind0,BindT),
      eq_classes( NV, NVT, ClassesUniq)
  ).



sd( [],       Class0,Class0, Preference0,Preference0, NV0,NV0).
sd( [X*K|Xs], Class0,ClassN, Preference0,PreferenceN, NV0,NVt) :-
  ( get_atts( X, class(Xc)) ->						% old
	NV0 = NV1,
	ord_add_element( Class0, Xc, Class1),
	( get_atts( X, type(t_none)) ->
	    preference( Preference0, 2-X-K, Preference1)
	;
	    preference( Preference0, 4-X-K, Preference1)
	)
  ;									% new
	Class1 = Class0,
	'C'( NV0, X, NV1),
	( get_atts( X, type(t_none)) ->
	    preference( Preference0, 1-X-K, Preference1)
	;
	    preference( Preference0, 3-X-K, Preference1)
	)
  ),
  sd( Xs, Class1,ClassN, Preference1,PreferenceN, NV1,NVt).

%
% A is best sofar, B is current
%
preference( A, B, Pref) :-
  A = Px-_-_,
  B = Py-_-_,
  compare( Rel, Px, Py),
  %vsc: added -> (01/06/06)
  ( Rel = = -> Pref = B
	     % ( arith_eval(abs(Ka)=<abs(Kb)) -> Pref=A ; Pref=B )
  ; Rel = < -> Pref = A
  ; Rel = > -> Pref = B
  ).

%
% equate after attach_class because other classes may contribute
% nonvars and will bind the tail of NV
%
eq_classes( NV, _,   Cs) :- var( NV), !,
  equate( Cs, _).
eq_classes( NV, NVT, Cs) :-
  class_new( Su, NV,NVT, []),
  attach_class( NV, Su),
  equate( Cs, Su).

equate( [],	_).
equate( [X|Xs], X) :- equate( Xs, X).

%
% assert: none of the Vars has a class attribute yet
%
attach_class( Xs,     _) :- var( Xs), !.
attach_class( [X|Xs], Class) :-
  put_atts( X, class(Class)),
  attach_class( Xs, Class).

/**
unconstrained( [X*K|Xs], Uc,Kuc, Rest) :-
  ( get_atts( X, type(t_none)) ->
      Uc = X,
      Kuc = K,
      Rest = Xs
  ;
      Rest = [X*K|Tail],
      unconstrained( Xs, Uc,Kuc, Tail)
  ).
**/
/**/
unconstrained( Lin, Uc,Kuc, Rest) :-
  decompose( Lin, H, _, _),
  sd( H, [],_, 9-9-0,Category-Uc-_, _,_),
  Category =< 2,
  delete_factor( Uc, Lin, Rest, Kuc).
/**/

%
% point the vars in Lin into the same equivalence class
% maybe join some global data
%
same_class( [],       _).
same_class( [X*_|Xs], Class) :-
  get_or_add_class( X, Class),
  same_class( Xs, Class).

get_or_add_class( X, Class) :-
  get_atts( X, class(ClassX)),
  !,
  ClassX = Class.				% explicit =/2 because of cut
get_or_add_class( X, Class) :-
  put_atts( X, class(Class)),
  class_new( Class, [X|Tail],Tail, []). 	% initial class atts

allvars( X, Allvars) :-
  get_atts( X, class(C)),
  class_allvars( C, Allvars).

deactivate_bound( t_l(_),    _).
deactivate_bound( t_u(_),    _).
deactivate_bound( t_lu(_,_), _).
deactivate_bound( t_L(L),    X) :- put_atts( X, type(t_l(L))).
deactivate_bound( t_Lu(L,U), X) :- put_atts( X, type(t_lu(L,U))).
deactivate_bound( t_U(U),    X) :- put_atts( X, type(t_u(U))).
deactivate_bound( t_lU(L,U), X) :- put_atts( X, type(t_lu(L,U))).

intro_at( X, Value, Type) :-
  put_atts( X, type(Type)),
  ( arith_eval( Value =:= 0) ->
      true
  ;
      backsubst_delta( X, Value)
  ).


%
% The choice t_lu -> t_Lu is arbitrary
%
undet_active( Lin) :-
  decompose( Lin, Lin1, _, _),
  undet_active_h( Lin1).

undet_active_h( []).
undet_active_h( [X*_|Xs]) :-
  get_atts( X, type(Type)),
  undet_active( Type, X),
  undet_active_h( Xs).

undet_active( t_none,	 _).			% type_activity
undet_active( t_L(_),	 _).
undet_active( t_Lu(_,_), _).
undet_active( t_U(_),	 _).
undet_active( t_lU(_,_), _).
undet_active( t_l(L),	 X) :- intro_at( X, L, t_L(L)).
undet_active( t_u(U),	 X) :- intro_at( X, U, t_U(U)).
undet_active( t_lu(L,U), X) :- intro_at( X, L, t_Lu(L,U)).

determine_active_dec( Lin) :-
  decompose( Lin, Lin1, _, _),
  arith_eval( -1, Mone),
  determine_active( Lin1, Mone).

determine_active_inc( Lin) :-
  decompose( Lin, Lin1, _, _),
  arith_eval( 1, One),
  determine_active( Lin1, One).

determine_active( [],	    _).
determine_active( [X*K|Xs], S) :-
  get_atts( X, type(Type)),
  determine_active( Type, X, K, S),
  determine_active( Xs, S).

determine_active( t_L(_),    _, _, _).
determine_active( t_Lu(_,_), _, _, _).
determine_active( t_U(_),    _, _, _).
determine_active( t_lU(_,_), _, _, _).
determine_active( t_l(L),    X, _, _) :- intro_at( X, L, t_L(L)).
determine_active( t_u(U),    X, _, _) :- intro_at( X, U, t_U(U)).
determine_active( t_lu(L,U), X, K, S) :-
  case_signum( K*S,
    intro_at( X, L, t_Lu(L,U)),
    fail,
    intro_at( X, U, t_lU(L,U))).

%
% Careful when an indep turns into t_none !!!
%
detach_bounds( V) :-
  get_atts( V, lin(Lin)),
  put_atts( V, [type(t_none),strictness(2'00)]),
  ( indep( Lin, V) ->
      ( ub( V, Vub-Vb-_) ->			% exchange against thightest
	  basis_drop( Vub),
	  pivot( Vub, V, Vb)
      ; lb( V, Vlb-Vb-_) ->
	  basis_drop( Vlb),
	  pivot( Vlb, V, Vb)
      ;
	  true
      )
  ;
      basis_drop( V)
  ).

% ----------------------------- manipulate the basis --------------------------- %

basis_drop( X) :-
  get_atts( X, class(Cv)),
  class_basis_drop( Cv, X).

basis( X, Basis) :-
  get_atts( X, class(Cv)),
  class_basis( Cv, Basis).

basis_add( X, NewBasis) :-
  get_atts( X, class(Cv)),
  class_basis_add( Cv, X, NewBasis).

basis_pivot( Leave, Enter) :-
  get_atts( Leave, class(Cv)),
  class_basis_pivot( Cv, Enter, Leave).

% ----------------------------------- pivot ------------------------------------ %

%
% Pivot ignoring rhs and active states
%
pivot( Dep, Indep) :-
  get_atts( Dep, lin(H)),
  delete_factor( Indep, H, H0, Coeff),
  arith_eval( -1/Coeff, K),
  arith_eval( -1, Mone),
  arith_eval(  0, Z),
  add_linear_ff( H0, K, [Z,Z,Dep*Mone], K, Lin),
  backsubst( Indep, Lin).


pivot_a( Dep, Indep, Vb,Wd) :-
  basis_pivot( Dep, Indep),
  pivot( Dep, Indep, Vb),
  put_atts( Indep, type(Wd)).

pivot_b( Vub, V, Vb, Wd) :-
  ( Vub == V ->
      put_atts( V, type(Vb)),
      pivot_b_delta( Vb, Delta),		% nonzero(Delta)
      backsubst_delta( V, Delta)
  ;
      pivot_a( Vub, V, Vb,Wd)
  ).

pivot_b_delta( t_Lu(L,U), Delta) :- arith_eval( L-U, Delta).
pivot_b_delta( t_lU(L,U), Delta) :- arith_eval( U-L, Delta).

select_active_bound( t_L(L),	L).
select_active_bound( t_Lu(L,_), L).
select_active_bound( t_U(U),	U).
select_active_bound( t_lU(_,U), U).
select_active_bound( t_none,	Z) :- arith_eval( 0, Z).
%
% for project.pl
%
select_active_bound( t_l(_),	Z) :- arith_eval( 0, Z).
select_active_bound( t_u(_),	Z) :- arith_eval( 0, Z).
select_active_bound( t_lu(_,_), Z) :- arith_eval( 0, Z).


%
% Pivot taking care of rhs and active states
%
pivot( Dep, Indep, IndAct) :-
  get_atts( Dep, lin(H)),
  put_atts( Dep, type(IndAct)),
  select_active_bound( IndAct, Abv),			% Dep or Indep
  delete_factor( Indep, H, H0, Coeff),
  arith_eval( -1/Coeff, K),
  arith_eval( 0, Z),
  arith_eval( -1, Mone),
  arith_eval( -Abv, Abvm),
  add_linear_ff( H0, K, [Z,Abvm,Dep*Mone], K, Lin),
  backsubst( Indep, Lin).

backsubst_delta( X, Delta) :-
  arith_eval( 1, One),
  arith_eval( 0, Z),
  backsubst( X, [Z,Delta,X*One]).

backsubst( X, Lin) :-
  allvars( X, Allvars),
  bs( Allvars, X, Lin).
%
% valid if nothing will go ground
%
bs( Xs,     _, _) :- var( Xs), !.
bs( [X|Xs], V, Lin) :-
  ( get_atts( X, lin(LinX)),
    nf_substitute( V, Lin, LinX, LinX1) ->
      put_atts( X, lin(LinX1)),
      bs( Xs, V, Lin)
  ;
      bs( Xs, V, Lin)
  ).


%
% rank increasing backsubstitution
%
bs_collect_bindings( Xs,     _, _,   Bind0,BindT) :- var( Xs), !, Bind0=BindT.
bs_collect_bindings( [X|Xs], V, Lin, Bind0,BindT) :-
  ( get_atts( X, lin(LinX)),
    nf_substitute( V, Lin, LinX, LinX1) ->
      put_atts( X, lin(LinX1)),
      decompose( LinX1, Hom, _, Inhom),
      bs_collect_binding( Hom, X, Inhom, Bind0,Bind1),
      bs_collect_bindings( Xs, V, Lin, Bind1,BindT)
  ;
      bs_collect_bindings( Xs, V, Lin, Bind0,BindT)
  ).

%
% The first clause exports bindings,
% the second (no longer) aliasings
%
bs_collect_binding( [],       X, Inhom) --> [ X-Inhom ].
bs_collect_binding( [_|_],    _, _)     --> [].
/*
bs_collect_binding( [Y*K|Ys], X, Inhom) -->
  ( { Ys = [],
      Y \== X,
      arith_eval( K=:=1),
      arith_eval( Inhom=:=0)
    } ->
      [ X-Y ]
  ;
      []
  ).
*/

%
% reconsider the basis
%
rcbl( [],		Bind0,Bind0).
rcbl( [X|Continuation], Bind0,BindT) :-
  ( rcb( X, Status, Violated) ->			% have a culprit
      rcbl_status( Status, X, Continuation, Bind0,BindT, Violated)
  ;
      rcbl( Continuation, Bind0,BindT)
  ).

%
% reconsider one element of the basis
% later: lift the binds
%
reconsider( X) :-
  rcb( X, Status, Violated),
  !,
  rcbl_status( Status, X, [], Binds,[], Violated),
  export_binding( Binds).
reconsider( _).

%
% Find a basis variable out of its bound or at its bound
% Try to move it into whithin its bound
%   a) impossible -> fail
%   b) optimum at the bound -> implied value
%   c) else look at the remaining basis variables
%
rcb( X, Status, Violated) :-
  get_atts( X, [lin(Lin),type(Type)]),
  decompose( Lin, H, R, I),
  ( Type = t_l(L),
		 arith_eval( R+I =< L),
		 Violated = l(L),
		 inc_step( H, Status)

  ; Type = t_u(U),
		 arith_eval( R+I >= U),
		 Violated = u(U),
		 dec_step( H, Status)

  ; Type = t_lu(L,U),
		 arith_eval( R+I, At),
		 (
		   arith_eval( At =< L),
		   Violated = l(L),
		   inc_step( H, Status)
		 ;
		   arith_eval( At >= U),
		   Violated = u(U),
		   dec_step( H, Status)
		 )
  %
  % don't care for other types
  %
  ).

rcbl_status( optimum,		    X, Cont, B0,Bt, Violated) :- rcbl_opt( Violated, X, Cont, B0,Bt).
rcbl_status( applied,		    X, Cont, B0,Bt, Violated) :- rcbl_app( Violated, X, Cont, B0,Bt).
rcbl_status( unlimited(Indep,DepT), X, Cont, B0,Bt, Violated) :- rcbl_unl( Violated, X, Cont, B0,Bt, Indep, DepT).

%
% Might reach optimum immediately without changing the basis,
% but in general we must assume that there were pivots.
% If the optimum meets the bound, we backsubstitute the implied
% value, solve will call us again to check for further implied
% values or unsatisfiability in the rank increased system.
%
rcbl_opt( l(L), X, Continuation, B0,B1) :-
  get_atts( X, [lin(Lin),strictness(Strict),type(Type)]),
  decompose( Lin, _, R, I),
  arith_eval( R+I, Opt),
  case_signum( L-Opt,
   (
    narrow_u( Type, X, Opt),			% { X =< Opt }
    rcbl( Continuation, B0,B1)
   ),
   (
    Strict /\ 2'10 =:= 0,			% meets lower
    arith_eval( -Opt, Mop),
    normalize_scalar( Mop, MopN),
    add_linear_11( MopN, Lin, Lin1),
    decompose( Lin1, Hom, _, Inhom),
  %vsc: added -> (01/06/06)
    ( Hom = [] ->    rcbl( Continuation, B0,B1)	% would not callback
    ; Hom = [_|_] -> solve( Hom, Lin1, Inhom, B0,B1)
    )
   ),
    fail
   ).
rcbl_opt( u(U), X, Continuation, B0,B1) :-
  get_atts( X, [lin(Lin),strictness(Strict),type(Type)]),
  decompose( Lin, _, R, I),
  arith_eval( R+I, Opt),
  case_signum( U-Opt,
    fail,
   (
    Strict /\ 2'01 =:= 0,			% meets upper
    arith_eval( -Opt, Mop),
    normalize_scalar( Mop, MopN),
    add_linear_11( MopN, Lin, Lin1),
    decompose( Lin1, Hom, _, Inhom),
  %vsc: added -> (01/06/06)
    ( Hom = []    ->  rcbl( Continuation, B0,B1)	% would not callback
    ; Hom = [_|_] -> solve( Hom, Lin1, Inhom, B0,B1)
    )
   ),
   (
    narrow_l( Type, X, Opt),			% { X >= Opt }
    rcbl( Continuation, B0,B1)
   )).

%
% Basis has already changed when this is called
%
rcbl_app( l(L), X, Continuation, B0,B1) :-
  get_atts( X, lin(Lin)),
  decompose( Lin, H, R, I),
  ( arith_eval( R+I > L) ->			% within bound now
      rcbl( Continuation, B0,B1)
  ;
      % arith_eval( R+I, Val), print( rcbl_app(X:L:Val)), nl,
      inc_step( H, Status),
      rcbl_status( Status, X, Continuation, B0,B1, l(L))
  ).
rcbl_app( u(U), X, Continuation, B0,B1) :-
  get_atts( X, lin(Lin)),
  decompose( Lin, H, R, I),
  ( arith_eval( R+I < U) ->			% within bound now
      rcbl( Continuation, B0,B1)
  ;
      dec_step( H, Status),
      rcbl_status( Status, X, Continuation, B0,B1, u(U))
  ).

%
% This is never called for a t_lu culprit
%
rcbl_unl( l(L), X, Continuation, B0,B1, Indep, DepT) :-
  pivot_a( X, Indep, t_L(L), DepT),		% changes the basis
  rcbl( Continuation, B0,B1).
rcbl_unl( u(U), X, Continuation, B0,B1, Indep, DepT) :-
  pivot_a( X, Indep, t_U(U), DepT),		% changes the basis
  rcbl( Continuation, B0,B1).

narrow_u( t_u(_),    X, U) :- put_atts( X, type(t_u(U))).
narrow_u( t_lu(L,_), X, U) :- put_atts( X, type(t_lu(L,U))).

narrow_l( t_l(_),    X, L) :- put_atts( X, type(t_l(L))).
narrow_l( t_lu(_,U), X, L) :- put_atts( X, type(t_lu(L,U))).

% ----------------------------------- dump -------------------------------------

dump_var( t_none, V, I,H) --> !,
  ( { H=[W*K],V==W,arith_eval(I=:=0),arith_eval(K=:=1) } ->	% indep var
      []
  ;
      { nf2sum( H, I, Sum) },
      [ V = Sum ]
  ).
dump_var( t_L(L), V, I,H) --> !, dump_var( t_l(L),    V, I,H).
dump_var( t_l(L), V, I,H) --> !,
  {
    H= [_*K|_], 			% avoid 1 >= 0
    get_atts( V, strictness(Strict)),
    Sm is Strict /\ 2'10,
    arith_eval( 1/K, Kr),
    arith_eval( Kr*(L-I), Li),
    mult_hom( H, Kr, H1),
    arith_eval( 0, Z), nf2sum( H1, Z, Sum),
    ( arith_eval( K > 0) ->
	dump_strict( Sm, Sum >= Li, Sum > Li, Result)
    ;
	dump_strict( Sm, Sum =< Li, Sum < Li, Result)
    )
  },
  [ Result ].
dump_var( t_U(U), V, I,H) --> !, dump_var( t_u(U),    V, I,H).
dump_var( t_u(U), V, I,H) --> !,
  {
    H= [_*K|_], 			% avoid 0 =< 1
    get_atts( V, strictness(Strict)),
    Sm is Strict /\ 2'01,
    arith_eval( 1/K, Kr),
    arith_eval( Kr*(U-I), Ui),
    mult_hom( H, Kr, H1),
    arith_eval( 0, Z), nf2sum( H1, Z, Sum),
    ( arith_eval( K > 0) ->
	dump_strict( Sm, Sum =< Ui, Sum < Ui, Result)
    ;
	dump_strict( Sm, Sum >= Ui, Sum > Ui, Result)
    )
  },
  [ Result ].
dump_var( t_Lu(L,U), V, I,H) --> !, dump_var( t_l(L), V,I,H),
				    dump_var( t_U(U), V,I,H).
dump_var( t_lU(L,U), V, I,H) --> !, dump_var( t_l(L), V,I,H),
				    dump_var( t_U(U), V,I,H).
dump_var( t_lu(L,U), V, I,H) --> !, dump_var( t_l(L), V,I,H),
				    dump_var( t_U(U), V,I,H).
dump_var( T,	  V, I,H) -->
  [ V:T:I+H ].

dump_strict( 0, Result, _, Result).
dump_strict( 1, _, Result, Result).
dump_strict( 2, _, Result, Result).

dump_nz( _, H, I) -->
  {
    H = [_*K|_],
    arith_eval( 1/K, Kr),
    arith_eval( -Kr*I, I1),
    mult_hom( H, Kr, H1),
    arith_eval( 0, Z), nf2sum( H1, Z, Sum)
  },
  [ Sum =\= I1 ].