1/*
    2
    3    Part of CLP(Q) (Constraint Logic Programming over Rationals)
    4
    5    Author:        Leslie De Koninck
    6    E-mail:        Leslie.DeKoninck@cs.kuleuven.be
    7    WWW:           http://www.swi-prolog.org
    8		   http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09
    9    Copyright (C): 2006, K.U. Leuven and
   10		   1992-1995, Austrian Research Institute for
   11		              Artificial Intelligence (OFAI),
   12			      Vienna, Austria
   13
   14    This software is based on CLP(Q,R) by Christian Holzbaur for SICStus
   15    Prolog and distributed under the license details below with permission from
   16    all mentioned authors.
   17
   18    This program is free software; you can redistribute it and/or
   19    modify it under the terms of the GNU General Public License
   20    as published by the Free Software Foundation; either version 2
   21    of the License, or (at your option) any later version.
   22
   23    This program is distributed in the hope that it will be useful,
   24    but WITHOUT ANY WARRANTY; without even the implied warranty of
   25    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   26    GNU General Public License for more details.
   27
   28    You should have received a copy of the GNU Lesser General Public
   29    License along with this library; if not, write to the Free Software
   30    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
   31
   32    As a special exception, if you link this library with other files,
   33    compiled with a Free Software compiler, to produce an executable, this
   34    library does not by itself cause the resulting executable to be covered
   35    by the GNU General Public License. This exception does not however
   36    invalidate any other reasons why the executable file might be covered by
   37    the GNU General Public License.
   38*/
   39
   40:- module(clpcd_bv,
   41          [
   42              allvars/2,
   43              detach_bounds/2,
   44              detach_bounds_vlv/6,
   45              iterate_dec/3,
   46              var_with_def_assign/3,
   47              vertex_value/3,
   48              maximize/2,
   49              minimize/2,
   50              sup/3,
   51              sup/5,
   52              inf/3,
   53              inf/5
   54          ]).   55
   56:- use_module(library(clpcd/class)).   57:- use_module(library(clpcd/domain_ops)).   58:- use_module(library(clpcd/indep)).   59:- use_module(library(clpcd/nf)).   60:- use_module(library(clpcd/solve)).   61:- use_module(library(clpcd/store)).   62:- use_module(library(clpcd/detact)).   63
   64% For the rhs maint. the following events are important:
   65%
   66%	-) introduction of an indep var at active bound B
   67%	-) narrowing of active bound
   68%	-) swap active bound
   69%	-) pivot
   70%
   71
   72% a variables bound (L/U) can have the states:
   73%
   74%	-) t_none	no bounds
   75%	-) t_l		inactive lower bound
   76%	-) t_u		inactive upper bound
   77%	-) t_L		active lower bound
   78%	-) t_U		active upper bound
   79%	-) t_lu		inactive lower and upper bound
   80%	-) t_Lu		active lower bound and inactive upper bound
   81%	-) t_lU		inactive lower bound and active upper bound
   82
   83% TODO
   84%
   85%
   86
   87var_with_def_assign(CLP,Var,Lin) :-
   88	Lin = [I,_|Hom],
   89	(   Hom = []
   90	->  % X=k
   91	    Var = I
   92	;   Hom = [l(V*K,_)|Cs]
   93	->  (   Cs = [],
   94		compare_d(CLP, =, K, 1),	% K =:= 1
   95		compare_d(CLP, =, 0, I)
   96	    ->	% X=Y
   97		Var = V
   98	    ;	% general case
   99		var_with_def_intern(t_none,CLP,Var,Lin,0)
  100	    )
  101	).
  102
  103maximize(CLP, Term) :-
  104	minimize(CLP, -Term).
  105
  106%
  107% This is NOT coded as minimize(Expr) :- inf(Expr,Expr).
  108%
  109% because the new version of inf/2 only visits
  110% the vertex where the infimum is assumed and returns
  111% to the 'current' vertex via backtracking.
  112% The rationale behind this construction is to eliminate
  113% all garbage in the solver data structures produced by
  114% the pivots on the way to the extremal point caused by
  115% {inf,sup}/{2,4}.
  116%
  117% If we are after the infimum/supremum for minimizing/maximizing,
  118% this strategy may have adverse effects on performance because
  119% the simplex algorithm is forced to re-discover the
  120% extremal vertex through the equation {Inf =:= Expr}.
  121%
  122% Thus the extra code for {minimize,maximize}/1.
  123%
  124% In case someone comes up with an example where
  125%
  126%   inf(Expr,Expr)
  127%
  128% outperforms the provided formulation for minimize - so be it.
  129% Both forms are available to the user.
  130%
  131minimize(CLP,Term) :-
  132	wait_linear(CLP, Term, Nf, minimize_lin(CLP,Nf)).
  133
  134% minimize_lin(CLP,Lin)
  135%
  136% Minimizes the linear expression Lin. It does so by making a new
  137% variable Dep and minimizes its value.
  138
  139minimize_lin(CLP,Lin) :-
  140	deref(CLP,Lin,Lind),
  141	var_with_def_intern(t_none,CLP,Dep,Lind,0),
  142	determine_active_dec(CLP, Lind),
  143	iterate_dec(CLP, Dep, Inf),
  144	add_constraint(Dep =:= Inf, CLP).
  145
  146sup(CLP,Expression,Sup) :-
  147	sup(CLP,Expression,Sup,[],[]).
  148
  149sup(CLP,Expression,Sup,Vector,Vertex) :-
  150	inf(CLP,-Expression,-Sup,Vector,Vertex).
  151
  152inf(CLP,Expression,Inf) :-
  153	inf(CLP,Expression,Inf,[],[]).
  154
  155inf(CLP,Expression,Inf,Vector,Vertex) :-
  156	% wait until Expression becomes linear, Nf contains linear Expression
  157	% in normal form
  158	wait_linear(CLP, Expression, Nf, inf_lin(CLP,Nf,Inf,Vector,Vertex)).
  159
  160inf_lin(CLP, Lin,_,Vector,_) :-
  161	deref(CLP,Lin,Lind),
  162	var_with_def_intern(t_none,CLP,Dep,Lind,0),	% make new variable Dep = Lind
  163	determine_active_dec(CLP, Lind),	% minimizes Lind
  164	iterate_dec(CLP, Dep, Inf),
  165	vertex_value(Vector, CLP, Values),
  166	nb_setval(inf,[Inf|Values]),
  167	fail.
  168inf_lin(CLP, _,Infimum,_,Vertex) :-
  169	nb_current(inf,L),
  170	nb_delete(inf),
  171	assign([Infimum|Vertex],CLP,L).
  172
  173% vertex_value(Vars,Values)
  174%
  175% Returns in <Values> the current values of the variables in <Vars>.
  176
  177vertex_value([], _, []).
  178vertex_value([X|Xs], CLP, [V|Vs]) :-
  179	rhs_value(CLP, X, V),
  180	vertex_value(Xs, CLP, Vs).
  181
  182% rhs_value(X,Value)
  183%
  184% Returns in <Value> the current value of variable <X>.
  185
  186rhs_value(CLP, Xn, Value) :-
  187	(   nonvar(Xn)
  188	->  Value = Xn
  189	;   var(Xn)
  190	->  deref_var(CLP, Xn, Xd),
  191	    Xd = [I,R|_],
  192	    eval_d(CLP, R+I, Value)
  193	).
  194
  195% assign(L1,CLP,L2)
  196%
  197% The elements of L1 are pairwise assigned to the elements of L2
  198% by means of asserting {X =:= Y} where X is an element of L1 and Y
  199% is the corresponding element of L2.
  200
  201assign([],_,[]).
  202assign([X|Xs],CLP,[Y|Ys]) :-
  203	add_constraint(X =:= Y, CLP), % more defensive/expressive than X=Y
  204	assign(Xs,CLP,Ys).
  205
  206% --------------------------------- optimization ------------------------------
  207%
  208% The _sn(S) =< 0 row might be temporarily infeasible.
  209% We use reconsider/2 to fix this.
  210%
  211%   s(S) e [_,0] = d +xi ... -xj, Rhs > 0 so we want to decrease s(S)
  212%
  213%   positive xi would have to be moved towards their lower bound,
  214%   negative xj would have to be moved towards their upper bound,
  215%
  216%   the row s(S) does not limit the lower bound of xi
  217%   the row s(S) does not limit the upper bound of xj
  218%
  219%   a) if some other row R is limiting xk, we pivot(R,xk),
  220%      s(S) will decrease and get more feasible until (b)
  221%   b) if there is no limiting row for some xi: we pivot(s(S),xi)
  222%					    xj: we pivot(s(S),xj)
  223%      which cures the infeasibility in one step
  224%
  225
  226
  227% iterate_dec(OptVar,Opt)
  228%
  229% Decreases the bound on the variables of the linear equation of OptVar as much
  230% as possible and returns the resulting optimal bound in Opt. Fails if for some
  231% variable, a status of unlimited is found.
  232
  233iterate_dec(CLP, OptVar, Opt) :-
  234	get_attr(OptVar,clpcd_itf,Att),
  235	arg(4,Att,lin([I,R|H])),
  236	dec_step(H, CLP, Status),
  237	(   Status = applied
  238	->  iterate_dec(CLP, OptVar, Opt)
  239	;   Status = optimum,
  240	    eval_d(CLP, R + I, Opt)
  241	).
  242
  243% allvars(X,Allvars)
  244%
  245% Allvars is a list of all variables in the class to which X belongs.
  246
  247allvars(X,Allvars) :-
  248	get_attr(X,clpcd_itf,Att),
  249	arg(6,Att,class(C)),
  250	class_allvars(C,Allvars).
  251
  252%
  253% Careful when an indep turns into t_none !!!
  254%
  255
  256detach_bounds(CLP, V) :-
  257	get_attr(V,clpcd_itf,Att),
  258	arg(2,Att,type(Type)),
  259	arg(4,Att,lin(Lin)),
  260	arg(5,Att,order(OrdV)),
  261	arg(6,Att,class(Class)),
  262	setarg(2,Att,type(t_none)),
  263	setarg(3,Att,strictness(0)),
  264	(   indep(CLP, Lin, OrdV)
  265	->  (   ub(Class,CLP,OrdV,Vub-Vb-_)
  266	    ->	% exchange against thightest
  267		class_basis_drop(Class,Vub),
  268		pivot(CLP, Vub, Class, OrdV, Vb, Type)
  269	    ;   lb(Class,CLP,OrdV,Vlb-Vb-_)
  270	    ->  class_basis_drop(Class,Vlb),
  271		pivot(CLP, Vlb, Class, OrdV, Vb, Type)
  272	    ;   true
  273	    )
  274	;   class_basis_drop(Class,V)
  275	).
  276
  277detach_bounds_vlv(CLP,OrdV,Lin,Class,Var,NewLin) :-
  278	(   indep(CLP, Lin, OrdV)
  279	->  Lin = [_,R|_],
  280	    (   ub(Class,CLP,OrdV,Vub-Vb-_)
  281	    ->  % in verify_lin, class might contain two occurrences of Var,
  282		% but it doesn't matter which one we delete
  283		class_basis_drop(Class,Var),
  284		pivot_vlv(CLP, Vub, Class, OrdV, Vb, R, NewLin)
  285	    ;   lb(Class,CLP,OrdV,Vlb-Vb-_)
  286	    ->  class_basis_drop(Class,Var),
  287		pivot_vlv(CLP, Vlb, Class, OrdV, Vb, R, NewLin)
  288	    ;   NewLin = Lin
  289	    )
  290	;   NewLin = Lin,
  291	    class_basis_drop(Class,Var)
  292	).
  293
  294% Rewrite Dep = ... + Coeff*Indep + ...
  295% into Indep = ... + -1/Coeff*Dep + ...
  296%
  297% For backsubstitution, old current value of Indep must be removed from RHS
  298% New current value of Dep must be added to RHS
  299% For solving: old current value of Indep should be out of RHS
  300
  301pivot_vlv(CLP,Dep,Class,IndepOrd,DepAct,AbvI,Lin) :-
  302	get_attr(Dep,clpcd_itf,Att),
  303	arg(4,Att,lin(H)),
  304	arg(5,Att,order(DepOrd)),
  305	setarg(2,Att,type(DepAct)),
  306	select_active_bound(DepAct,AbvD), % New current value for Dep
  307	delete_factor(IndepOrd,H,H0,Coeff), % Dep = ... + Coeff*Indep + ...
  308	eval_d(CLP, -AbvD, AbvDm),
  309	eval_d(CLP, -AbvI, AbvIm),
  310	add_linear_f1(CLP, [0,AbvIm], Coeff, H0, H1),
  311	div_d(CLP, -1, Coeff, K),
  312	add_linear_ff(CLP, H1, K, [0,AbvDm,l(Dep* -1,DepOrd)], K, Lin),
  313	    % Indep = -1/Coeff*... + 1/Coeff*Dep
  314	add_linear_11(CLP, Lin, [0,AbvIm], SubstLin),
  315	backsubst(CLP, Class, IndepOrd, SubstLin)