This repository has been archived on 2023-08-20. You can view files and clone it, but cannot push or open issues or pull requests.
vsc e5f4633c39 This commit was generated by cvs2svn to compensate for changes in r4,
which included commits to RCS files with non-trunk default branches.


git-svn-id: https://yap.svn.sf.net/svnroot/yap/trunk@5 b08c6af1-5177-4d33-ba66-4b1c6b8b522a
2001-04-09 19:54:03 +00:00

1257 lines
32 KiB
Prolog

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.
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) :-
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),
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),
( 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)),
( 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)),
( 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
solve( H, Lin, _, Bind0,BindT) :-
H = [_|_], % indexing
%
% [] is an empty ord_set, anything will be preferred
% over 9-9
%
sd( H, [],ClassesUniq, 9-9-0,Category-Selected-_, NV,NVT),
isolate( Selected, Lin, Lin1),
( 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),
( 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),
( 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),
( 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 ].