diff --git a/LGPL/clpr/Makefile.in b/LGPL/clpr/Makefile.in new file mode 100644 index 000000000..2c059c801 --- /dev/null +++ b/LGPL/clpr/Makefile.in @@ -0,0 +1,71 @@ +# +# default base directory for YAP installation +# (EROOT for architecture-dependent files) +# +prefix = @prefix@ +ROOTDIR = $(prefix) +EROOTDIR = @exec_prefix@ + +SHELL=@SHELL@ +PLBASE=@PLBASE@ +PLARCH=@PLARCH@ +PL=@PL@ +XPCEBASE=$(PLBASE)/xpce +PKGDOC=$(PLBASE)/doc/packages +PCEHOME=../../xpce +LIBDIR=$(PLBASE)/library +SHAREDIR=$(ROOTDIR)/share/Yap +CLPRDIR=$(SHAREDIR)/clpr +EXDIR=$(PKGDOC)/examples/clpr +DESTDIR= +srcdir=@srcdir@ +CLPRSOURCEDIR=$(srcdir)/clpr + +INSTALL=@INSTALL@ +INSTALL_PROGRAM=@INSTALL_PROGRAM@ +INSTALL_DATA=@INSTALL_DATA@ + +CLPRPRIV= $(CLPRSOURCEDIR)/arith_r.pl $(CLPRSOURCEDIR)/bb.pl $(CLPRSOURCEDIR)/bv.pl $(CLPRSOURCEDIR)/class.pl $(CLPRSOURCEDIR)/dump.pl $(CLPRSOURCEDIR)/fourmotz.pl \ + $(CLPRSOURCEDIR)/geler.pl $(CLPRSOURCEDIR)/ineq.pl $(CLPRSOURCEDIR)/itf3.pl $(CLPRSOURCEDIR)/nf.pl $(CLPRSOURCEDIR)/nfr.pl $(CLPRSOURCEDIR)/ordering.pl \ + $(CLPRSOURCEDIR)/project.pl $(CLPRSOURCEDIR)/redund.pl $(CLPRSOURCEDIR)/store.pl $(CLPRSOURCEDIR)/ugraphs.pl +LIBPL= $(srcdir)/clpr.yap $(srcdir)/clpr.pl +EXAMPLES= + +all:: + @echo "Nothing to be done for this package" + +install: $(LIBPL) + mkdir -p $(DESTDIR)$(CLPRDIR) + $(INSTALL_DATA) $(LIBPL) $(DESTDIR)$(SHAREDIR) + $(INSTALL_DATA) $(CLPRPRIV) $(DESTDIR)$(CLPRDIR) + $(INSTALL_DATA) $(srcdir)/README $(DESTDIR)$(CLPRDIR) + +rpm-install: install + +pdf-install: install-examples + +html-install: install-examples + +install-examples:: +# mkdir -p $(DESTDIR)$(EXDIR) +# (cd Examples && $(INSTALL_DATA) $(EXAMPLES) $(DESTDIR)$(EXDIR)) + +uninstall: + (cd $(CLPDIR) && rm -f $(LIBPL)) + rm -rf $(CLPRDIR) + +check:: +# $(PL) -q -f clpr_test.pl -g test,halt -t 'halt(1)' + + +################################################################ +# Clean +################################################################ + +clean: + rm -f *~ *% config.log + +distclean: clean + rm -f config.h config.cache config.status Makefile + rm -rf autom4te.cache + diff --git a/LGPL/clpr/README b/LGPL/clpr/README new file mode 100644 index 000000000..0e9046fb6 --- /dev/null +++ b/LGPL/clpr/README @@ -0,0 +1,22 @@ + SWI-Prolog CLP(R) + ----------------- + + +Author: Leslie De Koninck, K.U. Leuven as part of a thesis with supervisor + Bart Demoen and daily advisor Tom Schrijvers. Integrated into + the main SWI-Prolog source by Jan Wielemaker. + +This software is based on the CLP(Q,R) implementation by Christian +Holzbauer and released with permission from all above mentioned authors +and Christian Holzbauer under the standard SWI-Prolog license schema: +GPL-2 + statement to allow linking with proprietary software. + +The sources of this package are maintained in packages/clpr in the +SWI-Prolog source distribution. The documentation source is in +man/lib/clpr.doc as part of the overall SWI-Prolog documentation. + +Full documentation on CLP(Q,R) can be found at + + http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 + + diff --git a/LGPL/clpr/clpr.pl b/LGPL/clpr/clpr.pl new file mode 100644 index 000000000..742126574 --- /dev/null +++ b/LGPL/clpr/clpr.pl @@ -0,0 +1,10 @@ + +:- load_files(library(swi),[silent(true)]). + +:- yap_flag(unknown,error). + +:- include('clpr.pl'). + + + + diff --git a/LGPL/clpr/clpr.yap b/LGPL/clpr/clpr.yap new file mode 100644 index 000000000..742126574 --- /dev/null +++ b/LGPL/clpr/clpr.yap @@ -0,0 +1,10 @@ + +:- load_files(library(swi),[silent(true)]). + +:- yap_flag(unknown,error). + +:- include('clpr.pl'). + + + + diff --git a/LGPL/clpr/clpr/arith_r.pl b/LGPL/clpr/clpr/arith_r.pl new file mode 100644 index 000000000..ce4a33ab8 --- /dev/null +++ b/LGPL/clpr/clpr/arith_r.pl @@ -0,0 +1,85 @@ +/* $Id: arith_r.pl,v 1.1 2005-10-28 17:51:01 vsc Exp $ + + Part of CPL(R) (Constraint Logic Programming over Reals) + + Author: Leslie De Koninck + E-mail: Tom.Schrijvers@cs.kuleuven.ac.be + WWW: http://www.swi-prolog.org + http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 + Copyright (C): 2004, K.U. Leuven and + 1992-1995, Austrian Research Institute for + Artificial Intelligence (OFAI), + Vienna, Austria + + This software is part of Leslie De Koninck's master thesis, supervised + by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) + by Christian Holzbaur for SICStus Prolog and distributed under the + license details below with permission from all mentioned authors. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 + of the License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + As a special exception, if you link this library with other files, + compiled with a Free Software compiler, to produce an executable, this + library does not by itself cause the resulting executable to be covered + by the GNU General Public License. This exception does not however + invalidate any other reasons why the executable file might be covered by + the GNU General Public License. +*/ + +:- module(arith_r, + [ + arith_eps/1, + arith_normalize/2, + integerp/1, + integerp/2 + ]). + +arith_eps(1.0e-10). % for Monash #zero expansion 1.0e-12 +eps(1.0e-10,-1.0e-10). + +arith_normalize(X,Norm) :- + var(X), + !, + raise_exception(instantiation_error(arith_normalize(X,Norm),1)). +arith_normalize(rat(N,D),Norm) :- rat_float(N,D,Norm). +arith_normalize(X,Norm) :- + number(X), + Norm is float(X). + +integerp(X) :- + floor(/*float?*/X)=:=X. + +integerp(X,I) :- + floor(/*float?*/X)=:=X, + I is integer(X). + +% copied from arith.pl. + +rat_float(Nx,Dx,F) :- + limit_encoding_length(Nx,Dx,1023,Nxl,Dxl), + F is Nxl / Dxl. + +limit_encoding_length(0,D,_,0,D) :- !. % msb ... +limit_encoding_length(N,D,Bits,Nl,Dl) :- + Shift is min(max(msb(abs(N)),msb(D))-Bits, + min(msb(abs(N)),msb(D))), + Shift > 0, + !, + Ns is N>>Shift, + Ds is D>>Shift, + Gcd is gcd(Ns,Ds), + Nl is Ns//Gcd, + Dl is Ds//Gcd. +limit_encoding_length(N,D,_,N,D). diff --git a/LGPL/clpr/clpr/bb.pl b/LGPL/clpr/clpr/bb.pl new file mode 100644 index 000000000..4446c5f9c --- /dev/null +++ b/LGPL/clpr/clpr/bb.pl @@ -0,0 +1,301 @@ +/* $Id: bb.pl,v 1.1 2005-10-28 17:51:01 vsc Exp $ + + Part of CPL(R) (Constraint Logic Programming over Reals) + + Author: Leslie De Koninck + E-mail: Tom.Schrijvers@cs.kuleuven.ac.be + WWW: http://www.swi-prolog.org + http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 + Copyright (C): 2004, K.U. Leuven and + 1992-1995, Austrian Research Institute for + Artificial Intelligence (OFAI), + Vienna, Austria + + This software is part of Leslie De Koninck's master thesis, supervised + by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) + by Christian Holzbaur for SICStus Prolog and distributed under the + license details below with permission from all mentioned authors. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 + of the License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + As a special exception, if you link this library with other files, + compiled with a Free Software compiler, to produce an executable, this + library does not by itself cause the resulting executable to be covered + by the GNU General Public License. This exception does not however + invalidate any other reasons why the executable file might be covered by + the GNU General Public License. +*/ + + +:- module(bb, + [ + bb_inf/3, + bb_inf/5, + vertex_value/2 + ]). +:- use_module(bv, + [ + deref/2, + deref_var/2, + determine_active_dec/1, + inf/2, + iterate_dec/2, + sup/2, + var_with_def_assign/2 + ]). +:- use_module(nf, + [ {}/1, + entailed/1, + nf/2, + nf_constant/2, + repair/2, + wait_linear/3 + ]). + + + +:- dynamic blackboard_incumbent/1. + +% bb_inf(Ints,Term,Inf) +% +% Finds the infimum of Term where the variables Ints are to be integers. +% The infimum is stored in Inf. + +bb_inf(Is,Term,Inf) :- + bb_inf(Is,Term,Inf,_,0.001). + +bb_inf(Is,Term,Inf,Vertex,Eps) :- + nf(Eps,ENf), + nf_constant(ENf,EpsN), + wait_linear(Term,Nf,bb_inf_internal(Is,Nf,EpsN,Inf,Vertex)). + +% --------------------------------------------------------------------- + +% bb_inf_internal(Is,Lin,Eps,Inf,Vertex) +% +% Finds an infimum Inf for linear expression in normal form Lin, where +% all variables in Is are to be integers. Eps denotes the margin in which +% we accept a number as an integer (to deal with rounding errors etc.). + +bb_inf_internal(Is,Lin,Eps,_,_) :- + bb_intern(Is,IsNf,Eps), + ( + retract(blackboard_incumbent(_)) -> + + true + ; + true + ), + repair(Lin,LinR), % bb_narrow ... + deref(LinR,Lind), + var_with_def_assign(Dep,Lind), + determine_active_dec(Lind), + bb_loop(Dep,IsNf,Eps), + fail. +bb_inf_internal(_,_,_,Inf,Vertex) :- + retract(blackboard_incumbent(InfVal-Vertex)), % GC + { Inf =:= InfVal }. + +% bb_loop(Opt,Is,Eps) +% +% Minimizes the value of Opt where variables Is have to be integer values. +% Eps denotes the rounding error that is acceptable. This predicate can be +% backtracked to try different strategies. + +bb_loop(Opt,Is,Eps) :- + bb_reoptimize(Opt,Inf), + bb_better_bound(Inf), + vertex_value(Is,Ivs), + ( + bb_first_nonint(Is,Ivs,Eps,Viol,Floor,Ceiling) -> + + bb_branch(Viol,Floor,Ceiling), + bb_loop(Opt,Is,Eps) + ; + round_values(Ivs,RoundVertex), + ( + retract(blackboard_incumbent(_)) -> + + assert(blackboard_incumbent(Inf-RoundVertex)) + ; + assert(blackboard_incumbent(Inf-RoundVertex)) + ) + ). + +% bb_reoptimize(Obj,Inf) +% +% Minimizes the value of Obj and puts the result in Inf. +% This new minimization is necessary as making a bound integer may yield a +% different optimum. The added inequalities may also have led to binding. + +bb_reoptimize(Obj,Inf) :- + var(Obj), + iterate_dec(Obj,Inf). +bb_reoptimize(Obj,Inf) :- + nonvar(Obj), + Inf = Obj. + +% bb_better_bound(Inf) +% +% Checks if the new infimum Inf is better than the previous one (if such exists). + +bb_better_bound(Inf) :- + blackboard_incumbent(Inc-_), + !, + Inf - Inc < -1e-010. +bb_better_bound(_). + +% bb_branch(V,U,L) +% +% Stores that V =< U or V >= L, can be used for different strategies within bb_loop/3. + +bb_branch(V,U,_) :- {V =< U}. +bb_branch(V,_,L) :- {V >= L}. + +% vertex_value(Vars,Rhss) +% +% Returns in Rhss, the Rhs values corresponding to the Vars via rhs_value/2. + +vertex_value([],[]). +vertex_value([X|Xs],[V|Vs]) :- + rhs_value(X,V), + vertex_value(Xs,Vs). + +% rhs_value(X,Rhs) +% +% Returns the Rhs value of variable X. This is the value by which the +% variable is actively bounded. If X is a nonvar, the value of X is returned. + +rhs_value(Xn,Value) :- + ( + nonvar(Xn) -> + + Value = Xn + ; + var(Xn) -> + + deref_var(Xn,Xd), + Xd = [I,R|_], + Value is R+I + ). + +% bb_first_nonint(Ints,Rhss,Eps,Viol,Floor,Ceiling) +% +% Finds the first variable in Ints which doesn't have an active integer bound. +% Rhss contain the Rhs (R + I) values corresponding to the variables. +% The first variable that hasn't got an active integer bound, is returned in +% Viol. The floor and ceiling of its actual bound is returned in Floor and Ceiling. + +bb_first_nonint([I|Is],[Rhs|Rhss],Eps,Viol,F,C) :- + ( + Floor is float(floor(Rhs)), + Ceiling is float(ceiling(Rhs)), + Eps - min(Rhs-Floor,Ceiling-Rhs) < -1e-010 -> + + Viol = I, + F = Floor, + C = Ceiling + ; + bb_first_nonint(Is,Rhss,Eps,Viol,F,C) + ). + +% round_values([X|Xs],[Xr|Xrs]) +% +% Rounds of the values of the first list into the second list. + +round_values([],[]). +round_values([X|Xs],[Y|Ys]) :- + Y = float(round(X)), + round_values(Xs,Ys). + +% bb_intern([X|Xs],[Xi|Xis],Eps) +% +% Turns the elements of the first list into integers into the second +% list via bb_intern/4. + +bb_intern([],[],_). +bb_intern([X|Xs],[Xi|Xis],Eps) :- + nf(X,Xnf), + bb_intern(Xnf,Xi,X,Eps), + bb_intern(Xs,Xis,Eps). + + +% bb_intern(Nf,X,Term,Eps) +% +% Makes sure that Term which is normalized into Nf, is integer. +% X contains the possibly changed Term. If Term is a variable, +% then its bounds are hightened or lowered to the next integer. +% Otherwise, it is checked it Term is integer. + +bb_intern([],X,_,_) :- + !, + X = 0.0. +bb_intern([v(I,[])],X,_,Eps) :- + !, + X = I, + min(I-float(floor(I)),float(ceiling(I))-I) - Eps < 1e-010. +bb_intern([v(One,[V^1])],X,_,_) :- + Test is One - 1.0, + Test =< 1e-010, + Test >= -1e-010, + !, + V = X, + bb_narrow_lower(X), + bb_narrow_upper(X). +bb_intern(_,_,Term,_) :- + raise_exception(instantiation_error(bb_inf(Term,_,_),1)). + +% bb_narrow_lower(X) +% +% Narrows the lower bound so that it is an integer bound. +% We do this by finding the infimum of X and asserting that X +% is larger than the first integer larger or equal to the infimum +% (second integer if X is to be strict larger than the first integer). + +bb_narrow_lower(X) :- + ( + inf(X,Inf) -> + + Bound is float(ceiling(Inf)), + ( + entailed(X > Bound) -> + + {X >= Bound+1} + ; + {X >= Bound} + ) + ; + true + ). + +% bb_narrow_upper(X) +% +% See bb_narrow_lower/1. This predicate handles the upper bound. + +bb_narrow_upper(X) :- + ( + sup(X,Sup) -> + + Bound is float(floor(Sup)), + ( + entailed(X < Bound) -> + + {X =< Bound-1} + ; + {X =< Bound} + ) + ; + true + ). diff --git a/LGPL/clpr/clpr/bv.pl b/LGPL/clpr/clpr/bv.pl new file mode 100644 index 000000000..7ec07e986 --- /dev/null +++ b/LGPL/clpr/clpr/bv.pl @@ -0,0 +1,1805 @@ +/* $Id: bv.pl,v 1.1 2005-10-28 17:51:01 vsc Exp $ + + Part of CPL(R) (Constraint Logic Programming over Reals) + + Author: Leslie De Koninck + E-mail: Tom.Schrijvers@cs.kuleuven.ac.be + WWW: http://www.swi-prolog.org + http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 + Copyright (C): 2004, K.U. Leuven and + 1992-1995, Austrian Research Institute for + Artificial Intelligence (OFAI), + Vienna, Austria + + This software is part of Leslie De Koninck's master thesis, supervised + by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) + by Christian Holzbaur for SICStus Prolog and distributed under the + license details below with permission from all mentioned authors. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 + of the License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + As a special exception, if you link this library with other files, + compiled with a Free Software compiler, to produce an executable, this + library does not by itself cause the resulting executable to be covered + by the GNU General Public License. This exception does not however + invalidate any other reasons why the executable file might be covered by + the GNU General Public License. +*/ + + +:- module(bv, + [ + allvars/2, + backsubst/3, + backsubst_delta/4, + basis_add/2, + dec_step/2, + deref/2, + deref_var/2, + detach_bounds/1, + detach_bounds_vlv/5, + determine_active_dec/1, + determine_active_inc/1, + dump_var/6, + dump_nz/5, + export_binding/1, + export_binding/2, + get_or_add_class/2, + inc_step/2, + intro_at/3, + iterate_dec/2, + lb/3, + pivot_a/4, + pivot/4, + rcbl_status/6, + reconsider/1, + same_class/2, + solve/1, + solve_ord_x/3, + ub/3, + unconstrained/4, + var_intern/2, + var_intern/3, + var_with_def_assign/2, + var_with_def_intern/4, + maximize/1, + minimize/1, + sup/2, + sup/4, + inf/2, + inf/4, + 'solve_<'/1, + 'solve_=<'/1, + 'solve_=\\='/1, + log_deref/4 + ]). +:- use_module(store, + [ + add_linear_11/3, + add_linear_ff/5, + delete_factor/4, + indep/2, + isolate/3, + nf2sum/3, + nf_rhs_x/4, + nf_substitute/4, + normalize_scalar/2, + mult_hom/3, + mult_linear_factor/3 + ]). +:- use_module(class, + [ + class_allvars/2, + class_basis/2, + class_basis_add/3, + class_basis_drop/2, + class_basis_pivot/3, + class_new/4 + ]). +:- use_module(ineq, + [ + ineq/4 + ]). +:- use_module(nf, + [ {}/1, + split/3, + wait_linear/3 + ]). +:- use_module(bb, + [ + vertex_value/2 + ]). +:- use_module(library(ordsets), + [ + ord_add_element/3 + ]). + +:- dynamic blackboard_inf/1. + +% 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 ------------------------------------ % + +% deref(Lin,Lind) +% +% Makes a linear equation of the form [v(I,[])|H] into a solvable linear equation. +% If the variables are new, they are initialized with the linear equation X=X. + +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). + +% log_deref(Len,[Vs|VsTail],VsTail,Res) +% +% Logarithmically converts a linear equation in normal form ([v(_,_)|_]) into a +% linear equation in solver form ([I,R,K*X|_]). Res contains the result, Len is +% the length of the part to convert and [Vs|VsTail] is a difference list containing +% the equation in normal form. + +log_deref(0,Vs,Vs,Lin) :- + !, + Lin = [0.0,0.0]. +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). + +% deref_var(X,Lin) +% +% Returns the equation of variable X. If X is a new variable, a new equation X = X is made. + +deref_var(X,Lin) :- + get_attr(X,itf3,(_,_,lin(Lin),_)), + !. +deref_var(X,Lin) :- % create a linear var + Lin = [0.0,0.0,l(X*1.0,Ord)], + put_attr(X,itf3,(type(t_none),strictness(0),lin(Lin),order(Ord),n,n,n,n,n,n)). + +% TODO +% +% + +var_with_def_assign(Var,Lin) :- + Lin = [I,_|Hom], + ( + Hom = [] -> % X=k + + Var = I + ; + Hom = [l(V*K,_)|Cs] -> + + ( + Cs = [], + TestK is K - 1.0, % K =:= 1 + TestK =< 1e-010, + TestK >= -1e-010, + I >= -1e-010, % I =:= 0 + I =< 1e-010 -> % X=Y + + Var = V + ; % general case + var_with_def_intern(t_none,Var,Lin,0) + ) + ). + +% var_with_def_intern(Type,Var,Lin,Strictness) +% +% Makes Lin the linear equation of new variable Var, makes all variables of Lin, and Var +% of the same class and bounds Var by type(Type) and strictness(Strictness) + +var_with_def_intern(Type,Var,Lin,Strict) :- + put_attr(Var,itf3,(type(Type),strictness(Strict),lin(Lin),order(_),n,n,n,n,n,n)), % check uses + Lin = [_,_|Hom], + get_or_add_class(Var,Class), + same_class(Hom,Class). + +% TODO +% +% + +var_intern(Type,Var,Strict) :- + Lin = [0.0,0.0,l(Var*1.0,Ord)], + put_attr(Var,itf3,(type(Type),strictness(Strict),lin(Lin),order(Ord),n,n,n,n,n,n)), + get_or_add_class(Var,_Class). + +% TODO +% +% + +var_intern(Var,Class) :- % for ordered/1 but otherwise free vars + get_attr(Var,itf3,(type(_),_,lin(_),_)), + !, + get_or_add_class(Var,Class). +var_intern(Var,Class) :- + put_attr(Var,itf3,(type(t_none),strictness(0),lin([0.0,0.0,l(Var*1.0,Ord)]),order(Ord),n,n,n,n,n,n)), + get_or_add_class(Var,Class). + +% ------------------------------------------------------------------------------ + +% export_binding(Lst) +% +% Binds variables X to Y where Lst contains elements of the form [X-Y]. + +export_binding([]). +export_binding([X-Y|Gs]) :- + export_binding(Y,X), + export_binding(Gs). + +% export_binding(Y,X) +% +% Binds variable X to Y. If Y is a nonvar and equals 0, then X is set to 0 +% (numerically more stable) + +export_binding(Y,X) :- + var(Y), + Y = X. +export_binding(Y,X) :- + nonvar(Y), + ( + Y >= -1e-010, % Y =:= 0 + Y =< 1e-010 -> + + X = 0.0 + ; + Y = X + ). + +% 'solve_='(Nf) +% +% Solves linear equation Nf = 0 where Nf is in normal form. + +'solve_='(Nf) :- + deref(Nf,Nfd), % dereferences and turns Nf into solvable form Nfd + solve(Nfd). + +% 'solve_=\\='(Nf) +% +% Solves linear inequality Nf =\= 0 where Nf is in normal form. + +'solve_=\\='(Nf) :- + deref(Nf,Lind), % dereferences and turns Nf into solvable form Lind + Lind = [Inhom,_|Hom], + ( + Hom = [] -> % Lind = Inhom => check Inhom =\= 0 + + ( % Inhom =\= 0 + Inhom < -1e-010 -> + + true + ; + Inhom > 1e-010 + ) + ; + Hom = [_|_] -> + + var_with_def_intern(t_none,Nz,Lind,0), % make new variable Nz = Lind + get_attr(Nz,itf3,(Ty,St,Li,Or,Cl,Fo,_,RAtt)), + put_attr(Nz,itf3,(Ty,St,Li,Or,Cl,Fo,nonzero,RAtt)) % make Nz nonzero + ). + +% 'solve_<'(Nf) +% +% Solves linear inequality Nf < 0 where Nf is in normal form. + +'solve_<'(Nf) :- + split(Nf,H,I), + ineq(H,I,Nf,strict). + +% 'solve_=<'(Nf) +% +% Solves linear inequality Nf =< 0 where Nf is in normal form. + +'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) +% +% Minimizes the linear expression Lin. It does so by making a new +% variable Dep and minimizes its value. + +minimize_lin(Lin) :- + deref(Lin,Lind), + var_with_def_intern(t_none,Dep,Lind,0), + 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 until Expression becomes linear, Nf contains linear Expression in normal form + 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,0), % make new variable Dep = Lind + determine_active_dec(Lind), % minimizes Lind + iterate_dec(Dep,Inf), + vertex_value(Vector,Values), + ( + blackboard_inf(_) -> + + retract(blackboard_inf(_)), + assert(blackboard_inf([Inf|Values])) + ; + assert(blackboard_inf([Inf|Values])) + ), + fail. +inf_lin(_,Infimum,_,Vertex) :- + retract(blackboard_inf(L)), + assign([Infimum|Vertex],L). + +% assign(L1,L2) +% +% The elements of L1 are pairwise assigned to the elements of L2 +% by means of asserting {X =:= Y} where X is an element of L1 and Y +% is the corresponding element of L2. + +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 +% + + +% iterate_dec(OptVar,Opt) +% +% Decreases the bound on the variables of the linear equation of OptVar as much +% as possible and returns the resulting optimal bound in Opt. Fails if for some +% variable, a status of unlimited is found. + +iterate_dec(OptVar,Opt) :- + get_attr(OptVar,itf3,(_,_,lin(Lin),_)), + Lin = [I,R|H], + dec_step(H,Status), + ( + Status = applied -> + + iterate_dec(OptVar,Opt) + ; + Status = optimum -> + + Opt is R+I + ). + +% iterate_inc(OptVar,Opt) +% +% Increases the bound on the variables of the linear equation of OptVar as much +% as possible and returns the resulting optimal bound in Opt. Fails if for some +% variable, a status of unlimited is found. + +iterate_inc(OptVar,Opt) :- + get_attr(OptVar,itf3,(_,_,lin(Lin),_)), + Lin = [I,R|H], + inc_step(H,Status), + ( + Status = applied -> + + iterate_inc(OptVar,Opt) + ; + Status = optimum -> + + Opt is R+I + ). + +% +% 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([l(V*K,OrdV)|Vs],Status) :- + get_attr(V,itf3,(type(W),_,_,_,class(Class),_)), + ( + W = t_U(U) -> + + ( + K > 1e-010 -> % K > 0 + + ( + lb(Class,OrdV,Vub-Vb-_) -> % found a lower bound + + 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) -> + + ( + K > 1e-010 -> % K > 0 + + Status = applied, + Init is L-U, + class_basis(Class,Deps), + lb(Deps,OrdV,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) -> + + ( + K < -1e-010 -> % K < 0 + + ( + ub(Class,OrdV,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) -> + + ( + K < -1e-010 -> % K < 0 + + Status = applied, + Init is U-L, + class_basis(Class,Deps), + ub(Deps,OrdV,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). % if status has not been set yet: no changes +inc_step([l(V*K,OrdV)|Vs],Status) :- + get_attr(V,itf3,(type(W),_,_,_,class(Class),_)), + ( + W = t_U(U) -> % active upperbound + + ( + K < -1e-010 -> % K < 0: try to find lower bound for V + + ( + lb(Class,OrdV,Vub-Vb-_) -> + + Status = applied, + pivot_a(Vub,V,Vb,t_u(U)) + ; + Status = unlimited(V,t_u(U)) + ) + ; % K > 0: no way to increment + inc_step(Vs,Status) + ) + ; + W = t_lU(L,U) -> % active upperbound, passive lowerbound + + ( + K < -1e-010 -> % K < 0, try to find lower bound + + Status = applied, + Init is L-U, + class_basis(Class,Deps), + lb(Deps,OrdV,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) -> + + ( + K > 1e-010 -> % K > 0 + + ( + ub(Class,OrdV,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) -> + + ( + K > 1e-010 -> % K > 0 + + Status = applied, + Init is U-L, + class_basis(Class,Deps), + ub(Deps,OrdV,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) + ). + +% ------------------------- 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(Class,OrdX,Ub) +% +% See lb/3: this is similar + +ub(Class,OrdX,Ub) :- + class_basis(Class,Deps), + ub_first(Deps,OrdX,Ub). + +% ub_first(Deps,X,Dep-W-Ub) +% +% Finds the tightest upperbound for variable X from the linear equations of +% basis variables Deps, and puts the resulting bound in Ub. Dep is the basis +% variable that generates the bound, and W is bound of that variable that has +% to be activated to achieve this. + +ub_first([Dep|Deps],OrdX,Tightest) :- + ( + get_attr(Dep,itf3,(type(Type),_,lin(Lin),_)), + ub_inner(Type,OrdX,Lin,W,Ub), + Ub > -1e-010 -> % Ub >= 0 + + ub(Deps,OrdX,Dep-W-Ub,Tightest) + ; + ub_first(Deps,OrdX,Tightest) + ). + +% ub(Deps,OrdX,TightestIn,TightestOut) +% +% See lb/4: this is similar + +ub([],_,T0,T0). +ub([Dep|Deps],OrdX,T0,T1) :- + ( + get_attr(Dep,itf3,(type(Type),_,lin(Lin),_)), + ub_inner(Type,OrdX,Lin,W,Ub), + T0 = _-Ubb, + Ub - Ubb < -1e-010, % Ub < Ubb: tighter upper bound is a smaller one + Ub > -1e-010 -> % Ub >= 0: upperbound should be larger than 0; rare failure + + ub(Deps,OrdX,Dep-W-Ub,T1) % tighter bound, use new bound + ; + ub(Deps,OrdX,T0,T1) % no tighter bound, keep current one + ). + +% ub_inner(Type,OrdX,Lin,W,Ub) +% +% See lb_inner/5: this is similar + +ub_inner(t_l(L),OrdX,Lin,t_L(L),Ub) :- + nf_rhs_x(Lin,OrdX,Rhs,K), % Rhs is right hand side of lin. eq. Lin containing term X*K + K < -1e-010, % K < 0 + Ub is (L-Rhs)/K. +ub_inner(t_u(U),OrdX,Lin,t_U(U),Ub) :- + nf_rhs_x(Lin,OrdX,Rhs,K), + K > 1e-010, % K > 0 + Ub is (U-Rhs)/K. +ub_inner(t_lu(L,U),OrdX,Lin,W,Ub) :- + nf_rhs_x(Lin,OrdX,Rhs,K), + ( + K < -1e-010 -> % K < 0, use lowerbound + + W = t_Lu(L,U), + Ub = (L-Rhs)/K + + ; + K > 1e-010 -> % K > 0, use upperbound + + W = t_lU(L,U), + Ub = (U-Rhs)/K + ). + +% lb(Class,OrdX,Lb) +% +% Returns in Lb how much we can lower the upperbound of X without violating +% a bound of the basisvariables. +% Lb has the form Dep-W-Lb with Dep the variable whose bound is violated when +% lowering the bound for X more, W the actual bound that has to be activated and +% Lb the amount that the upperbound can be lowered. +% X has ordering OrdX and class Class. + +lb(Class,OrdX,Lb) :- + class_basis(Class,Deps), + lb_first(Deps,OrdX,Lb). + +% lb_first(Deps,OrdX,Tightest) +% +% Returns in Tightest how much we can lower the upperbound of X without violating +% a bound of Deps. +% Tightest has the form Dep-W-Lb with Dep the variable whose bound is violated when +% lowering the bound for X more, W the actual bound that has to be activated and +% Lb the amount that the upperbound can be lowered. X has ordering attribute OrdX. + +lb_first([Dep|Deps],OrdX,Tightest) :- + ( + get_attr(Dep,itf3,(type(Type),_,lin(Lin),_)), + lb_inner(Type,OrdX,Lin,W,Lb), + Lb < 1e-010 -> % Lb =< 0: Lb > 0 means a violated bound + + lb(Deps,OrdX,Dep-W-Lb,Tightest) + ; + lb_first(Deps,OrdX,Tightest) + ). + +% lb(Deps,OrdX,TightestIn,TightestOut) +% +% See lb_first/3: this one does the same thing, but is used for the steps after the first one +% and remembers the tightest bound so far. + +lb([],_,T0,T0). +lb([Dep|Deps],OrdX,T0,T1) :- + ( + get_attr(Dep,itf3,(type(Type),_,lin(Lin),_)), + lb_inner(Type,OrdX,Lin,W,Lb), + T0 = _-Lbb, + Lb - Lbb > 1e-010, % Lb > Lbb: choose the least lowering, others might violate bounds + Lb < 1e-010 -> % Lb =< 0: violation of a bound (without lowering) + + lb(Deps,OrdX,Dep-W-Lb,T1) + ; + lb(Deps,OrdX,T0,T1) + ). + +% lb_inner(Type,X,Lin,W,Lb) +% +% Returns in Lb how much lower we can make X without violating a bound +% by using the linear equation Lin of basis variable B which has type +% Type and which has to activate a bound (type W) to do so. +% +% E.g. when B has a lowerbound L, then L should always be smaller than I + R. +% So a lowerbound of X (which has scalar K in Lin), could be at most (L-(I+R))/K +% lower than its upperbound (if K is positive). +% Also note that Lb should always be smaller than 0, otherwise the row isn't feasible. +% X has ordering attribute OrdX. + +lb_inner(t_l(L),OrdX,Lin,t_L(L),Lb) :- + nf_rhs_x(Lin,OrdX,Rhs,K), % if linear equation Lin contains the term X*K, + % Rhs is the right hand side of that equation + K > 1e-010, % K > 0 + Lb is (L-Rhs)/K. +lb_inner(t_u(U),OrdX,Lin,t_U(U),Lb) :- + nf_rhs_x(Lin,OrdX,Rhs,K), + K < -1e-010, % K < 0 + Lb is (U-Rhs)/K. +lb_inner(t_lu(L,U),OrdX,Lin,W,Lb) :- + nf_rhs_x(Lin,OrdX,Rhs,K), + ( + K < -1e-010 -> + + W = t_lU(L,U), + Lb is (U-Rhs)/K + ; + K > 1e-010 -> + + W = t_Lu(L,U), + Lb is (L-Rhs)/K + ). + + +% ---------------------------------- 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) :- + Lin = [I,_|H], + solve(H,Lin,I,Bindings,[]), + export_binding(Bindings). + +% solve(Hom,Lin,I,Bind,BindT) +% +% Solves a linear equation Lin = [I,_|H] = 0 and exports the generated bindings. + +solve([],_,I,Bind0,Bind0) :- + !, + I >= -1e-010, % I =:= 0: redundant or trivially unsat + I =< 1e-010. +solve(H,Lin,_,Bind0,BindT) :- + sd(H,[],ClassesUniq,9-9-0,Category-Selected-_,NV,NVT), + get_attr(Selected,itf3,(_,_,_,order(Ord),_)), + isolate(Ord,Lin,Lin1), % Lin = 0 => Selected = Lin1 + ( + Category = 1 -> % classless variable, no bounds + + get_attr(Selected,itf3,(Ty,St,_,RAtt)), + put_attr(Selected,itf3,(Ty,St,lin(Lin1),RAtt)), + Lin1 = [Inhom,_|Hom], + bs_collect_binding(Hom,Selected,Inhom,Bind0,BindT), + eq_classes(NV,NVT,ClassesUniq) + ; + Category = 2 -> % class variable, no bounds + + get_attr(Selected,itf3,(_,_,_,_,class(NewC),_)), + class_allvars(NewC,Deps), + ( + ClassesUniq = [_] -> % rank increasing + + bs_collect_bindings(Deps,Ord,Lin1,Bind0,BindT) + ; + Bind0 = BindT, + bs(Deps,Ord,Lin1) + ), + eq_classes(NV,NVT,ClassesUniq) + ; + Category = 3 -> % classless variable, all variables in Lin and Selected are bounded + + get_attr(Selected,itf3,(type(Type),St,_,RAtt)), + put_attr(Selected,itf3,(type(Type),St,lin(Lin1),RAtt)), + deactivate_bound(Type,Selected), + eq_classes(NV,NVT,ClassesUniq), + basis_add(Selected,Basis), + undet_active(Lin1), % we can't tell which bound will likely be a problem at this point + Lin1 = [Inhom,_|Hom], + bs_collect_binding(Hom,Selected,Inhom,Bind0,Bind1), % only if Hom = [] + rcbl(Basis,Bind1,BindT) % reconsider entire basis + ; + Category = 4 -> % class variable, all variables in Lin and Selected are bounded + + get_attr(Selected,itf3,(type(Type),_,_,_,class(NewC),_)), + class_allvars(NewC,Deps), + ( + ClassesUniq = [_] -> % rank increasing + + bs_collect_bindings(Deps,Ord,Lin1,Bind0,Bind1) + ; + Bind0 = Bind1, + bs(Deps,Ord,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(H,Lin,I,X,[Bind|BindT],BindT) +% +% + +solve_x(Lin,X) :- + Lin = [I,_|H], + solve_x(H,Lin,I,X,Bindings,[]), + export_binding(Bindings). + +solve_x([],_,I,_,Bind0,Bind0) :- + !, + I >= -1e-010, % I =:= 0: redundant or trivially unsat + I =< 1e-010. +solve_x(H,Lin,_,X,Bind0,BindT) :- + sd(H,[],ClassesUniq,9-9-0,_,NV,NVT), + get_attr(X,itf3,(_,_,_,order(OrdX),_)), + isolate(OrdX,Lin,Lin1), + ( + get_attr(X,itf3,(_,_,_,_,class(NewC),_)) -> + + class_allvars(NewC,Deps), + ( + ClassesUniq = [_] -> % rank increasing + + bs_collect_bindings(Deps,OrdX,Lin1,Bind0,BindT) + ; + Bind0 = BindT, + bs(Deps,OrdX,Lin1) + ), + eq_classes(NV,NVT,ClassesUniq) + ; + get_attr(X,itf3,(Ty,St,_,RAtt)), + put_attr(X,itf3,(Ty,St,lin(Lin1),RAtt)), + Lin1 = [Inhom,_|Hom], + bs_collect_binding(Hom,X,Inhom,Bind0,BindT), + eq_classes(NV,NVT,ClassesUniq) + ). + +% solve_ord_x(Lin,OrdX,ClassX) +% +% Does the same thing as solve_x/2, but has the ordering of X and its class as input +% This also means that X has a class which is not sure in solve_x/2 + +solve_ord_x(Lin,OrdX,ClassX) :- + Lin = [I,_|H], + solve_ord_x(H,Lin,I,OrdX,ClassX,Bindings,[]), + export_binding(Bindings). + +solve_ord_x([],_,I,_,_,Bind0,Bind0) :- + I >= -1e-010, % I =:= 0 + I =< 1e-010. +solve_ord_x([_|_],Lin,_,OrdX,ClassX,Bind0,BindT) :- + isolate(OrdX,Lin,Lin1), + Lin1 = [_,_|H1], + sd(H1,[],ClassesUniq1,9-9-0,_,NV,NVT), % do sd on Lin without X, then add class of X + ord_add_element(ClassesUniq1,ClassX,ClassesUniq), + class_allvars(ClassX,Deps), + ( + ClassesUniq = [_] -> % rank increasing + + bs_collect_bindings(Deps,OrdX,Lin1,Bind0,BindT) + ; + Bind0 = BindT, + bs(Deps,OrdX,Lin1) + ), + eq_classes(NV,NVT,ClassesUniq). + + +% sd(H,[],ClassesUniq,9-9-0,Category-Selected-_,NV,NVT) + +% sd(Hom,ClassesIn,ClassesOut,PreferenceIn,PreferenceOut,[NV|NVTail],NVTail) +% +% ClassesOut is a sorted list of the different classes that are either in ClassesIn or +% that are the classes of the variables in Hom. Variables that do not belong to a class +% yet, are put in the difference list NV. + +sd([],Class0,Class0,Preference0,Preference0,NV0,NV0). +sd([l(X*K,_)|Xs],Class0,ClassN,Preference0,PreferenceN,NV0,NVt) :- + ( + get_attr(X,itf3,(_,_,_,_,class(Xc),_)) -> % old: has class + + NV0 = NV1, + ord_add_element(Class0,Xc,Class1), + ( + get_attr(X,itf3,(type(t_none),_)) -> + + preference(Preference0,2-X-K,Preference1) % has class, no bounds => category 2 + ; + preference(Preference0,4-X-K,Preference1) % has class, is bounded => category 4 + ) + ; % new: has no class + Class1 = Class0, + 'C'(NV0,X,NV1), % X has no class yet, add to list of new variables + ( + get_attr(X,itf3,(type(t_none),_)) -> + + preference(Preference0,1-X-K,Preference1) % no class, no bounds => category 1 + ; + preference(Preference0,3-X-K,Preference1) % no class, is bounded => category 3 + ) + ), + sd(Xs,Class1,ClassN,Preference1,PreferenceN,NV1,NVt). + +% +% A is best sofar, B is current +% smallest prefered +preference(A,B,Pref) :- + A = Px-_-_, + B = Py-_-_, + compare(Rel,Px,Py), + ( + Rel = (=) -> + + Pref = B % ( abs(Ka)= Pref=A ; Pref=B ) + ; + Rel = (<) -> + + Pref = A + ; + Rel = (>) -> + + Pref = B + ). + +% eq_classes(NV,NVTail,Cs) +% +% Attaches all classless variables NV to a new class and equates all other classes +% with this class. The equate operation only happens after attach_class because the +% unification of classes can bind the tail of the AllVars attribute to a nonvar and +% then the attach_class operation wouldn't work. + +eq_classes(NV,_,Cs) :- + var(NV), + !, + equate(Cs,_). +eq_classes(NV,NVT,Cs) :- + class_new(Su,NV,NVT,[]), % make a new class Su with NV as the variables + attach_class(NV,Su), % attach the variables NV to 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), % Tail + !. +attach_class([X|Xs],Class) :- + get_attr(X,itf3,(Ty,St,Li,Or,_,RAtt)), + put_attr(X,itf3,(Ty,St,Li,Or,class(Class),RAtt)), + attach_class(Xs,Class). + +% unconstrained(Lin,Uc,Kuc,Rest) +% +% Finds an unconstrained variable Uc (type(t_none)) in Lin with scalar Kuc and removes it +% from Lin to return Rest. + +unconstrained(Lin,Uc,Kuc,Rest) :- + Lin = [_,_|H], + sd(H,[],_,9-9-0,Category-Uc-_,_,_), + Category =< 2, + get_attr(Uc,itf3,(_,_,_,order(OrdUc),_)), + delete_factor(OrdUc,Lin,Rest,Kuc). + + +% +% point the vars in Lin into the same equivalence class +% maybe join some global data +% +same_class([],_). +same_class([l(X*_,_)|Xs],Class) :- + get_or_add_class(X,Class), + same_class(Xs,Class). + +% get_or_add_class(X,Class) +% +% Returns in Class the class of X if X has one, or a new class where X now belongs to if X didn't have one. + +get_or_add_class(X,Class) :- + get_attr(X,itf3,(_,_,_,_,class(ClassX),_)), + !, + ClassX = Class. % explicit =/2 because of cut +get_or_add_class(X,Class) :- + get_attr(X,itf3,(Ty,St,Li,Or,_,RAtt)), + put_attr(X,itf3,(Ty,St,Li,Or,class(Class),RAtt)), + class_new(Class,[X|Tail],Tail,[]). % initial class atts + +% allvars(X,Allvars) +% +% Allvars is a list of all variables in the class to which X belongs. + +allvars(X,Allvars) :- + get_attr(X,itf3,(_,_,_,_,class(C),_)), + class_allvars(C,Allvars). + +% deactivate_bound(Type,Variable) +% +% The Type of the variable is changed to reflect the deactivation of its bounds. +% t_L(_) becomes t_l(_), t_lU(_,_) becomes t_lu(_,_) and so on. + +deactivate_bound(t_l(_),_). +deactivate_bound(t_u(_),_). +deactivate_bound(t_lu(_,_),_). +deactivate_bound(t_L(L),X) :- + get_attr(X,itf3,(_,RAtt)), + put_attr(X,itf3,(type(t_l(L)),RAtt)). +deactivate_bound(t_Lu(L,U),X) :- + get_attr(X,itf3,(_,RAtt)), + put_attr(X,itf3,(type(t_lu(L,U)),RAtt)). +deactivate_bound(t_U(U),X) :- + get_attr(X,itf3,(_,RAtt)), + put_attr(X,itf3,(type(t_u(U)),RAtt)). +deactivate_bound(t_lU(L,U),X) :- + get_attr(X,itf3,(_,RAtt)), + put_attr(X,itf3,(type(t_lu(L,U)),RAtt)). + +% intro_at(X,Value,Type) +% +% Variable X gets new type Type which reflects the activation of a bound with value +% Value. In the linear equations of all the variables belonging to the same class as +% X, X is substituted by [0,Value,X] to reflect the new active bound. + +intro_at(X,Value,Type) :- + get_attr(X,itf3,(_,St,Li,order(Ord),class(Class),RAtt)), + put_attr(X,itf3,(type(Type),St,Li,order(Ord),class(Class),RAtt)), + ( + Value >= -1e-010, % Value =:= 0 + Value =< 1e-010 -> + + true + ; + backsubst_delta(Class,Ord,X,Value) + ). + +% undet_active(Lin) +% +% For each variable in the homogene part of Lin, a bound is activated +% if an inactive bound exists. (t_l(L) becomes t_L(L) and so on) + +undet_active(Lin) :- + Lin = [_,_|Lin1], + undet_active_h(Lin1). + +% undet_active_h(Hom) +% +% For each variable in homogene part Hom, a bound is activated if an +% inactive bound exists (t_l(L) becomes t_L(L) and so on) + +undet_active_h([]). +undet_active_h([l(X*_,_)|Xs]) :- + get_attr(X,itf3,(type(Type),_)), + undet_active(Type,X), + undet_active_h(Xs). + +% undet_active(Type,Var) +% +% An inactive bound of Var is activated if such exists +% t_lu(L,U) is arbitrarily chosen to become t_Lu(L,U) + +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) +% +% Activates inactive bounds on the variables of Lin if such bounds exist. +% If the type of a variable is t_none, this fails. This version is aimed +% to make the R component of Lin as small as possible in order not to violate +% an upperbound (see reconsider/1) + +determine_active_dec(Lin) :- + Lin = [_,_|Lin1], + determine_active(Lin1,-1.0). + +% determine_active_inc(Lin) +% +% Activates inactive bounds on the variables of Lin if such bounds exist. +% If the type of a variable is t_none, this fails. This version is aimed +% to make the R component of Lin as large as possible in order not to violate +% a lowerbound (see reconsider/1) + +determine_active_inc(Lin) :- + Lin = [_,_|Lin1], + determine_active(Lin1,1.0). + +% determine_active(Hom,S) +% +% For each variable in Hom, activates its bound if it is not yet activated. +% For the case of t_lu(_,_) the lower or upper bound is activated depending on K and S: +% If sign of K*S is negative, then lowerbound, otherwise upperbound. + +determine_active([],_). +determine_active([l(X*K,_)|Xs],S) :- + get_attr(X,itf3,(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) :- + TestKs is K*S, + ( + TestKs < -1e-010 -> + + intro_at(X,L,t_Lu(L,U)) % K*S < 0 + ; + TestKs > 1e-010 -> + + intro_at(X,U,t_lU(L,U)) + ). + +% +% Careful when an indep turns into t_none !!! +% +detach_bounds(V) :- + get_attr(V,itf3,(_,_,lin(Lin),order(OrdV),class(Class),RAtt)), + put_attr(V,itf3,(type(t_none),strictness(0),lin(Lin),order(OrdV),class(Class),RAtt)), + ( + indep(Lin,OrdV) -> + + ( + ub(Class,OrdV,Vub-Vb-_) -> % exchange against thightest + + class_basis_drop(Class,Vub), + pivot(Vub,Class,OrdV,Vb) + ; + lb(Class,OrdV,Vlb-Vb-_) -> + + class_basis_drop(Class,Vlb), + pivot(Vlb,Class,OrdV,Vb) + ; + true + ) + ; + class_basis_drop(Class,V) + ). + +detach_bounds_vlv(OrdV,Lin,Class,Var,NewLin) :- + ( + indep(Lin,OrdV) -> + + ( + ub(Class,OrdV,Vub-Vb-_) -> + + % in verify_lin, class might contain two occurrences of Var, but it + % doesn't matter which one we delete + class_basis_drop(Class,Var), + pivot_vlv(Vub,Class,OrdV,Vb,NewLin) + ; + lb(Class,OrdV,Vlb-Vb-_) -> + + class_basis_drop(Class,Var), + pivot_vlv(Vlb,Class,OrdV,Vb,NewLin) + ; + NewLin = Lin + ) + ; + NewLin = Lin, + class_basis_drop(Class,Var) + ). + + + +% ----------------------------- manipulate the basis --------------------------- % + +% basis_drop(X) +% +% Removes X from the basis of the class to which X belongs. + +basis_drop(X) :- + get_attr(X,itf3,(_,_,_,_,class(Cv),_)), + class_basis_drop(Cv,X). + +% basis(X,Basis) +% +% Basis is the basis of the class to which X belongs. + +basis(X,Basis) :- + get_attr(X,itf3,(_,_,_,_,class(Cv),_)), + class_basis(Cv,Basis). + +% basis_add(X,NewBasis) +% +% NewBasis is the result of adding X to the basis of the class to which X belongs. + +basis_add(X,NewBasis) :- + get_attr(X,itf3,(_,_,_,_,class(Cv),_)), + class_basis_add(Cv,X,NewBasis). + +% basis_pivot(Leave,Enter) +% +% Removes Leave from the basis of the class to which it belongs, and adds Enter to that basis. + +basis_pivot(Leave,Enter) :- + get_attr(Leave,itf3,(_,_,_,_,class(Cv),_)), + class_basis_pivot(Cv,Enter,Leave). + +% ----------------------------------- pivot ------------------------------------ % + +% pivot(Dep,Indep) +% +% The linear equation of variable Dep, is transformed into one of variable Indep, +% containing Dep. Then, all occurrences of Indep in linear equations are substituted +% by this new definition + +% +% Pivot ignoring rhs and active states +% + +pivot(Dep,Indep) :- + get_attr(Dep,itf3,(_,_,lin(H),order(OrdDep),_)), + get_attr(Indep,itf3,(_,_,_,order(Ord),class(Class),_)), + delete_factor(Ord,H,H0,Coeff), + K is -1.0/Coeff, + add_linear_ff(H0,K,[0.0,0.0,l(Dep* -1.0,OrdDep)],K,Lin), + backsubst(Class,Ord,Lin). + +% pivot_a(Dep,Indep,IndepT,DepT) +% +% Removes Dep from the basis, puts Indep in, and pivots the equation of +% Dep to become one of Indep. The type of Dep becomes DepT (which means +% it gets deactivated), the type of Indep becomes IndepT (which means it +% gets activated) + + +pivot_a(Dep,Indep,Vb,Wd) :- + basis_pivot(Dep,Indep), + get_attr(Indep,itf3,(_,_,_,order(Ord),class(Class),_)), + pivot(Dep,Class,Ord,Vb), + get_attr(Indep,itf3,(_,RAtt)), + put_attr(Indep,itf3,(type(Wd),RAtt)). + +pivot_b(Vub,V,Vb,Wd) :- + ( + Vub == V -> + + get_attr(V,itf3,(_,St,Li,order(Ord),class(Class),RAtt)), + put_attr(V,itf3,(type(Vb),St,Li,order(Ord),class(Class),RAtt)), + pivot_b_delta(Vb,Delta), % nonzero(Delta) + backsubst_delta(Class,Ord,V,Delta) + ; + pivot_a(Vub,V,Vb,Wd) + ). + +pivot_b_delta(t_Lu(L,U),Delta) :- Delta is L-U. +pivot_b_delta(t_lU(L,U),Delta) :- Delta is U-L. + +% select_active_bound(Type,Bound) +% +% Returns the bound that is active in Type (if such exists, 0 otherwise) + +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,0.0). +% +% for project.pl +% +select_active_bound(t_l(_),0.0). +select_active_bound(t_u(_),0.0). +select_active_bound(t_lu(_,_),0.0). + + +% pivot(Dep,Class,IndepOrd,IndAct) +% +% See pivot/2. +% In addition, variable Indep with ordering IndepOrd has an active bound IndAct. + +% +% +% Pivot taking care of rhs and active states +% +pivot(Dep,Class,IndepOrd,IndAct) :- + get_attr(Dep,itf3,(_,St,lin(H),order(DepOrd),RAtt)), + put_attr(Dep,itf3,(type(IndAct),St,lin(H),order(DepOrd),RAtt)), + select_active_bound(IndAct,Abv), % Dep or Indep + delete_factor(IndepOrd,H,H0,Coeff), % Dep = ... + Coeff*Indep + ... + K is -1.0/Coeff, + Abvm is -Abv, + add_linear_ff(H0,K,[0.0,Abvm,l(Dep* -1.0,DepOrd)],K,Lin), % Indep = -1/Coeff*... + 1/Coeff*Dep + backsubst(Class,IndepOrd,Lin). + +pivot_vlv(Dep,Class,IndepOrd,IndAct,Lin) :- + get_attr(Dep,itf3,(_,St,lin(H),order(DepOrd),RAtt)), + put_attr(Dep,itf3,(type(IndAct),St,lin(H),order(DepOrd),RAtt)), + select_active_bound(IndAct,Abv), % Dep or Indep + delete_factor(IndepOrd,H,H0,Coeff), % Dep = ... + Coeff*Indep + ... + K is -1.0/Coeff, + Abvm is -Abv, + add_linear_ff(H0,K,[0.0,Abvm,l(Dep* -1.0,DepOrd)],K,Lin), % Indep = -1/Coeff*... + 1/Coeff*Dep + backsubst(Class,IndepOrd,Lin). + + + +% backsubst_delta(Class,OrdX,X,Delta) +% +% X with ordering attribute OrdX, is substituted in all linear equations of variables in the class Class, +% by linear equation [0,Delta,l(X*1,OrdX)]. This reflects the activation of a bound. + +backsubst_delta(Class,OrdX,X,Delta) :- + backsubst(Class,OrdX,[0.0,Delta,l(X*1.0,OrdX)]). + +% backsubst(Class,OrdX,Lin) +% +% X with ordering OrdX is substituted in all linear equations of variables in the class Class, +% by linear equation Lin + +backsubst(Class,OrdX,Lin) :- + class_allvars(Class,Allvars), + bs(Allvars,OrdX,Lin). + +% bs(Vars,OrdV,Lin) +% +% In all linear equations of the variables Vars, variable V with ordering attribute OrdV +% is substituted by linear equation Lin. +% +% valid if nothing will go ground +% + +bs(Xs,_,_) :- + var(Xs), + !. +bs([X|Xs],OrdV,Lin) :- + ( + get_attr(X,itf3,(Ty,St,lin(LinX),RAtt)), + nf_substitute(OrdV,Lin,LinX,LinX1) -> + + put_attr(X,itf3,(Ty,St,lin(LinX1),RAtt)), + bs(Xs,OrdV,Lin) + ; + bs(Xs,OrdV,Lin) + ). + + +% +% rank increasing backsubstitution +% + +% bs_collect_bindings(Deps,SelectedOrd,Lin,Bind,BindT) +% +% Collects bindings (of the form [X-I] where X = I is the binding) by substituting +% Selected in all linear equations of the variables Deps (which are of the same class), +% by Lin. Selected has ordering attribute SelectedOrd. +% +% E.g. when V = 2X + 3Y + 4, X = 3V + 2Z and Y = 4X + 3 +% we can substitute V in the linear equation of X: X = 6X + 9Y + 2Z + 12 +% we can't substitute V in the linear equation of Y of course. + +bs_collect_bindings(Xs,_,_,Bind0,BindT) :- + var(Xs), + !, + Bind0 = BindT. +bs_collect_bindings([X|Xs],OrdV,Lin,Bind0,BindT) :- + ( + get_attr(X,itf3,(Ty,St,lin(LinX),RAtt)), + nf_substitute(OrdV,Lin,LinX,LinX1) -> + + put_attr(X,itf3,(Ty,St,lin(LinX1),RAtt)), + LinX1 = [Inhom,_|Hom], + bs_collect_binding(Hom,X,Inhom,Bind0,Bind1), + bs_collect_bindings(Xs,OrdV,Lin,Bind1,BindT) + ; + bs_collect_bindings(Xs,OrdV,Lin,Bind0,BindT) + ). + +% bs_collect_binding(Hom,Selected,Inhom,Bind,BindT) +% +% Collects binding following from Selected = Hom + Inhom. +% If Hom = [], returns the binding Selected-Inhom (=0) +% +bs_collect_binding([],X,Inhom) --> [X-Inhom]. +bs_collect_binding([_|_],_,_) --> []. + +% +% reconsider the basis +% + +% rcbl(Basis,Bind,BindT) +% +% + +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 +% +% +% Idea: consider a variable V with linear equation Lin. +% When a bound on a variable X of Lin gets activated, its value, multiplied with the scalar of X, +% is added to the R component of Lin. +% When we consider the lowerbound of V, it must be smaller than R + I, since R contains at best the +% lowerbounds of the variables in Lin (but could contain upperbounds, which are of course larger). +% So checking this can show the violation of a bound of V. A similar case works for the upperbound. + +rcb(X,Status,Violated) :- + get_attr(X,itf3,(type(Type),_,lin(Lin),_)), + Lin = [I,R|H], + ( + % case 1: lowerbound: R + I should always be larger than the lowerbound + Type = t_l(L) -> + + R + I - L < 1e-010, % R + I =< L + Violated = l(L), + inc_step(H,Status) + ; + % case 2: upperbound: R + I should always be smaller than the upperbound + Type = t_u(U) -> + + R + I - U > -1e-010, % R+I >= U + Violated = u(U), + dec_step(H,Status) + + ; + % case 3: check both + Type = t_lu(L,U) -> + + At is R + I, + ( + At - L < 1e-010 -> % At =< L + + Violated = l(L), + inc_step(H,Status) + ; + At - U > -1e-010 -> % At >= U + + Violated = u(U), + dec_step( H, Status) + ) + % don't care for other types + ). + +% rcbl_status(Status,X,Continuation,[Bind|BindT],BindT,Violated) +% +% + +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_attr(X,itf3,(type(Type),strictness(Strict),lin(Lin),_)), + Lin = [I,R|_], + Opt is R+I, + TestLO is L - Opt, + ( + TestLO < -1e-010 -> + + narrow_u(Type,X,Opt), % { X =< Opt } + rcbl(Continuation,B0,B1) + ; + TestLO > 1e-010 -> + + fail + ; + Strict /\ 2 =:= 0, % meets lower + Mop is -Opt, + normalize_scalar(Mop,MopN), + add_linear_11(MopN,Lin,Lin1), + Lin1 = [Inhom,_|Hom], + ( + Hom = [] -> + + rcbl(Continuation,B0,B1) % would not callback + ; + Hom = [_|_] -> + + solve(Hom,Lin1,Inhom,B0,B1) + ) + ). +rcbl_opt(u(U),X,Continuation,B0,B1) :- + get_attr(X,itf3,(type(Type),strictness(Strict),lin(Lin),_)), + Lin = [I,R|_], + Opt is R+I, + TestUO is U-Opt, + ( + TestUO < -1e-010 -> + + fail + ; + TestUO > 1e-010 -> + + narrow_l(Type,X,Opt), % { X >= Opt } + rcbl(Continuation,B0,B1) + ; + Strict /\ 1 =:= 0, % meets upper + Mop is -Opt, + normalize_scalar(Mop,MopN), + add_linear_11(MopN,Lin,Lin1), + Lin1 = [Inhom,_|Hom], + ( + Hom = [] -> + + rcbl(Continuation,B0,B1) % would not callback + ; + Hom = [_|_] -> + + solve(Hom,Lin1,Inhom,B0,B1) + ) + ). + +% +% Basis has already changed when this is called +% +rcbl_app(l(L),X,Continuation,B0,B1) :- + get_attr(X,itf3,(_,_,lin(Lin),_)), + Lin = [I,R|H], + ( + R + I - L > 1e-010 -> % R+I > L: within bound now + + rcbl(Continuation,B0,B1) + ; + inc_step(H,Status), + rcbl_status(Status,X,Continuation,B0,B1,l(L)) + ). +rcbl_app(u(U),X,Continuation,B0,B1) :- + get_attr(X,itf3,(_,_,lin(Lin),_)), + Lin = [I,R|H], + ( + R + I - U < -1e-010 -> % 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(Type,X,U) +% +% Narrows down the upperbound of X (type Type) to U. +% Fails if Type is not t_u(_) or t_lu(_) + +narrow_u(t_u(_),X,U) :- + get_attr(X,itf3,(_,RAtt)), + put_attr(X,itf3,(type(t_u(U)),RAtt)). +narrow_u(t_lu(L,_),X,U) :- + get_attr(X,itf3,(_,RAtt)), + put_attr(X,itf3,(type(t_lu(L,U)),RAtt)). + +% narrow_l(Type,X,L) +% +% Narrows down the lowerbound of X (type Type) to L. +% Fails if Type is not t_l(_) or t_lu(_) + +narrow_l( t_l(_), X, L) :- + get_attr(X,itf3,(_,RAtt)), + put_attr(X,itf3,(type(t_l(L)),RAtt)). + +narrow_l( t_lu(_,U), X, L) :- + get_attr(X,itf3,(_,RAtt)), + put_attr(X,itf3,(type(t_lu(L,U)),RAtt)). + + +% ----------------------------------- dump ------------------------------------- + +% dump_var(Type,Var,I,H,Dump,DumpTail) +% +% Returns in Dump a representation of the linear constraint on variable +% Var which has linear equation H + I and has type Type. + +dump_var(t_none,V,I,H) --> !, + ( + { + H=[l(W*K,_)], + V==W, + I >= -1e-010, % I=:=0 + I =< 1e-010, + TestK is K - 1.0, % K=:=1 + TestK =< 1e-010, + TestK >= -1e-010 + } -> % indep var + + [] + ; + { nf2sum(H,I,Sum) }, + [V = Sum] + ). +dump_var(t_L(L),V,I,H) --> + !, + dump_var(t_l(L),V,I,H). +% case lowerbound: V >= L or V > L +% say V >= L, and V = K*V1 + ... + I, then K*V1 + ... + I >= L +% and K*V1 + ... >= L-I and V1 + .../K = (L-I)/K +dump_var(t_l(L),V,I,H) --> + !, + { + H = [l(_*K,_)|_], % avoid 1 >= 0 + get_attr(V,itf3,(_,strictness(Strict),_)), + Sm is Strict /\ 2, + Kr is 1.0/K, + Li is Kr*(L-I), + mult_hom(H,Kr,H1), + nf2sum(H1,0.0,Sum), + ( + K > 1e-010 -> % 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 = [l(_*K,_)|_], % avoid 0 =< 1 + get_attr(V,itf3,(_,strictness(Strict),_)), + Sm is Strict /\ 1, + Kr is 1.0/K, + Ui is Kr*(U-I), + mult_hom(H,Kr,H1), + nf2sum(H1,0.0,Sum), + ( + K > 1e-010 -> % 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) --> % shouldn't happen + [ V:T:I+H ]. + +% dump_strict(FilteredStrictness,Nonstrict,Strict,Res) +% +% Unifies Res with either Nonstrict or Strict depending on FilteredStrictness. +% FilteredStrictness is the component of strictness related to the bound: 0 means nonstrict, +% 1 means strict upperbound, 2 means strict lowerbound, 3 is filtered out to either 1 or 2. + +dump_strict(0,Result,_,Result). +dump_strict(1,_,Result,Result). +dump_strict(2,_,Result,Result). + +% dump_nz(V,H,I,Dump,DumpTail) +% +% Returns in Dump a representation of the nonzero constraint of variable V which has linear +% equation H + I. + +dump_nz(_,H,I) --> + { + H = [l(_*K,_)|_], + Kr is 1.0/K, + I1 is -Kr*I, + mult_hom(H,Kr,H1), + nf2sum(H1,0.0,Sum) + }, + [ Sum =\= I1 ]. diff --git a/LGPL/clpr/clpr/class.pl b/LGPL/clpr/clpr/class.pl new file mode 100644 index 000000000..c5b8a9071 --- /dev/null +++ b/LGPL/clpr/clpr/class.pl @@ -0,0 +1,154 @@ +/* $Id: class.pl,v 1.1 2005-10-28 17:51:01 vsc Exp $ + + Part of CPL(R) (Constraint Logic Programming over Reals) + + Author: Leslie De Koninck + E-mail: Tom.Schrijvers@cs.kuleuven.ac.be + WWW: http://www.swi-prolog.org + http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 + Copyright (C): 2004, K.U. Leuven and + 1992-1995, Austrian Research Institute for + Artificial Intelligence (OFAI), + Vienna, Austria + + This software is part of Leslie De Koninck's master thesis, supervised + by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) + by Christian Holzbaur for SICStus Prolog and distributed under the + license details below with permission from all mentioned authors. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 + of the License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + As a special exception, if you link this library with other files, + compiled with a Free Software compiler, to produce an executable, this + library does not by itself cause the resulting executable to be covered + by the GNU General Public License. This exception does not however + invalidate any other reasons why the executable file might be covered by + the GNU General Public License. +*/ + + +:- module(class, + [ + class_allvars/2, + class_new/4, + class_drop/2, + class_basis/2, + class_basis_add/3, + class_basis_drop/2, + class_basis_pivot/3, + class_get_prio/2, + class_put_prio/2, + + ordering/1, + arrangement/2 + + ]). + +:- use_module(ordering, + [ + combine/3, + ordering/1, + arrangement/2 + ]). + +% called when two classes are unified: the allvars lists are appended to eachother, as well as the basis +% lists. +% +% note: La=[A,B,...,C|Lat], Lb=[D,E,...,F|Lbt], so new La = [A,B,...,C,D,E,...,F|Lbt] + +attr_unify_hook(class(La,Lat,ABasis,PrioA),Y) :- + !, + var(Y), + get_attr(Y,class,class(Lb,Lbt,BBasis,PrioB)), + Lat = Lb, + append(ABasis,BBasis,CBasis), + combine(PrioA,PrioB,PrioC), + put_attr(Y,class,class(La,Lbt,CBasis,PrioC)). +attr_unify_hook(_,_). + +class_new(Class,All,AllT,Basis) :- + put_attr(Su,class,class(All,AllT,Basis,[])), + Su = Class. + +class_get_prio(Class,Priority) :- get_attr(Class,class,class(_,_,_,Priority)). + +class_put_prio(Class,Priority) :- + get_attr(Class,class,class(All,AllT,Basis,_)), + put_attr(Class,class,class(All,AllT,Basis,Priority)). + +class_drop(Class,X) :- + get_attr(Class,class,class(Allvars,Tail,Basis,Priority)), + delete_first(Allvars,X,NewAllvars), + delete_first(Basis,X,NewBasis), + put_attr(Class,class,class(NewAllvars,Tail,NewBasis,Priority)). + +class_allvars(Class,All) :- get_attr(Class,class,class(All,_,_,_)). + +% class_basis(Class,Basis) +% +% Returns the basis of class Class. + +class_basis(Class,Basis) :- get_attr(Class,class,class(_,_,Basis,_)). + +% class_basis_add(Class,X,NewBasis) +% +% adds X in front of the basis and returns the new basis + +class_basis_add(Class,X,NewBasis) :- + NewBasis = [X|Basis], + get_attr(Class,class,class(All,AllT,Basis,Priority)), + put_attr(Class,class,class(All,AllT,NewBasis,Priority)). + +% class_basis_drop(Class,X) +% +% removes the first occurence of X from the basis (if exists) + +class_basis_drop(Class,X) :- + get_attr(Class,class,class(All,AllT,Basis0,Priority)), + delete_first(Basis0,X,Basis), + Basis0 \== Basis, % anything deleted ? + !, + put_attr(Class,class,class(All,AllT,Basis,Priority)). +class_basis_drop(_,_). + +% class_basis_pivot(Class,Enter,Leave) +% +% removes first occurence of Leave from the basis and adds Enter in front of the basis + +class_basis_pivot(Class,Enter,Leave) :- + get_attr(Class,class,class(All,AllT,Basis0,Priority)), + delete_first(Basis0,Leave,Basis1), + put_attr(Class,class,class(All,AllT,[Enter|Basis1],Priority)). + +% delete_first(Old,Element,New) +% +% removes the first occurence of Element from Old and returns the result in New +% +% note: test via syntactic equality, not unifiability + +delete_first(L,_,Res) :- + var(L), + !, + Res = L. +delete_first([],_,[]). +delete_first([Y|Ys],X,Res) :- + ( + X==Y -> + Res = Ys + ; + Res = [Y|Tail], + delete_first( Ys, X, Tail) + ). + diff --git a/LGPL/clpr/clpr/dump.pl b/LGPL/clpr/clpr/dump.pl new file mode 100644 index 000000000..0dbd5b021 --- /dev/null +++ b/LGPL/clpr/clpr/dump.pl @@ -0,0 +1,341 @@ +/* $Id: dump.pl,v 1.1 2005-10-28 17:51:01 vsc Exp $ + + Part of CPL(R) (Constraint Logic Programming over Reals) + + Author: Leslie De Koninck + E-mail: Tom.Schrijvers@cs.kuleuven.ac.be + WWW: http://www.swi-prolog.org + http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 + Copyright (C): 2004, K.U. Leuven and + 1992-1995, Austrian Research Institute for + Artificial Intelligence (OFAI), + Vienna, Austria + + This software is part of Leslie De Koninck's master thesis, supervised + by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) + by Christian Holzbaur for SICStus Prolog and distributed under the + license details below with permission from all mentioned authors. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 + of the License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + As a special exception, if you link this library with other files, + compiled with a Free Software compiler, to produce an executable, this + library does not by itself cause the resulting executable to be covered + by the GNU General Public License. This exception does not however + invalidate any other reasons why the executable file might be covered by + the GNU General Public License. +*/ + + +:- module(dump, + [ + dump/3, + projecting_assert/1 + ]). +:- use_module(class, + [ + class_allvars/2 + ]). +:- use_module(geler, + [ + collect_nonlin/3 + ]). +:- use_module(library(assoc), + [ + empty_assoc/1, + get_assoc/3, + put_assoc/4, + assoc_to_list/2 + ]). +:- use_module(itf3, + [ + dump_linear/4, + dump_nonzero/4 + ]). +:- use_module(project, + [ + project_attributes/2 + ]). +:- use_module(ordering, + [ ordering/1 + ]). +:- dynamic blackboard_copy/1. % replacement of bb_put and bb_delete by dynamic predicate + +this_linear_solver(clpr). + +% dump(Target,NewVars,Constraints) +% +% Returns in Constraints, the constraints that currently hold on Target where all variables +% in Target are copied to new variables in NewVars and the constraints are given on these new +% variables. In short, you can safely manipulate NewVars and Constraints without changing the +% constraints on Target. + +dump([],[],[]). +dump(Target,NewVars,Constraints) :- + ( + ( + proper_varlist(Target) -> + + true + ; + % Target is not a list of variables + raise_exception(instantiation_error(dump(Target,NewVars,Constraints),1)) + ), + ordering(Target), + related_linear_vars(Target,All), % All contains all variables of the classes of Target variables. + nonlin_crux(All,Nonlin), + project_attributes(Target,All), + related_linear_vars(Target,Again), % project drops/adds vars + all_attribute_goals(Again,Gs,Nonlin), + empty_assoc(D0), + mapping(Target,NewVars,D0,D1), % late (AVL suffers from put_atts) + copy(Gs,Copy,D1,_), % strip constraints + ( + blackboard_copy(_) -> + + retract(blackboard_copy(_)), + assert(blackboard_copy(NewVars/Copy)) + ; + assert(blackboard_copy(NewVars/Copy)) + ), + fail % undo projection + ; + retract(blackboard_copy(NewVars/Constraints)) % garbage collect + ). + +:- meta_predicate projecting_assert(:). + +projecting_assert(QClause) :- + strip_module(QClause, Module, Clause), % JW: SWI-Prolog not always qualifies the term! + copy_term(Clause,Copy,Constraints), + l2c(Constraints,Conj), % fails for [] + this_linear_solver(Sm), % proper module for {}/1 + !, + ( + Copy = (H:-B) -> % former rule + + Module:assert((H:-Sm:{Conj},B)) + ; % former fact + Module:assert((Copy:-Sm:{Conj})) + ). +projecting_assert(Clause) :- % not our business + assert(Clause). + +copy_term(Term,Copy,Constraints) :- + ( + term_variables(Term,Target), % get all variables in Term + related_linear_vars(Target,All), % get all variables of the classes of the variables in Term + nonlin_crux(All,Nonlin), % get a list of all the nonlinear goals of these variables + project_attributes(Target,All), + related_linear_vars(Target,Again), % project drops/adds vars + all_attribute_goals(Again,Gs,Nonlin), + empty_assoc(D0), + copy(Term/Gs,TmpCopy,D0,_), % strip constraints + ( + blackboard_copy(_) -> + + retract(blackboard_copy(_)), + assert(blackboard_copy(TmpCopy)) + ; + assert(blackboard_copy(TmpCopy)) + ), + fail % undo projection + ; + retract(blackboard_copy(Copy/Constraints)) % garbage collect + ). + +% l2c(Lst,Conj) +% +% converts a list to a round list: [a,b,c] -> (a,b,c) and [a] becomes a + +l2c([X|Xs],Conj) :- + ( + Xs = [] -> + + Conj = X + ; + Conj = (X,Xc), + l2c(Xs,Xc) + ). + +% proper_varlist(List) +% +% Returns whether Lst is a list of variables. +% First predicate is to avoid unification of a variable with a list. + +proper_varlist(X) :- + var(X), + !, + fail. +proper_varlist([]). +proper_varlist([X|Xs]) :- + var(X), + proper_varlist(Xs). + +% related_linear_vars(Vs,All) +% +% Generates a list of all variables that are in the classes of the variables in Vs. + +related_linear_vars(Vs,All) :- + empty_assoc(S0), + related_linear_sys(Vs,S0,Sys), + related_linear_vars(Sys,All,[]). + +% related_linear_sys(Vars,Assoc,List) +% +% Generates in List, a list of all to classes to which variables in Vars belong. +% Assoc should be an empty association list and is used internally. +% List contains elements of the form C-C where C is a class and both C's are equal. + +related_linear_sys([],S0,L0) :- assoc_to_list(S0,L0). +related_linear_sys([V|Vs],S0,S2) :- + ( + get_attr(V,itf3,(_,_,_,_,class(C),_)) -> + + put_assoc(C,S0,C,S1) + ; + S1 = S0 + ), + related_linear_sys(Vs,S1,S2). + +% related_linear_vars(Classes,[Vars|VarsTail],VarsTail) +% +% Generates a difference list of all variables in the classes in Classes. +% Classes contains elements of the form C-C where C is a class and both C's are equal. + +related_linear_vars([]) --> []. +related_linear_vars([S-_|Ss]) --> + { + class_allvars(S,Otl) + }, + cpvars(Otl), + related_linear_vars(Ss). + +% cpvars(Vars,Out,OutTail) +% +% Makes a new difference list of the difference list Vars. +% All nonvars are removed. + +cpvars(Xs) --> {var(Xs)}, !. +cpvars([X|Xs]) --> + ( + {var(X)} -> + + [X] + ; + [] + ), + cpvars(Xs). + +% nonlin_crux(All,Gss) +% +% Collects all pending non-linear constraints of variables in All. +% This marks all nonlinear goals of the variables as run and cannot +% be reversed manually. + +nonlin_crux(All,Gss) :- + collect_nonlin(All,Gs,[]), % collect the nonlinear goals of variables All + % this marks the goals as run and cannot be reversed manually + this_linear_solver(Solver), + nonlin_strip(Gs,Solver,Gss). + +% nonlin_strip(Gs,Solver,Res) +% +% Removes the goals from Gs that are not from solver Solver. + +nonlin_strip([],_,[]). +nonlin_strip([M:What|Gs],Solver,Res) :- + ( + M == Solver -> + + ( + What = {G} -> + + Res = [G|Gss] + ; + Res = [What|Gss] + ) + ; + Res = Gss + ), + nonlin_strip(Gs,Solver,Gss). + +all_attribute_goals([]) --> []. +all_attribute_goals([V|Vs]) --> + dump_linear(V,toplevel), + dump_nonzero(V,toplevel), + all_attribute_goals(Vs). + +% mapping(L1,L2,AssocIn,AssocOut) +% +% Makes an association mapping of lists L1 and L2: +% L1 = [L1H|L1T] and L2 = [L2H|L2T] then the association L1H-L2H is formed +% and the tails are mapped similarly. + +mapping([],[],D0,D0). +mapping([T|Ts],[N|Ns],D0,D2) :- + put_assoc(T,D0,N,D1), + mapping(Ts,Ns,D1,D2). + +% copy(Term,Copy,AssocIn,AssocOut) +% +% Makes a copy of Term by changing all variables in it to new ones and +% building an association between original variables and the new ones. +% E.g. when Term = test(A,B,C), Copy = test(D,E,F) and an association between +% A and D, B and E and C and F is formed in AssocOut. AssocIn is input association. + +copy(Term,Copy,D0,D1) :- + var(Term), + ( + get_assoc(Term,D0,New) -> + + Copy = New, + D1 = D0 + ; + put_assoc(Term,D0,Copy,D1) + ). +copy(Term,Copy,D0,D1) :- + nonvar(Term), % Term is a functor + functor(Term,N,A), + functor(Copy,N,A), % Copy is new functor with the same name and arity as Term + copy(A,Term,Copy,D0,D1). + +% copy(Nb,Term,Copy,AssocIn,AssocOut) +% +% Makes a copy of the Nb arguments of Term by changing all variables in it to +% new ones and building an association between original variables and the new ones. +% See also copy/4 + +copy(0,_,_,D0,D0) :- !. +copy(1,T,C,D0,D1) :- + !, + arg(1,T,At1), + arg(1,C,Ac1), + copy(At1,Ac1,D0,D1). +copy(2,T,C,D0,D2) :- !, + arg(1,T,At1), + arg(1,C,Ac1), + copy(At1,Ac1,D0,D1), + arg(2,T,At2), + arg(2,C,Ac2), + copy(At2,Ac2,D1,D2). +copy(N,T,C,D0,D2) :- + arg(N,T,At), + arg(N,C,Ac), + copy(At,Ac,D0,D1), + N1 is N-1, + copy(N1,T,C,D1,D2). + +end_of_file. diff --git a/LGPL/clpr/clpr/fourmotz.pl b/LGPL/clpr/clpr/fourmotz.pl new file mode 100644 index 000000000..6ca6c0aeb --- /dev/null +++ b/LGPL/clpr/clpr/fourmotz.pl @@ -0,0 +1,499 @@ +/* $Id: fourmotz.pl,v 1.1 2005-10-28 17:51:01 vsc Exp $ + + Part of CPL(R) (Constraint Logic Programming over Reals) + + Author: Leslie De Koninck + E-mail: Tom.Schrijvers@cs.kuleuven.ac.be + WWW: http://www.swi-prolog.org + http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 + Copyright (C): 2004, K.U. Leuven and + 1992-1995, Austrian Research Institute for + Artificial Intelligence (OFAI), + Vienna, Austria + + This software is part of Leslie De Koninck's master thesis, supervised + by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) + by Christian Holzbaur for SICStus Prolog and distributed under the + license details below with permission from all mentioned authors. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 + of the License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + As a special exception, if you link this library with other files, + compiled with a Free Software compiler, to produce an executable, this + library does not by itself cause the resulting executable to be covered + by the GNU General Public License. This exception does not however + invalidate any other reasons why the executable file might be covered by + the GNU General Public License. +*/ + + +:- module(fourmotz, + [ + fm_elim/3 + ]). +:- use_module(bv, + [ + allvars/2, + basis_add/2, + detach_bounds/1, + pivot/4, + var_with_def_intern/4 + ]). +:- use_module(class, + [ + class_allvars/2 + ]). +:- use_module(project, + [ + drop_dep/1, + drop_dep_one/1, + make_target_indep/2 + ]). +:- use_module(redund, + [ + redundancy_vars/1 + ]). +:- use_module(store, + [ + add_linear_11/3, + add_linear_f1/4, + indep/2, + nf_coeff_of/3, + normalize_scalar/2 + ]). + + + +fm_elim(Vs,Target,Pivots) :- + prefilter(Vs,Vsf), + fm_elim_int(Vsf,Target,Pivots). + +% prefilter(Vars,Res) +% +% filters out target variables and variables that do not occur in bounded linear equations. +% Stores that the variables in Res are to be kept independent. + +prefilter([],[]). +prefilter([V|Vs],Res) :- + ( + get_attr(V,itf3,(Ty,St,Li,Or,Cl,Fo,No,n,_,Ke)), + occurs(V) -> % V is a nontarget variable that occurs in a bounded linear equation + + Res = [V|Tail], + put_attr(V,itf3,(Ty,St,Li,Or,Cl,Fo,No,n,keep_indep,Ke)), + prefilter(Vs,Tail) + ; + prefilter(Vs,Res) + ). + +% +% the target variables are marked with an attribute, and we get a list +% of them as an argument too +% +fm_elim_int([],_,Pivots) :- % done + unkeep(Pivots). +fm_elim_int(Vs,Target,Pivots) :- + Vs = [_|_], + ( + best(Vs,Best,Rest) -> + + occurences(Best,Occ), + elim_min(Best,Occ,Target,Pivots,NewPivots) + ; % give up + NewPivots = Pivots, + Rest = [] + ), + fm_elim_int(Rest,Target,NewPivots). + +% best(Vs,Best,Rest) +% +% Finds the variable with the best result (lowest Delta) in fm_cp_filter +% and returns the other variables in Rest. + +best(Vs,Best,Rest) :- + findall(Delta-N,fm_cp_filter(Vs,Delta,N),Deltas), + keysort(Deltas,[_-N|_]), + select_nth(Vs,N,Best,Rest). + +% fm_cp_filter(Vs,Delta,N) +% +% For an indepenent variable V in Vs, which is the N'th element in Vs, +% find how many inequalities are generated when this variable is eliminated. +% Note that target variables and variables that only occur in unbounded equations +% should have been removed from Vs via prefilter/2 + +fm_cp_filter(Vs,Delta,N) :- + length(Vs,Len), % Len = number of variables in Vs + mem(Vs,X,Vst), % Selects a variable X in Vs, Vst is the list of elements after X in Vs + get_attr(X,itf3,(_,_,lin(Lin),order(OrdX),_,_,_,n,_)), % no target variable + indep(Lin,OrdX), % X is an independent variable + occurences(X,Occ), + Occ = [_|_], + cp_card(Occ,0,Lnew), + length(Occ,Locc), + Delta is Lnew-Locc, + length(Vst,Vstl), + N is Len-Vstl. % X is the Nth element in Vs + +% mem(Xs,X,XsT) +% +% If X is a member of Xs, XsT is the list of elements after X in Xs. + +mem([X|Xs],X,Xs). +mem([_|Ys],X,Xs) :- mem(Ys,X,Xs). + +% select_nth(List,N,Nth,Others) +% +% Selects the N th element of List, stores it in Nth and returns the rest of the list in Others. + +select_nth(List,N,Nth,Others) :- + select_nth(List,1,N,Nth,Others). + +select_nth([X|Xs],N,N,X,Xs) :- !. +select_nth([Y|Ys],M,N,X,[Y|Xs]) :- + M1 is M+1, + select_nth(Ys,M1,N,X,Xs). + +% +% fm_detach + reverse_pivot introduce indep t_none, which +% invalidates the invariants +% +elim_min(V,Occ,Target,Pivots,NewPivots) :- + crossproduct(Occ,New,[]), + activate_crossproduct(New), + reverse_pivot(Pivots), + fm_detach(Occ), + allvars(V,All), + redundancy_vars(All), % only for New \== [] + make_target_indep(Target,NewPivots), + drop_dep(All). + +% +% restore NF by reverse pivoting +% +reverse_pivot([]). +reverse_pivot([I:D|Ps]) :- + get_attr(D,itf3,(type(Dt),St,Li,Or,Cl,Fo,No,Ta,KI,_)), + put_attr(D,itf3,(type(Dt),St,Li,Or,Cl,Fo,No,Ta,KI,n)), % no longer + get_attr(I,itf3,(_,_,_,order(OrdI),class(ClI),_)), + pivot(D,ClI,OrdI,Dt), + reverse_pivot(Ps). + +% unkeep(Pivots) +% +% + +unkeep([]). +unkeep([_:D|Ps]) :- + get_attr(D,itf3,(Ty,St,Li,Or,Cl,Fo,No,Ta,KI,_)), + put_attr(D,itf3,(Ty,St,Li,Or,Cl,Fo,No,Ta,KI,n)), + drop_dep_one(D), + unkeep(Ps). + + +% +% All we drop are bounds +% +fm_detach( []). +fm_detach([V:_|Vs]) :- + detach_bounds(V), + fm_detach(Vs). + +% activate_crossproduct(Lst) +% +% For each inequality Lin =< 0 (or Lin < 0) in Lst, a new variable is created: +% Var = Lin and Var =< 0 (or Var < 0). Var is added to the basis. + +activate_crossproduct([]). +activate_crossproduct([lez(Strict,Lin)|News]) :- + Z = 0.0, + var_with_def_intern(t_u(Z),Var,Lin,Strict), + % Var belongs to same class as elements in Lin + basis_add(Var,_), + activate_crossproduct(News). + +% ------------------------------------------------------------------------------ + +% crossproduct(Lst,Res,ResTail) +% +% See crossproduct/4 +% This predicate each time puts the next element of Lst as First in crossproduct/4 +% and lets the rest be Next. + +crossproduct([]) --> []. +crossproduct([A|As]) --> + crossproduct(As,A), + crossproduct(As). + +% crossproduct(Next,First,Res,ResTail) +% +% Eliminates a variable in linear equations First + Next and stores the generated +% inequalities in Res. +% Let's say A:K1 = First and B:K2 = first equation in Next. +% A = ... + K1*V + ... +% B = ... + K2*V + ... +% Let K = -K2/K1 +% then K*A + B = ... + 0*V + ... +% from the bounds of A and B, via cross_lower/7 and cross_upper/7, new inequalities +% are generated. Then the same is done for B:K2 = next element in Next. + +crossproduct([],_) --> []. +crossproduct([B:Kb|Bs],A:Ka) --> + { + get_attr(A,itf3,(type(Ta),strictness(Sa),lin(LinA),_)), + get_attr(B,itf3,(type(Tb),strictness(Sb),lin(LinB),_)), + K is -Kb/Ka, + add_linear_f1(LinA,K,LinB,Lin) % Lin doesn't contain the target variable anymore + }, + ( + { 0.0 - K < -1e-010 } -> % K > 0: signs were opposite + + { Strict is Sa \/ Sb }, + cross_lower(Ta,Tb,K,Lin,Strict), + cross_upper(Ta,Tb,K,Lin,Strict) + ; % La =< A =< Ua -> -Ua =< -A =< -La + + { + flip(Ta,Taf), + flip_strict(Sa,Saf), + Strict is Saf \/ Sb + }, + cross_lower(Taf,Tb,K,Lin,Strict), + cross_upper(Taf,Tb,K,Lin,Strict) + ), + crossproduct(Bs,A:Ka). + +% cross_lower(Ta,Tb,K,Lin,Strict,Res,ResTail) +% +% Generates a constraint following from the bounds of A and B. +% When A = LinA and B = LinB then Lin = K*LinA + LinB. Ta is the type +% of A and Tb is the type of B. Strict is the union of the strictness +% of A and B. If K is negative, then Ta should have been flipped (flip/2). +% The idea is that if La =< A =< Ua and Lb =< B =< Ub (=< can also be <) +% then if K is positive, K*La + Lb =< K*A + B =< K*Ua + Ub. +% if K is negative, K*Ua + Lb =< K*A + B =< K*La + Ub. +% This predicate handles the first inequality and adds it to Res in the form +% lez(Sl,Lhs) meaning K*La + Lb - (K*A + B) =< 0 or K*Ua + Lb - (K*A + B) =< 0 +% with Sl being the strictness and Lhs the lefthandside of the equation. +% See also cross_upper/7 + +cross_lower(Ta,Tb,K,Lin,Strict) --> + { + lower(Ta,La), + lower(Tb,Lb), + !, + L is K*La+Lb, + normalize_scalar(L,Ln), + Mone = -1.0, + add_linear_f1(Lin,Mone,Ln,Lhs), + Sl is Strict >> 1 % normalize to upper bound + }, + [ lez(Sl,Lhs) ]. +cross_lower(_,_,_,_,_) --> []. + +% cross_upper(Ta,Tb,K,Lin,Strict,Res,ResTail) +% +% See cross_lower/7 +% This predicate handles the second inequality: +% -(K*Ua + Ub) + K*A + B =< 0 or -(K*La + Ub) + K*A + B =< 0 + +cross_upper(Ta,Tb,K,Lin,Strict) --> + { + upper(Ta,Ua), + upper(Tb,Ub), + !, + U is -(K*Ua+Ub), + normalize_scalar(U,Un), + add_linear_11(Un,Lin,Lhs), + Su is Strict /\ 1 % normalize to upper bound + }, + [ lez(Su,Lhs) ]. +cross_upper(_,_,_,_,_) --> []. + +% lower(Type,Lowerbound) +% +% Returns the lowerbound of type Type if it has one. +% E.g. if type = t_l(L) then Lowerbound is L, +% if type = t_lU(L,U) then Lowerbound is L, +% if type = t_u(U) then fails + +lower(t_l(L),L). +lower(t_lu(L,_),L). +lower(t_L(L),L). +lower(t_Lu(L,_),L). +lower(t_lU(L,_),L). + +% upper(Type,Upperbound) +% +% Returns the upperbound of type Type if it has one. +% See lower/2 + +upper(t_u(U),U). +upper(t_lu(_,U),U). +upper(t_U(U),U). +upper(t_Lu(_,U),U). +upper(t_lU(_,U),U). + +% flip(Type,FlippedType) +% +% Flips the lower and upperbound, so the old lowerbound becomes the new upperbound and +% vice versa. + +flip(t_l(X),t_u(X)). +flip(t_u(X),t_l(X)). +flip(t_lu(X,Y),t_lu(Y,X)). +flip(t_L(X),t_u(X)). +flip(t_U(X),t_l(X)). +flip(t_lU(X,Y),t_lu(Y,X)). +flip(t_Lu(X,Y),t_lu(Y,X)). + +% flip_strict(Strict,FlippedStrict) +% +% Does what flip/2 does, but for the strictness. + +flip_strict(0,0). +flip_strict(1,2). +flip_strict(2,1). +flip_strict(3,3). + +% cp_card(Lst,CountIn,CountOut) +% +% Counts the number of bounds that may generate an inequality in +% crossproduct/3 + +cp_card([],Ci,Ci). +cp_card([A|As],Ci,Co) :- + cp_card(As,A,Ci,Cii), + cp_card(As,Cii,Co). + +% cp_card(Next,First,CountIn,CountOut) +% +% Counts the number of bounds that may generate an inequality in +% crossproduct/4. + +cp_card([],_,Ci,Ci). +cp_card([B:Kb|Bs],A:Ka,Ci,Co) :- + get_attr(A,itf3,(type(Ta),_)), + get_attr(B,itf3,(type(Tb),_)), + K is -Kb/Ka, + ( + 0.0- K < -1e-010 -> % K > 0: signs were opposite + + cp_card_lower(Ta,Tb,Ci,Cii), + cp_card_upper(Ta,Tb,Cii,Ciii) + ; + flip(Ta,Taf), + cp_card_lower(Taf,Tb,Ci,Cii), + cp_card_upper(Taf,Tb,Cii,Ciii) + ), + cp_card(Bs,A:Ka,Ciii,Co). + +% cp_card_lower(TypeA,TypeB,SIn,SOut) +% +% SOut = SIn + 1 if both TypeA and TypeB have a lowerbound. + +cp_card_lower(Ta,Tb,Si,So) :- + lower(Ta,_), + lower(Tb,_), + !, + So is Si+1. +cp_card_lower(_,_,Si,Si). + +% cp_card_upper(TypeA,TypeB,SIn,SOut) +% +% SOut = SIn + 1 if both TypeA and TypeB have an upperbound. + +cp_card_upper(Ta,Tb,Si,So) :- + upper(Ta,_), + upper(Tb,_), + !, + So is Si+1. +cp_card_upper(_,_,Si,Si). + +% ------------------------------------------------------------------------------ + +% occurences(V,Occ) +% +% Returns in Occ the occurrences of variable V in the linear equations of dependent variables +% with bound =\= t_none in the form of D:K where D is a dependent variable and K is the scalar +% of V in the linear equation of D. + +occurences(V,Occ) :- + get_attr(V,itf3,(_,_,_,order(OrdV),class(C),_)), + class_allvars(C,All), + occurences(All,OrdV,Occ). + +% occurences(De,OrdV,Occ) +% +% Returns in Occ the occurrences of variable V with order OrdV in the linear equations of +% dependent variables De with bound =\= t_none in the form of D:K where D is a dependent +% variable and K is the scalar of V in the linear equation of D. + +occurences(De,_,[]) :- + var(De), + !. +occurences([D|De],OrdV,Occ) :- + ( + get_attr(D,itf3,(type(Type),_,lin(Lin),_)), + occ_type_filter(Type), + nf_coeff_of(Lin,OrdV,K) -> + + Occ = [D:K|Occt], + occurences(De,OrdV,Occt) + ; + occurences(De,OrdV,Occ) + ). + +% occ_type_filter(Type) +% +% Succeeds when Type is any other type than t_none. Is used in occurences/3 and occurs/2 + +occ_type_filter(t_l(_)). +occ_type_filter(t_u(_)). +occ_type_filter(t_lu(_,_)). +occ_type_filter(t_L(_)). +occ_type_filter(t_U(_)). +occ_type_filter(t_lU(_,_)). +occ_type_filter(t_Lu(_,_)). + +% occurs(V) +% +% Checks whether variable V occurs in a linear equation of a dependent variable with a bound +% =\= t_none. + +occurs(V) :- + get_attr(V,itf3,(_,_,_,order(OrdV),class(C),_)), + class_allvars(C,All), + occurs(All,OrdV). + +% occurs(De,OrdV) +% +% Checks whether variable V with order OrdV occurs in a linear equation of any dependent variable +% in De with a bound =\= t_none. + +occurs(De,_) :- + var(De), + !, + fail. +occurs([D|De],OrdV) :- + ( + get_attr(D,itf3,(type(Type),_,lin(Lin),_)), + occ_type_filter(Type), + nf_coeff_of(Lin,OrdV,_) -> + + true + ; + occurs(De,OrdV) + ). diff --git a/LGPL/clpr/clpr/geler.pl b/LGPL/clpr/clpr/geler.pl new file mode 100644 index 000000000..a86472bbf --- /dev/null +++ b/LGPL/clpr/clpr/geler.pl @@ -0,0 +1,230 @@ +/* $Id: geler.pl,v 1.1 2005-10-28 17:51:01 vsc Exp $ + + Part of CPL(R) (Constraint Logic Programming over Reals) + + Author: Leslie De Koninck + E-mail: Tom.Schrijvers@cs.kuleuven.ac.be + WWW: http://www.swi-prolog.org + http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 + Copyright (C): 2004, K.U. Leuven and + 1992-1995, Austrian Research Institute for + Artificial Intelligence (OFAI), + Vienna, Austria + + This software is part of Leslie De Koninck's master thesis, supervised + by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) + by Christian Holzbaur for SICStus Prolog and distributed under the + license details below with permission from all mentioned authors. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 + of the License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + As a special exception, if you link this library with other files, + compiled with a Free Software compiler, to produce an executable, this + library does not by itself cause the resulting executable to be covered + by the GNU General Public License. This exception does not however + invalidate any other reasons why the executable file might be covered by + the GNU General Public License. +*/ + + +:- module(geler_r, + [ + geler/2, + project_nonlin/3, + collect_nonlin/3 + ]). + +%:- attribute goals/1, all_nonlin/1. + +attribute_goal(X,Goals) :- + get_attr(X,geler_r,g(goals(Gs),_)), + nonexhausted(Gs,Goals,[]), + Goals = [_|_]. +attribute_goal(X,Conj) :- + get_attr(X,geler_r,g(_,all_nonlin(Goals))), + l2conj(Goals,Conj). + +% l2conj(List,Conj) +% +% turns a List into a conjunction of the form (El,Conj) where Conj +% is of the same form recursively and El is an element of the list + +l2conj([X|Xs],Conj) :- + ( + Xs = [], + Conj = X + ; + Xs = [_|_], + Conj = (X,Xc), + l2conj(Xs,Xc) + ). + +% nonexhausted(Goals,OutList,OutListTail) +% +% removes the goals that have already run from Goals +% and puts the result in the difference list OutList + +nonexhausted(run(Mutex,G)) --> + ( + {var(Mutex)} -> + + [G] + ; + [] + ). +nonexhausted((A,B)) --> + nonexhausted(A), + nonexhausted(B). + +attr_unify_hook(g(goals(Gx),_),Y) :- + !, + ( + var(Y), + ( + % possibly mutual goals. these need to be run. other goals are run + % as well to remove redundant goals. + get_attr(Y,geler_r,g(goals(Gy),Other)) -> + + Later = [Gx,Gy], + ( + Other = n -> + + del_attr(Y,geler_r) + ; + put_attr(Y,geler_r,g(n,Other)) + ) + ; + % no goals in Y, so no mutual goals of X and Y, store goals of X in Y + % no need to run any goal. + get_attr(Y,geler_r,g(n,Other)) -> + + Later = [], + put_attr(Y,geler_r,g(goals(Gx),Other)) + ; + Later = [], + put_attr(Y,geler_r,g(goals(Gx),n)) + ) + ; + nonvar(Y), + Later = [Gx] + ), + call_list(Later). +attr_unify_hook(_,_). % no goals in X + +% call_list(List) +% +% Calls all the goals in List. + +call_list([]). +call_list([G|Gs]) :- + call(G), + call_list(Gs). + + + + +% +% called from project.pl +% +project_nonlin(_,Cvas,Reachable) :- + collect_nonlin(Cvas,L,[]), + sort(L,Ls), + term_variables(Ls,Reachable). + %put_attr(_,all_nonlin(Ls)). + + +collect_nonlin([]) --> []. +collect_nonlin([X|Xs]) --> + ( + {get_attr(X,geler_r,g(goals(Gx),_))} -> + + trans(Gx), + collect_nonlin(Xs) + ; + collect_nonlin(Xs) + ). + +% trans(Goals,OutList,OutListTail) +% +% transforms the goals (of the form run(Mutex,Goal) +% that are in Goals (in the conjunction form, see also l2conj) +% that have not been run (var(Mutex) into a readable output format +% and notes that they're done (Mutex = done). Because of the Mutex +% variable, each goal is only added once (so not for each variable). + +trans((A,B)) --> + trans(A), + trans(B). +trans(run(Mutex,Gs)) --> + ( + {var(Mutex)} -> + {Mutex = done}, + transg(Gs) + ; + [] + ). + +transg((A,B)) --> + !, + transg(A), + transg(B). +transg(M:G) --> + !, + M:transg(G). +transg(G) --> [G]. + +% run(Mutex,G) +% +% Calls goal G if it has not yet run (Mutex is still variable) +% and stores that it has run (Mutex = done). This is done so +% that when X = Y and X and Y are in the same goal, that goal +% is called only once. + +run(Mutex,_) :- nonvar(Mutex). +run(Mutex,G) :- + var(Mutex), + Mutex = done, + call(G). + +% geler(Vars,Goal) +% +% called by nf.pl when an unsolvable non-linear expression is found +% Vars contain the variables of the expression, Goal contains the predicate of nf.pl to be called when +% the variables are bound. + +geler(Vars,Goal) :- + attach(Vars,run(_Mutex,Goal)). + % one goal gets the same mutex on every var, so it is run only once + +% attach(Vars,Goal) +% +% attaches a new goal to be awoken when the variables get bounded. +% when the old value of the attribute goals = OldGoal, then the new value = (Goal,OldGoal) + +attach([],_). +attach([V|Vs],Goal) :- + ( + var(V), + get_attr(V,geler_r,g(goals(Gv),Other)) -> + + put_attr(V,geler_r,g(goals((Goal,Gv)),Other)) + ; + get_attr(V,geler_r,(n,Other)) -> + + put_attr(V,geler_r,g(goals(Goal),Other)) + ; + put_attr(V,geler_r,g(goals(Goal),n)) + ), + attach(Vs,Goal). diff --git a/LGPL/clpr/clpr/ineq.pl b/LGPL/clpr/clpr/ineq.pl new file mode 100644 index 000000000..32926cc5a --- /dev/null +++ b/LGPL/clpr/clpr/ineq.pl @@ -0,0 +1,1530 @@ +/* $Id: ineq.pl,v 1.1 2005-10-28 17:51:01 vsc Exp $ + + Part of CPL(R) (Constraint Logic Programming over Reals) + + Author: Leslie De Koninck + E-mail: Tom.Schrijvers@cs.kuleuven.ac.be + WWW: http://www.swi-prolog.org + http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 + Copyright (C): 2004, K.U. Leuven and + 1992-1995, Austrian Research Institute for + Artificial Intelligence (OFAI), + Vienna, Austria + + This software is part of Leslie De Koninck's master thesis, supervised + by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) + by Christian Holzbaur for SICStus Prolog and distributed under the + license details below with permission from all mentioned authors. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 + of the License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + As a special exception, if you link this library with other files, + compiled with a Free Software compiler, to produce an executable, this + library does not by itself cause the resulting executable to be covered + by the GNU General Public License. This exception does not however + invalidate any other reasons why the executable file might be covered by + the GNU General Public License. +*/ + + +:- module(ineq, + [ + ineq/4, + ineq_one/4, + ineq_one_n_n_0/1, + ineq_one_n_p_0/1, + ineq_one_s_n_0/1, + ineq_one_s_p_0/1 + ]). +:- use_module(bv, + [ + backsubst/3, + backsubst_delta/4, + basis_add/2, + dec_step/2, + deref/2, + determine_active_dec/1, + determine_active_inc/1, + export_binding/1, + get_or_add_class/2, + inc_step/2, + lb/3, + pivot_a/4, + rcbl_status/6, + reconsider/1, + same_class/2, + solve/1, + ub/3, + unconstrained/4, + var_intern/3, + var_with_def_intern/4 + ]). +:- use_module(store, + [ + add_linear_11/3, + add_linear_ff/5, + normalize_scalar/2 + ]). + +% ineq(H,I,Nf,Strictness) +% +% Solves the inequality Nf < 0 or Nf =< 0 where Nf is in normal form +% and H and I are the homogene and inhomogene parts of Nf. + +ineq([],I,_,Strictness) :- ineq_ground(Strictness,I). +ineq([v(K,[X^1])|Tail],I,Lin,Strictness) :- + ineq_cases(Tail,I,Lin,Strictness,X,K). + +ineq_cases([],I,_,Strictness,X,K) :- % K*X + I < 0 or K*X + I =< 0 + ineq_one(Strictness,X,K,I). +ineq_cases([_|_],_,Lin,Strictness,_,_) :- + deref(Lin,Lind), % Id+Hd =< 0 + Lind = [Inhom,_|Hom], + ineq_more(Hom,Inhom,Lind,Strictness). + +% ineq_ground(Strictness,I) +% +% Checks whether a grounded inequality I < 0 or I =< 0 is satisfied. + +ineq_ground(strict,I) :- I - 0.0 < -1e-010. % I < 0 +ineq_ground(nonstrict,I) :- I - 0.0 < 1e-010. % I =< 0 + +% ineq_one(Strictness,X,K,I) +% +% Solves the inequality K*X + I < 0 or K*X + I =< 9 + +ineq_one(strict,X,K,I) :- + ( + K > 1e-010 -> % K > 0 + + ( + I >= -1e-010, % I =:= 0 + I =< 1e-010 -> + + ineq_one_s_p_0(X) % K*X < 0, K > 0 => X < 0 + ; + Inhom is I/K, + ineq_one_s_p_i(X,Inhom) % K*X + I < 0, K > 0 => X + I/K < 0 + ) + ; + ( + I >= -1e-010, % I =:= 0 + I =< 1e-010 -> + + ineq_one_s_n_0(X) % K*X < 0, K < 0 => -X < 0 + ; + Inhom is -I/K, + ineq_one_s_n_i(X,Inhom) % K*X + I < 0, K < 0 => -X - I/K < 0 + ) + ). +ineq_one(nonstrict,X,K,I) :- + ( + K > 1e-010 -> % K > 0 + + ( + I >= -1e-010, % I =:= 0 + I =< 1e-010 -> + + ineq_one_n_p_0(X) % K*X =< 0, K > 0 => X =< 0 + ; + Inhom is I/K, + ineq_one_n_p_i(X,Inhom) % K*X + I =< 0, K > 0 => X + I/K =< 0 + ) + ; + ( + I >= -1e-010, % I =:= 0 + I =< 1e-010 -> + + ineq_one_n_n_0(X) % K*X =< 0, K < 0 => -X =< 0 + ; + Inhom is -I/K, + ineq_one_n_n_i(X,Inhom) % K*X + I =< 0, K < 0 => -X - I/K =< 0 + ) + ). + +% --------------------------- strict ---------------------------- + +% ineq_one_s_p_0(X) +% +% Solves the inequality X < 0 + +ineq_one_s_p_0(X) :- + get_attr(X,itf3,(_,_,lin(LinX),_)), + !, % old variable, this is deref + LinX = [Ix,_|OrdX], + ineq_one_old_s_p_0(OrdX,X,Ix). +ineq_one_s_p_0(X) :- % new variable, nothing depends on it + var_intern(t_u(0.0),X,1). % put a strict inactive upperbound on the variable + +% ineq_one_s_n_0(X) +% +% Solves the inequality X > 0 + +ineq_one_s_n_0(X) :- + get_attr(X,itf3,(_,_,lin(LinX),_)), + !, + LinX = [Ix,_|OrdX], + ineq_one_old_s_n_0(OrdX,X,Ix). +ineq_one_s_n_0(X) :- + var_intern(t_l(0.0),X,2). % puts a strict inactive lowerbound on the variable + +% ineq_one_s_p_i(X,I) +% +% Solves the inequality X < -I + +ineq_one_s_p_i(X,I) :- + get_attr(X,itf3,(_,_,lin(LinX),_)), + !, + LinX = [Ix,_|OrdX], + ineq_one_old_s_p_i(OrdX,I,X,Ix). +ineq_one_s_p_i(X,I) :- + Bound is -I, + var_intern(t_u(Bound),X,1). % puts a strict inactive upperbound on the variable + +% ineq_one_s_n_i(X,I) +% +% Solves the inequality X > I + +ineq_one_s_n_i(X,I) :- + get_attr(X,itf3,(_,_,lin(LinX),_)), + !, + LinX = [Ix,_|OrdX], + ineq_one_old_s_n_i(OrdX,I,X,Ix). +ineq_one_s_n_i(X,I) :- var_intern(t_l(I),X,2). % puts a strict inactive lowerbound on the variable + +% ineq_one_old_s_p_0(Hom,X,Inhom) +% +% Solves the inequality X < 0 where X has linear equation Hom + Inhom + +ineq_one_old_s_p_0([],_,Ix) :- Ix < -1e-010. % X = I: Ix < 0 +ineq_one_old_s_p_0([l(Y*Ky,_)|Tail],X,Ix) :- + ( + Tail = [] -> % X = K*Y + I + + Bound is -Ix/Ky, + update_indep(strict,Y,Ky,Bound) % X < 0, X = K*Y + I => Y < -I/K or Y > -I/K (depending on K) + ; + Tail = [_|_] -> + + get_attr(X,itf3,(type(Type),strictness(Old),lin(Lin),_)), + udus(Type,X,Lin,0.0,Old) % update strict upperbound + ). + +% ineq_one_old_s_p_0(Hom,X,Inhom) +% +% Solves the inequality X > 0 where X has linear equation Hom + Inhom + +ineq_one_old_s_n_0([],_,Ix) :- Ix > 1e-010. % X = I: Ix > 0 +ineq_one_old_s_n_0([l(Y*Ky,_)|Tail], X, Ix) :- + ( + Tail = [] -> % X = K*Y + I + + Coeff is -Ky, + Bound is Ix/Coeff, + update_indep(strict,Y,Coeff,Bound) + ; + Tail = [_|_] -> + + get_attr(X,itf3,(type(Type),strictness(Old),lin(Lin),_)), + udls(Type,X,Lin,0.0,Old) % update strict lowerbound + ). + +% ineq_one_old_s_p_i(Hom,C,X,Inhom) +% +% Solves the inequality X + C < 0 where X has linear equation Hom + Inhom + +ineq_one_old_s_p_i([],I,_,Ix) :- Ix + I < -1e-010. % X = I +ineq_one_old_s_p_i([l(Y*Ky,_)|Tail],I,X,Ix) :- + ( + Tail = [] -> % X = K*Y + I + + Bound is -(Ix+I)/Ky, + update_indep(strict,Y,Ky,Bound) + ; + Tail = [_|_] -> + + Bound is -I, + get_attr(X,itf3,(type(Type),strictness(Old),lin(Lin),_)), + udus(Type,X,Lin,Bound,Old) % update strict upperbound + ). + +% ineq_one_old_s_n_i(Hom,C,X,Inhom) +% +% Solves the inequality X - C > 0 where X has linear equation Hom + Inhom + +ineq_one_old_s_n_i([],I,_,Ix) :- -Ix + I < -1e-010. % X = I +ineq_one_old_s_n_i([l(Y*Ky,_)|Tail],I,X,Ix) :- + ( + Tail = [] -> % X = K*Y + I + + Coeff is -Ky, + Bound is (Ix-I)/Coeff, + update_indep(strict,Y,Coeff,Bound) + ; + Tail = [_|_] -> + + get_attr(X,itf3,(type(Type),strictness(Old),lin(Lin),_)), + udls(Type,X,Lin,I,Old) % update strict lowerbound + ). + +% -------------------------- nonstrict -------------------------- + +% ineq_one_n_p_0(X) +% +% Solves the inequality X =< 0 + +ineq_one_n_p_0(X) :- + get_attr(X,itf3,(_,_,lin(LinX),_)), + !, % old variable, this is deref + LinX = [Ix,_|OrdX], + ineq_one_old_n_p_0(OrdX,X,Ix). +ineq_one_n_p_0(X) :- % new variable, nothing depends on it + var_intern(t_u(0.0),X,0). % nonstrict upperbound + +% ineq_one_n_n_0(X) +% +% Solves the inequality X >= 0 + +ineq_one_n_n_0(X) :- + get_attr(X,itf3,(_,_,lin(LinX),_)), + !, + LinX = [Ix,_|OrdX], + ineq_one_old_n_n_0(OrdX,X,Ix). +ineq_one_n_n_0(X) :- + var_intern(t_l(0.0),X,0). % nonstrict lowerbound + +% ineq_one_n_p_i(X,I) +% +% Solves the inequality X =< -I + +ineq_one_n_p_i(X,I) :- + get_attr(X,itf3,(_,_,lin(LinX),_)), + !, + LinX = [Ix,_|OrdX], + ineq_one_old_n_p_i(OrdX,I,X,Ix). +ineq_one_n_p_i(X,I) :- + Bound is -I, + var_intern(t_u(Bound),X,0). % nonstrict upperbound + +% ineq_one_n_n_i(X,I) +% +% Solves the inequality X >= I + +ineq_one_n_n_i(X,I) :- + get_attr(X,itf3,(_,_,lin(LinX),_)), + !, + LinX = [Ix,_|OrdX], + ineq_one_old_n_n_i(OrdX,I,X,Ix). +ineq_one_n_n_i(X,I) :- + var_intern(t_l(I),X,0). % nonstrict lowerbound + +% ineq_one_old_n_p_0(Hom,X,Inhom) +% +% Solves the inequality X =< 0 where X has linear equation Hom + Inhom + +ineq_one_old_n_p_0([],_,Ix) :- Ix < 1e-010. % X =I +ineq_one_old_n_p_0([l(Y*Ky,_)|Tail],X,Ix) :- + ( + Tail = [] -> % X = K*Y + I + + Bound is -Ix/Ky, + update_indep(nonstrict,Y,Ky,Bound) + ; + Tail = [_|_] -> + + get_attr(X,itf3,(type(Type),strictness(Old),lin(Lin),_)), + udu(Type,X,Lin,0.0,Old) % update nonstrict upperbound + ). + +% ineq_one_old_n_n_0(Hom,X,Inhom) +% +% Solves the inequality X >= 0 where X has linear equation Hom + Inhom + +ineq_one_old_n_n_0([],_,Ix) :- Ix > -1e-010. % X = I +ineq_one_old_n_n_0([l(Y*Ky,_)|Tail], X, Ix) :- + ( + Tail = [] -> % X = K*Y + I + + Coeff is -Ky, + Bound is Ix/Coeff, + update_indep(nonstrict,Y,Coeff,Bound) + ; + Tail = [_|_] -> + + get_attr(X,itf3,(type(Type),strictness(Old),lin(Lin),_)), + udl(Type,X,Lin,0.0,Old) % update nonstrict lowerbound + ). + +% ineq_one_old_n_p_i(Hom,C,X,Inhom) +% +% Solves the inequality X + C =< 0 where X has linear equation Hom + Inhom + +ineq_one_old_n_p_i([],I,_,Ix) :- Ix + I < 1e-010. % X = I +ineq_one_old_n_p_i([l(Y*Ky,_)|Tail],I,X,Ix) :- + ( + Tail = [] -> % X = K*Y + I + + Bound is -(Ix+I)/Ky, + update_indep(nonstrict,Y,Ky,Bound) + ; + Tail = [_|_] -> + + Bound is -I, + get_attr(X,itf3,(type(Type),strictness(Old),lin(Lin),_)), + udu(Type,X,Lin,Bound,Old) % update nonstrict upperbound + ). + +% ineq_one_old_n_n_i(Hom,C,X,Inhom) +% +% Solves the inequality X - C >= 0 where X has linear equation Hom + Inhom + +ineq_one_old_n_n_i([],I,_,Ix) :- -Ix + I < 1e-010. % X = I +ineq_one_old_n_n_i([l(Y*Ky,_)|Tail],I,X,Ix) :- + ( + Tail = [] -> + + Coeff is -Ky, + Bound is (Ix-I)/Coeff, + update_indep(nonstrict,Y,Coeff,Bound) + ; + Tail = [_|_] -> + + get_attr(X,itf3,(type(Type),strictness(Old),lin(Lin),_)), + udl(Type,X,Lin,I,Old) + ). + +% --------------------------------------------------------------- + +% ineq_more(Hom,Inhom,Lin,Strictness) +% +% Solves the inequality Lin < 0 or Lin =< 0 with Lin = Hom + Inhom + +ineq_more([],I,_,Strictness) :- ineq_ground(Strictness,I). % I < 0 or I =< 0 +ineq_more([l(X*K,_)|Tail],Id,Lind,Strictness) :- + ( % X*K < Id or X*K =< Id + Tail = [] -> % one var: update bound instead of slack introduction + + get_or_add_class(X,_), % makes sure X belongs to a class + Bound is -Id/K, + update_indep(Strictness,X,K,Bound) % new bound + ; + Tail = [_|_] -> + + ineq_more(Strictness,Lind) + ). + +% ineq_more(Strictness,Lin) +% +% Solves the inequality Lin < 0 or Lin =< 0 + +ineq_more(strict,Lind) :- + ( + unconstrained(Lind,U,K,Rest) -> % never fails, no implied value + % Lind < 0 => Rest < -K*U where U has no bounds + + var_intern(t_l(0.0),S,2), % create slack variable S + get_attr(S,itf3,(_,_,_,order(OrdS),_)), + Ki is -1.0/K, + add_linear_ff(Rest,Ki,[0.0,0.0,l(S*1.0,OrdS)],Ki,LinU), % U = (-1/K)*Rest + (-1/K)*S + LinU = [_,_|Hu], + get_or_add_class(U,Class), + same_class(Hu,Class), % put all variables of new lin. eq. of U in the same class + get_attr(U,itf3,(_,_,_,order(OrdU),class(ClassU),_)), + backsubst(ClassU,OrdU,LinU) % substitute U by new lin. eq. everywhere in the class + ; + var_with_def_intern(t_u(0.0),S,Lind,1), % Lind < 0 => Lind = S with S < 0 + basis_add(S,_), % adds S to the basis + determine_active_dec(Lind), % activate bounds + reconsider(S) % reconsider basis + ). +ineq_more(nonstrict,Lind) :- + ( + unconstrained(Lind,U,K,Rest) -> % never fails, no implied value + % Lind =< 0 => Rest =< -K*U where U has no bounds + var_intern(t_l(0.0),S,0), % create slack variable S + Ki is -1.0/K, + get_attr(S,itf3,(_,_,_,order(OrdS),_)), + add_linear_ff(Rest,Ki,[0.0,0.0,l(S*1.0,OrdS)],Ki,LinU), % U = (-1K)*Rest + (-1/K)*S + LinU = [_,_|Hu], + get_or_add_class(U,Class), + same_class(Hu,Class), % put all variables of new lin. eq of U in the same class + get_attr(U,itf3,(_,_,_,order(OrdU),class(ClassU),_)), + backsubst(ClassU,OrdU,LinU) % substitute U by new lin. eq. everywhere in the class + ; % all variables are constrained + var_with_def_intern(t_u(0.0),S,Lind,0), % Lind =< 0 => Lind = S with S =< 0 + basis_add(S,_), % adds S to the basis + determine_active_dec(Lind), + reconsider(S) + ). + + +% update_indep(Strictness,X,K,Bound) +% +% Updates the bound of independent variable X where X < Bound or X =< Bound +% or X > Bound or X >= Bound, depending on Strictness and K. + +update_indep(strict,X,K,Bound) :- + get_attr(X,itf3,(type(Type),strictness(Old),lin(Lin),_)), + ( + K < -1e-010 -> + + uils(Type,X,Lin,Bound,Old) % update independent lowerbound strict + ; + uius(Type,X,Lin,Bound,Old) % update independent upperbound strict + ). +update_indep(nonstrict,X,K,Bound) :- + get_attr(X,itf3,(type(Type),strictness(Old),lin(Lin),_)), + ( + K < -1e-010 -> + + uil(Type,X,Lin,Bound,Old) % update independent lowerbound nonstrict + ; + uiu(Type,X,Lin,Bound,Old) % update independent upperbound nonstrict + ). + + +% --------------------------------------------------------------------------------------- + +% +% Update a bound on a var xi +% +% a) independent variable +% +% a1) update inactive bound: done +% +% a2) update active bound: +% Determine [lu]b including most constraining row R +% If we are within: done +% else pivot(R,xi) and introduce bound via (b) +% +% a3) introduce a bound on an unconstrained var: +% All vars that depend on xi are unconstrained (invariant) -> +% the bound cannot invalidate any Lhs +% +% b) dependent variable +% +% repair upper or lower (maybe just swap with an unconstrained var from Rhs) +% + +% +% Sign = 1,0,-1 means inside,at,outside +% + +% Read following predicates as update (dependent/independent) (lowerbound/upperbound) (strict) + +% udl(Type,X,Lin,Bound,Strict) +% +% Updates lower bound of dependent variable X with linear equation +% Lin that had type Type and strictness Strict, to the new non-strict +% bound Bound. + +udl(t_none,X,Lin,Bound,_Sold) :- + get_attr(X,itf3,(_,_,Li,order(Ord),RAtt)), + put_attr(X,itf3,(type(t_l(Bound)),strictness(0),Li,order(Ord),RAtt)), + ( + unconstrained(Lin,Uc,Kuc,Rest) -> + + % X = Lin => -1/K*Rest + 1/K*X = U where U has no bounds + Ki is -1.0/Kuc, + add_linear_ff(Rest,Ki,[0.0,0.0,l(X* -1.0,Ord)],Ki,LinU), + get_attr(Uc,itf3,(_,_,_,order(OrdU),class(Class),_)), + backsubst(Class,OrdU,LinU) + ; + % no unconstrained variables in Lin: make X part of basis and reconsider + basis_add(X,_), + determine_active_inc(Lin), + reconsider(X) + ). +udl(t_l(L),X,Lin,Bound,Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + true % new bound is smaller than old one: keep old + ; + TestBL > 1e-010 -> + + % new bound is larger than old one: use new and reconsider basis + Strict is Sold /\ 1, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_l(Bound)),strictness(Strict),RAtt)), + reconsider_lower(X,Lin,Bound) % makes sure that Lin still satisfies lowerbound Bound + ; + true % new bound is equal to old one, new one is nonstrict: keep old + ). + +udl(t_u(U),X,Lin,Bound,_Sold) :- + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + fail % new bound is larger than upperbound: fail + ; + TestUB > 1e-010 -> + + % new bound is smaller than upperbound: add new and reconsider basis + get_attr(X,itf3,(_,RAtt)), + put_attr(X,itf3,(type(t_lu(Bound,U)),RAtt)), + reconsider_lower(X,Lin,Bound) % makes sure that Lin still satisfies lowerbound Bound + ; + solve_bound(Lin,Bound) % new bound is equal to upperbound: solve + ). +udl(t_lu(L,U),X,Lin,Bound,Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + true % smaller than lowerbound: keep + ; + TestBL > 1e-010 -> % larger than lowerbound: check upperbound + + TestUB is U - Bound, + ( + TestUB < -1e-010 -> + + fail % larger than upperbound: fail + ; + TestUB > 1e-010 -> + + % smaller than upperbound: use new and reconsider basis + Strict is Sold /\ 1, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lu(Bound,U)),strictness(Strict),RAtt)), + reconsider_lower(X,Lin,Bound) + ; + % equal to upperbound: if strictness matches => solve + Sold /\ 1 =:= 0, + solve_bound(Lin,Bound) + ) + ; + true % equal to lowerbound and nonstrict: keep + ). + +% udls(Type,X,Lin,Bound,Strict) +% +% Updates lower bound of dependent variable X with linear equation +% Lin that had type Type and strictness Strict, to the new strict +% bound Bound. + +udls(t_none,X,Lin,Bound,_Sold) :- + get_attr(X,itf3,(_,_,Li,order(Ord),RAtt)), + put_attr(X,itf3,(type(t_l(Bound)),strictness(2),Li,order(Ord),RAtt)), + ( + unconstrained(Lin,Uc,Kuc,Rest) -> + + % X = Lin => U = -1/K*Rest + 1/K*X with U an unconstrained variable + Ki is -1.0/Kuc, + add_linear_ff(Rest,Ki,[0.0,0.0,l(X* -1.0,Ord)],Ki,LinU), + get_attr(Uc,itf3,(_,_,_,order(OrdU),class(Class),_)), + backsubst(Class,OrdU,LinU) + ; + % no unconstrained variables: add X to basis and reconsider basis + basis_add(X,_), + determine_active_inc(Lin), + reconsider(X) + ). +udls(t_l(L),X,Lin,Bound,Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + true % smaller than lowerbound: keep + ; + TestBL > 1e-010 -> + + % larger than lowerbound: use new and reconsider basis + Strict is Sold \/ 2, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_l(Bound)),strictness(Strict),RAtt)), + reconsider_lower(X,Lin,Bound) + ; + % equal to lowerbound: check strictness + Strict is Sold \/ 2, + get_attr(X,itf3,(Ty,_,RAtt)), + put_attr(X,itf3,(Ty,strictness(Strict),RAtt)) + ). +udls(t_u(U),X,Lin,Bound,Sold) :- + U - Bound > 1e-010, % smaller than upperbound: set new bound + Strict is Sold \/ 2, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lu(Bound,U)),strictness(Strict),RAtt)), + reconsider_lower(X,Lin,Bound). +udls(t_lu(L,U),X,Lin,Bound,Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + true % smaller than lowerbound: keep + ; + TestBL > 1e-010 -> + + % larger than lowerbound: check upperbound and possibly use new and reconsider basis + U - Bound > 1e-010, + Strict is Sold \/ 2, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lu(Bound,U)),strictness(Strict),RAtt)), + reconsider_lower(X,Lin,Bound) + ; + % equal to lowerbound: put new strictness + Strict is Sold \/ 2, + get_attr(X,itf3,(Ty,_,RAtt)), + put_attr(X,itf3,(Ty,strictness(Strict),RAtt)) + ). + +% udu(Type,X,Lin,Bound,Strict) +% +% Updates upper bound of dependent variable X with linear equation +% Lin that had type Type and strictness Strict, to the new non-strict +% bound Bound. + +udu(t_none,X,Lin,Bound,_Sold) :- + get_attr(X,itf3,(_,_,Li,order(Ord),RAtt)), + put_attr(X,itf3,(type(t_u(Bound)),strictness(0),Li,order(Ord),RAtt)), + ( + unconstrained(Lin,Uc,Kuc,Rest) -> + + % X = Lin => U = -1/K*Rest + 1/K*X with U an unconstrained variable + Ki is -1.0/Kuc, + add_linear_ff(Rest,Ki,[0.0,0.0,l(X* -1.0,Ord)],Ki,LinU), + get_attr(Uc,itf3,(_,_,_,order(OrdU),class(Class),_)), + backsubst(Class,OrdU,LinU) + ; + % no unconstrained variables: add X to basis and reconsider basis + basis_add(X,_), + determine_active_dec(Lin), % try to lower R + reconsider(X) + ). +udu(t_u(U),X,Lin,Bound,Sold) :- + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + true % larger than upperbound: keep + ; + TestUB > 1e-010 -> + + % smaller than upperbound: update and reconsider basis + Strict is Sold /\ 2, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_u(Bound)),strictness(Strict),RAtt)), + reconsider_upper(X,Lin,Bound) + ; + true % equal to upperbound and nonstrict: keep + ). +udu(t_l(L),X,Lin,Bound,_Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + fail % smaller than lowerbound: fail + ; + TestBL > 1e-010 -> + + % larger than lowerbound: use new and reconsider basis + get_attr(X,itf3,(_,RAtt)), + put_attr(X,itf3,(type(t_lu(L,Bound)),RAtt)), + reconsider_upper(X,Lin,Bound) + ; + solve_bound(Lin,Bound) % equal to lowerbound: solve + ). +udu(t_lu(L,U),X,Lin,Bound,Sold) :- + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + true % larger than upperbound: keep + ; + TestUB > 1e-010 -> + + % smaller than upperbound: check lowerbound + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + fail % smaller than lowerbound: fail + ; + TestBL > 1e-010 -> + + % larger than lowerbound: update and reconsider basis + Strict is Sold /\ 2, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lu(L,Bound)),strictness(Strict),RAtt)), + reconsider_upper(X,Lin,Bound) + ; + % equal to lowerbound: check strictness and possibly solve + Sold /\ 2 =:= 0, + solve_bound(Lin,Bound) + ) + ; + true % equal to upperbound and nonstrict: keep + ). + +% udus(Type,X,Lin,Bound,Strict) +% +% Updates upper bound of dependent variable X with linear equation +% Lin that had type Type and strictness Strict, to the new strict +% bound Bound. + +udus(t_none,X,Lin,Bound,_Sold) :- + get_attr(X,itf3,(_,_,Li,order(Ord),RAtt)), + put_attr(X,itf3,(type(t_u(Bound)),strictness(1),Li,order(Ord),RAtt)), + ( + unconstrained(Lin,Uc,Kuc,Rest) -> + + % X = Lin => U = -1/K*Rest + 1/K*X with U an unconstrained variable + Ki is -1.0/Kuc, + add_linear_ff(Rest,Ki,[0.0,0.0,l(X* -1.0,Ord)],Ki,LinU), + get_attr(Uc,itf3,(_,_,_,order(OrdU),class(Class),_)), + backsubst(Class,OrdU,LinU) + ; + % no unconstrained variables: add X to basis and reconsider basis + basis_add(X,_), + determine_active_dec(Lin), + reconsider(X) + ). +udus(t_u(U),X,Lin,Bound,Sold) :- + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + true % larger than upperbound: keep + ; + TestUB > 1e-010 -> + + % smaller than upperbound: update bound and reconsider basis + Strict is Sold \/ 1, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_u(Bound)),strictness(Strict),RAtt)), + reconsider_upper(X,Lin,Bound) + ; + % equal to upperbound: set new strictness + Strict is Sold \/ 1, + get_attr(X,itf3,(Ty,_,RAtt)), + put_attr(X,itf3,(Ty,strictness(Strict),RAtt)) + ). +udus(t_l(L),X,Lin,Bound,Sold) :- + Bound - L > 1e-010, % larger than lowerbound: update and reconsider basis + Strict is Sold \/ 1, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lu(L,Bound)),strictness(Strict),RAtt)), + reconsider_upper(X,Lin,Bound). +udus(t_lu(L,U),X,Lin,Bound,Sold) :- + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + true % larger than upperbound: keep + ; + TestUB > 1e-010 -> + + % smaller than upperbound: check lowerbound, possibly update and reconsider basis + Bound - L > 1e-010, + Strict is Sold \/ 1, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lu(L,Bound)),strictness(Strict),RAtt)), + reconsider_upper(X,Lin,Bound) + ; + % equal to upperbound: update strictness + Strict is Sold \/ 1, + get_attr(X,itf3,(Ty,_,RAtt)), + put_attr(X,itf3,(Ty,strictness(Strict),RAtt)) + ). + +% uiu(Type,X,Lin,Bound,Strict) +% +% Updates upper bound of independent variable X with linear equation +% Lin that had type Type and strictness Strict, to the new non-strict +% bound Bound. + +uiu(t_none,X,_Lin,Bound,_) :- % X had no bounds + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_u(Bound)),strictness(0),RAtt)). +uiu(t_u(U),X,_Lin,Bound,Sold) :- + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + true % larger than upperbound: keep + ; + TestUB > 1e-010 -> + + % smaller than upperbound: update. + Strict is Sold /\ 2, % update strictness: strictness of lowerbound is kept, + % strictness of upperbound is set to non-strict + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_u(Bound)),strictness(Strict),RAtt)) + ; + true % equal to upperbound and nonstrict: keep + ). +uiu(t_l(L),X,Lin,Bound,_Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + fail % Lowerbound was smaller than new upperbound: fail + ; + TestBL > 1e-010 -> + + % Upperbound is larger than lowerbound: store new bound + get_attr(X,itf3,(_,RAtt)), + put_attr(X,itf3,(type(t_lu(L,Bound)),RAtt)) + ; + solve_bound(Lin,Bound) % Lowerbound was equal to new upperbound: solve + ). +uiu(t_L(L),X,Lin,Bound,_Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + fail % Same as for t_l + ; + TestBL > 1e-010 -> + + % Same as for t_l (new bound becomes t_Lu) + get_attr(X,itf3,(_,RAtt)), + put_attr(X,itf3,(type(t_Lu(L,Bound)),RAtt)) + ; + solve_bound(Lin,Bound) % Same as for t_l + ). +uiu(t_lu(L,U),X,Lin,Bound,Sold) :- + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + true % Upperbound was smaller than new bound: keep + ; + TestUB > 1e-010 -> + + TestBL is Bound-L, % Upperbound was larger than new bound: check lowerbound + ( + TestBL < -1e-010 -> + + fail % Lowerbound was larger than new bound: fail + ; + TestBL > 1e-010 -> + + % Lowerbound was smaller than new bound: store new bound + Strict is Sold /\ 2, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lu(L,Bound)),strictness(Strict),RAtt)) + ; + % Lowerbound was equal to new bound: solve + Sold /\ 2 =:= 0, % Only solve when strictness matches + solve_bound(Lin,Bound) + ) + ; + true % Upperbound was equal to new bound and new bound non-strict: keep + ). +uiu(t_Lu(L,U),X,Lin,Bound,Sold) :- % See t_lu case + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + true + ; + TestUB > 1e-010 -> + + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + fail + ; + TestBL > 1e-010 -> + + Strict is Sold /\ 2, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_Lu(L,Bound)),strictness(Strict),RAtt)) + ; + Sold /\ 2 =:= 0, + solve_bound(Lin,Bound) + ) + ; + true + ). +uiu(t_U(U),X,_Lin,Bound,Sold) :- + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + true % larger than upperbound: keep + ; + TestUB > 1e-010 -> + + % smaller than active upperbound: check how much active upperbound can be lowered. + % if enough, just lower bound, otherwise update the bound, make X dependent and reconsider basis + Strict is Sold /\ 2, + ( + get_attr(X,itf3,(_,_,_,order(OrdX),class(ClassX),_)), + lb(ClassX,OrdX,Vlb-Vb-Lb), + Bound - (Lb+U) < 1e-010 -> + + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_U(Bound)),strictness(Strict),RAtt)), + pivot_a(Vlb,X,Vb,t_u(Bound)), + reconsider(X) + ; + get_attr(X,itf3,(_,_,Li,order(OrdX),class(ClassX),RAtt)), + put_attr(X,itf3,(type(t_U(Bound)),strictness(Strict),Li,order(OrdX),class(ClassX),RAtt)), + Delta is Bound-U, + backsubst_delta(ClassX,OrdX,X,Delta) + ) + ; + true % equal to upperbound and non-strict: keep + ). +uiu(t_lU(L,U),X,Lin,Bound,Sold) :- + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + true % larger than upperbound: keep + ; + TestUB > 1e-010 -> + + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + fail % smaller than lowerbound: fail + ; + TestBL > 1e-010 -> + + % larger than lowerbound: see t_U case for rest + Strict is Sold /\ 2, + ( + get_attr(X,itf3,(_,_,_,order(OrdX),class(ClassX),_)), + lb(ClassX,OrdX,Vlb-Vb-Lb), + Bound - (Lb+U) < 1e-010 -> + + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lU(L,Bound)),strictness(Strict),RAtt)), + pivot_a(Vlb,X,Vb,t_lu(L,Bound)), + reconsider(X) + ; + get_attr(X,itf3,(_,_,Li,order(OrdX),class(ClassX),RAtt)), + put_attr(X,itf3,(type(t_lU(L,Bound)),strictness(Strict),Li,order(OrdX),class(ClassX),RAtt)), + Delta is Bound-U, + backsubst_delta(ClassX,OrdX,X,Delta) + ) + ; + % equal to lowerbound: check strictness and solve + Sold /\ 2 =:= 0, + solve_bound(Lin,Bound) + ) + ; + true % equal to upperbound and non-strict: keep + % smaller than upperbound: check lowerbound + ). + +% uius(Type,X,Lin,Bound,Strict) +% +% Updates upper bound of independent variable X with linear equation +% Lin that had type Type and strictness Strict, to the new strict +% bound Bound. (see also uiu/5) + +uius(t_none,X,_Lin,Bound,_Sold) :- + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_u(Bound)),strictness(1),RAtt)). +uius(t_u(U),X,_Lin,Bound,Sold) :- + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + true + ; + TestUB > 1e-010 -> + + Strict is Sold \/ 1, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_u(Bound)),strictness(Strict),RAtt)) + ; + Strict is Sold \/ 1, + get_attr(X,itf3,(Ty,_,RAtt)), + put_attr(X,itf3,(Ty,strictness(Strict),RAtt)) + ). +uius(t_l(L),X,_Lin,Bound,Sold) :- + Bound - L > 1e-010, + Strict is Sold \/ 1, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lu(L,Bound)),strictness(Strict),RAtt)). +uius(t_L(L),X,_Lin,Bound,Sold) :- + Bound - L > 1e-010, + Strict is Sold \/ 1, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_Lu(L,Bound)),strictness(Strict),RAtt)). +uius(t_lu(L,U),X,_Lin,Bound,Sold) :- + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + true + ; + TestUB > 1e-010 -> + + Bound - L > 1e-010, + Strict is Sold \/ 1, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lu(L,Bound)),strictness(Strict),RAtt)) + ; + Strict is Sold \/ 1, + get_attr(X,itf3,(Ty,_,RAtt)), + put_attr(X,itf3,(Ty,strictness(Strict),RAtt)) + ). +uius(t_Lu(L,U),X,_Lin,Bound,Sold) :- + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + true + ; + TestUB > 1e-010 -> + + Bound - L > 1e-010, + Strict is Sold \/ 1, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_Lu(L,Bound)),strictness(Strict),RAtt)) + ; + Strict is Sold \/ 1, + get_attr(X,itf3,(Ty,_,RAtt)), + put_attr(X,itf3,(Ty,strictness(Strict),RAtt)) + ). +uius(t_U(U),X,_Lin,Bound,Sold) :- + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + true + ; + TestUB > 1e-010 -> + + Strict is Sold \/ 1, + ( + get_attr(X,itf3,(_,_,_,order(OrdX),class(ClassX),_)), + lb(ClassX,OrdX,Vlb-Vb-Lb), + Bound - (Lb+U) < 1e-010 -> + + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_U(Bound)),strictness(Strict),RAtt)), + pivot_a(Vlb,X,Vb,t_u(Bound)), + reconsider(X) + ; + get_attr(X,itf3,(_,_,Li,order(OrdX),class(ClassX),RAtt)), + put_attr(X,itf3,(type(t_U(Bound)),strictness(Strict),Li,order(OrdX),class(ClassX),RAtt)), + Delta is Bound-U, + backsubst_delta(ClassX,OrdX,X,Delta) + ) + ; + Strict is Sold \/ 1, + get_attr(X,itf3,(Ty,_,RAtt)), + put_attr(X,itf3,(Ty,strictness(Strict),RAtt)) + ). +uius(t_lU(L,U),X,_Lin,Bound,Sold) :- + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + true + ; + TestUB > 1e-010 -> + + Bound - L > 1e-010, + Strict is Sold \/ 1, + ( + get_attr(X,itf3,(_,_,_,order(OrdX),class(ClassX),_)), + lb(ClassX,OrdX,Vlb-Vb-Lb), + Bound - (Lb+U) < 1e-010 -> + + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lU(L,Bound)),strictness(Strict),RAtt)), + pivot_a(Vlb,X,Vb,t_lu(L,Bound)), + reconsider(X) + ; + get_attr(X,itf3,(_,_,Li,order(OrdX),class(ClassX),RAtt)), + put_attr(X,itf3,(type(t_lU(L,Bound)),strictness(Strict),Li,order(OrdX),class(ClassX),RAtt)), + Delta is Bound-U, + backsubst_delta(ClassX,OrdX,X,Delta) + ) + ; + Strict is Sold \/ 1, + get_attr(X,itf3,(Ty,_,RAtt)), + put_attr(X,itf3,(Ty,strictness(Strict),RAtt)) + ). + +% uil(Type,X,Lin,Bound,Strict) +% +% Updates lower bound of independent variable X with linear equation +% Lin that had type Type and strictness Strict, to the new non-strict +% bound Bound. (see also uiu/5) + + +uil(t_none,X,_Lin,Bound,_Sold) :- + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_l(Bound)),strictness(0),RAtt)). +uil(t_l(L),X,_Lin,Bound,Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + true + ; + TestBL > 1e-010 -> + + Strict is Sold /\ 1, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_l(Bound)),strictness(Strict),RAtt)) + ; + true + ). +uil(t_u(U),X,Lin,Bound,_Sold) :- + TestUB is U - Bound, + ( + TestUB < -1e-010 -> + + fail + ; + TestUB > 1e-010 -> + + get_attr(X,itf3,(_,RAtt)), + put_attr(X,itf3,(type(t_lu(Bound,U)),RAtt)) + ; + solve_bound(Lin,Bound) + ). +uil(t_U(U),X,Lin,Bound,_Sold) :- + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + fail + ; + TestUB > 1e-010 -> + + get_attr(X,itf3,(_,RAtt)), + put_attr(X,itf3,(type(t_lU(Bound,U)),RAtt)) + ; + solve_bound(Lin,Bound) + ). +uil(t_lu(L,U),X,Lin,Bound,Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + true + ; + TestBL > 1e-010 -> + + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + fail + ; + TestUB > 1e-010 -> + + Strict is Sold /\ 1, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lu(Bound,U)),strictness(Strict),RAtt)) + ; + Sold /\ 1 =:= 0, + solve_bound(Lin,Bound) + ) + ; + true + ). +uil(t_lU(L,U),X,Lin,Bound,Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + true + ; + TestBL > 1e-010 -> + + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + fail + ; + TestUB > 1e-010 -> + + Strict is Sold /\ 1, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lU(Bound,U)),strictness(Strict),RAtt)) + ; + Sold /\ 1 =:= 0, + solve_bound(Lin,Bound) + ) + ; + true + ). +uil(t_L(L),X,_Lin,Bound,Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + true + ; + TestBL > 1e-010 -> + + Strict is Sold /\ 1, + ( + get_attr(X,itf3,(_,_,_,order(OrdX),class(ClassX),_)), + ub(ClassX,OrdX,Vub-Vb-Ub), + Bound - (Ub + L) > -1e-010 -> + + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_L(Bound)),strictness(Strict),RAtt)), + pivot_a(Vub,X,Vb,t_l(Bound)), + reconsider(X) + ; + get_attr(X,itf3,(_,_,Li,order(OrdX),class(ClassX),RAtt)), + put_attr(X,itf3,(type(t_L(Bound)),strictness(Strict),Li,order(OrdX),class(ClassX),RAtt)), + Delta is Bound-L, + backsubst_delta(ClassX,OrdX,X,Delta) + ) + ; + true + ). +uil(t_Lu(L,U),X,Lin,Bound,Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + true + ; + TestBL > 1e-010 -> + + TestUB is U-Bound, + ( + TestUB < -1e-010 -> + + fail + ; + TestUB > 1e-010 -> + + Strict is Sold /\ 1, + ( + get_attr(X,itf3,(_,_,_,order(OrdX),class(ClassX),_)), + ub(ClassX,OrdX,Vub-Vb-Ub), + Bound - (Ub + L) > -1e-010 -> + + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_Lu(Bound,U)),strictness(Strict),RAtt)), + pivot_a(Vub,X,Vb,t_lu(Bound,U)), + reconsider(X) + ; + get_attr(X,itf3,(_,_,Li,order(OrdX),class(ClassX),RAtt)), + put_attr(X,itf3,(type(t_Lu(Bound,U)),strictness(Strict),Li,order(OrdX),class(ClassX),RAtt)), + Delta is Bound-L, + backsubst_delta(ClassX,OrdX,X,Delta) + ) + ; + Sold /\ 1 =:= 0, + solve_bound(Lin,Bound) + ) + ; + true + ). + +% uils(Type,X,Lin,Bound,Strict) +% +% Updates lower bound of independent variable X with linear equation +% Lin that had type Type and strictness Strict, to the new strict +% bound Bound. (see also uiu/5) + +uils(t_none,X,_Lin,Bound,_Sold) :- + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_l(Bound)),strictness(2),RAtt)). +uils(t_l(L),X,_Lin,Bound,Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + true + ; + TestBL > 1e-010 -> + + Strict is Sold \/ 2, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_l(Bound)),strictness(Strict),RAtt)) + + ; + Strict is Sold \/ 2, + get_attr(X,itf3,(Ty,_,RAtt)), + put_attr(X,itf3,(Ty,strictness(Strict),RAtt)) + ). +uils(t_u(U),X,_Lin,Bound,Sold) :- + U - Bound > 1e-010, + Strict is Sold \/ 2, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lu(Bound,U)),strictness(Strict),RAtt)). +uils(t_U(U),X,_Lin,Bound,Sold) :- + U - Bound > 1e-010, + Strict is Sold \/ 2, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lU(Bound,U)),strictness(Strict),RAtt)). +uils(t_lu(L,U),X,_Lin,Bound,Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + true + ; + TestBL > 1e-010 -> + + U - Bound > 1e-010, + Strict is Sold \/ 2, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lu(Bound,U)),strictness(Strict),RAtt)) + ; + Strict is Sold \/ 2, + get_attr(X,itf3,(Ty,_,RAtt)), + put_attr(X,itf3,(Ty,strictness(Strict),RAtt)) + ). +uils(t_lU(L,U),X,_Lin,Bound,Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + true + ; + TestBL > 1e-010 -> + + U - Bound > 1e-010, + Strict is Sold \/ 2, + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_lU(Bound,U)),strictness(Strict),RAtt)) + ; + Strict is Sold \/ 2, + get_attr(X,itf3,(Ty,_,RAtt)), + put_attr(X,itf3,(Ty,strictness(Strict),RAtt)) + ). +uils(t_L(L),X,_Lin,Bound,Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + true + ; + TestBL > 1e-010 -> + + Strict is Sold \/ 2, + ( + get_attr(X,itf3,(_,_,_,order(OrdX),class(ClassX),_)), + ub(ClassX,OrdX,Vub-Vb-Ub), + Bound - (Ub + L) > -1e-010 -> + + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_L(Bound)),strictness(Strict),RAtt)), + pivot_a(Vub,X,Vb,t_l(Bound)), + reconsider(X) + ; + get_attr(X,itf3,(_,_,Li,order(OrdX),class(ClassX),RAtt)), + put_attr(X,itf3,(type(t_L(Bound)),strictness(Strict),Li,order(OrdX),class(ClassX),RAtt)), + Delta is Bound-L, + backsubst_delta(ClassX,OrdX,X,Delta) + ) + ; + Strict is Sold \/ 2, + get_attr(X,itf3,(Ty,_,RAtt)), + put_attr(X,itf3,(Ty,strictness(Strict),RAtt)) + ). +uils(t_Lu(L,U),X,_Lin,Bound,Sold) :- + TestBL is Bound-L, + ( + TestBL < -1e-010 -> + + true + ; + TestBL > 1e-010 -> + + U - Bound > 1e-010, + Strict is Sold \/ 2, + ( + get_attr(X,itf3,(_,_,_,order(OrdX),class(ClassX),_)), + ub(ClassX,OrdX,Vub-Vb-Ub), + Bound - (Ub + L) > -1e-010 -> + + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_Lu(Bound,U)),strictness(Strict),RAtt)), + pivot_a(Vub,X,Vb,t_lu(Bound,U)), + reconsider(X) + ; + get_attr(X,itf3,(_,_,Li,order(OrdX),class(ClassX),RAtt)), + put_attr(X,itf3,(type(t_Lu(Bound,U)),strictness(Strict),Li,order(OrdX),class(ClassX),RAtt)), + Delta is Bound-L, + backsubst_delta(ClassX,OrdX,X,Delta) + ) + ; + Strict is Sold \/ 2, + get_attr(X,itf3,(Ty,_,RAtt)), + put_attr(X,itf3,(Ty,strictness(Strict),RAtt)) + ). + +% reconsider_upper(X,Lin,U) +% +% Checks if the upperbound of X which is U, satisfies the bounds +% of the variables in Lin: let R be the sum of all the bounds on +% the variables in Lin, and I be the inhomogene part of Lin, then +% upperbound U should be larger or equal to R + I (R may contain +% lowerbounds). +% See also rcb/3 in bv.pl + +reconsider_upper(X,Lin,U) :- + Lin = [I,R|H], + R + I - U > -1e-010, % violation + !, + dec_step(H,Status), % we want to decrement R + rcbl_status(Status,X,[],Binds,[],u(U)), + export_binding(Binds). +reconsider_upper( _, _, _). + +% reconsider_lower(X,Lin,L) +% +% Checks if the lowerbound of X which is L, satisfies the bounds +% of the variables in Lin: let R be the sum of all the bounds on +% the variables in Lin, and I be the inhomogene part of Lin, then +% lowerbound L should be smaller or equal to R + I (R may contain +% upperbounds). +% See also rcb/3 in bv.pl + +reconsider_lower(X,Lin,L) :- + Lin = [I,R|H], + R+I - L < 1e-010, % violation + !, + inc_step(H,Status), % we want to increment R + rcbl_status(Status,X,[],Binds,[],l(L)), + export_binding(Binds). +reconsider_lower(_,_,_). + +% +% lin is dereferenced +% + +% solve_bound(Lin,Bound) +% +% Solves the linear equation Lin - Bound = 0 +% Lin is the linear equation of X, a variable whose bounds have narrowed to value Bound + +solve_bound(Lin,Bound) :- + Bound >= -1e-010, + Bound =< 1e-010, + !, + solve(Lin). +solve_bound(Lin,Bound) :- + Nb is -Bound, + normalize_scalar(Nb,Nbs), + add_linear_11(Nbs,Lin,Eq), + solve(Eq). diff --git a/LGPL/clpr/clpr/itf3.pl b/LGPL/clpr/clpr/itf3.pl new file mode 100644 index 000000000..2d00db996 --- /dev/null +++ b/LGPL/clpr/clpr/itf3.pl @@ -0,0 +1,399 @@ +/* $Id: itf3.pl,v 1.1 2005-10-28 17:51:01 vsc Exp $ + + Part of CPL(R) (Constraint Logic Programming over Reals) + + Author: Leslie De Koninck + E-mail: Tom.Schrijvers@cs.kuleuven.ac.be + WWW: http://www.swi-prolog.org + http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 + Copyright (C): 2004, K.U. Leuven and + 1992-1995, Austrian Research Institute for + Artificial Intelligence (OFAI), + Vienna, Austria + + This software is part of Leslie De Koninck's master thesis, supervised + by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) + by Christian Holzbaur for SICStus Prolog and distributed under the + license details below with permission from all mentioned authors. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 + of the License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + As a special exception, if you link this library with other files, + compiled with a Free Software compiler, to produce an executable, this + library does not by itself cause the resulting executable to be covered + by the GNU General Public License. This exception does not however + invalidate any other reasons why the executable file might be covered by + the GNU General Public License. +*/ + + +% attribute = (type(_),strictness(_),lin(_),order(_),class(_),forward(_),nonzero, +% target,keep_indep,keep) + +:- module(itf3, + [ + dump_linear/4, + dump_nonzero/4 + ]). +:- use_module(arith_r). +:- use_module(bv, + [ + deref/2, + detach_bounds_vlv/5, + dump_var/6, + dump_nz/5, + solve/1, + solve_ord_x/3 + ]). +:- use_module(nf, + [ + nf/2 + ]). +:- use_module(store, + [ + add_linear_11/3, + indep/2, + nf_coeff_of/3 + ]). +:- use_module(class, + [ + class_drop/2 + ]). + +this_linear_solver(clpr). + +% +% Parametrize the answer presentation mechanism +% (toplevel,compiler/debugger ...) +% +:- dynamic presentation_context/1. + +% replaces the old presentation context by a new one. +presentation_context(Old,New) :- + clause(presentation_context(Current), _), + !, + Current = Old, + retractall(presentation_context(_)), + assert(presentation_context(New)). +presentation_context(toplevel,New) :- assert(presentation_context(New)). %default + +% +% attribute_goal( V, V:Atts) :- get_atts( V, Atts). +% +attribute_goal(V,Goal) :- + + presentation_context(Cont,Cont), + dump_linear(V,Cont,Goals,Gtail), + dump_nonzero(V,Cont,Gtail,[]), + l2wrapped(Goals,Goal). + +l2wrapped([],true). +l2wrapped([X|Xs],Conj) :- + ( + Xs = [], + wrap(X,Conj) + ; + Xs = [_|_], + wrap(X,Xw), + Conj = (Xw,Xc), + l2wrapped(Xs,Xc) + ). + +% +% Tests should be pulled out of the loop ... +% +wrap(C,W) :- + prolog_flag(typein_module,Module), + this_linear_solver(Solver), + ( + Module == Solver -> + + W = {C} + ; + predicate_property(Module:{_},imported_from(Solver)) -> + + W = {C} + ; + W = Solver:{C} + ). + +dump_linear(V,Context) --> + { + get_attr(V,itf3,(type(Type),_,lin(Lin),_)), + !, + Lin = [I,_|H] + }, + % + % This happens if not all target variables can be made independend + % Example: examples/option.pl: + % | ?- go2(S,W). + % + % W = 21/4, + % S>=0, + % S<50 ? ; + % + % W>5, + % S=221/4-W, this line would be missing !!! + % W=<21/4 + % + ( + {Type=t_none; get_attr(V,itf3,(_,_,_,_,_,_,_,n,_))} -> + + [] + ; + dump_v(Context,t_none,V,I,H) + ), + % + ( + {Type=t_none, get_attr(V,itf3,(_,_,_,_,_,_,_,n,_)) } -> % nonzero produces such + + [] + ; + dump_v(Context,Type,V,I,H) + ). +dump_linear(_,_) --> []. + +dump_v(toplevel,Type,V,I,H) --> dump_var(Type,V,I,H). +% dump_v(compiler,Type,V,I,H) --> compiler_dump_var(Type,V,I,H). + +dump_nonzero(V,Cont) --> + { + get_attr(V,itf3,(_,_,lin(Lin),_,_,_,nonzero,_)), + !, + Lin = [I,_|H] + }, + dump_nz(Cont,V,H,I). +dump_nonzero(_,_) --> []. + +dump_nz(toplevel,V,H,I) --> dump_nz(V,H,I). +% dump_nz(compiler,V,H,I) --> compiler_dump_nz(V,H,I). + +numbers_only(Y) :- + var(Y), + !. +numbers_only(Y) :- + arith_normalize(Y,Y), + !. +numbers_only(Y) :- + this_linear_solver(Solver), + ( + Solver == clpr -> + + What = 'a real number' + ; + Solver == clpq -> + + What = 'a rational number' + ), + raise_exception(type_error(_X = Y,2,What,Y)). + +attr_unify_hook((n,n,n,n,n,n,n,_),_) :- !. +attr_unify_hook((_,_,_,_,_,forward(F),_),Y) :- + !, + fwd_deref(F,Y). +attr_unify_hook((Ty,St,Li,Or,Cl,_,No,_),Y) :- + % numbers_only(Y), + verify_nonzero(No,Y), + verify_type(Ty,St,Y,Later,[]), + verify_lin(Or,Cl,Li,Y), + call_list(Later). + +% call_list(List) +% +% Calls all the goals in List. + +call_list([]). +call_list([G|Gs]) :- + call(G), + call_list(Gs). + +%%%% +fwd_deref(X,Y) :- + nonvar(X), + X = Y. +fwd_deref(X,Y) :- + var(X), + ( + get_attr(X,itf3,itf(_,_,_,forward(F),_,_,_,_,_,_)) -> + + fwd_deref(F,Y) + ; + X = Y + ). + +% verify_nonzero(Nonzero,Y) +% +% if Nonzero = nonzero, then verify that Y is not zero +% (if possible, otherwise set Y to be nonzero) + +verify_nonzero(nonzero,Y) :- + !, + ( + var(Y) -> + ( + + get_attr(Y,itf3,itf(A,B,C,D,E,F,_,H,I,J)) -> + + put_attr(Y,itf3,itf(A,B,C,D,E,F,nonzero,H,I,J)) + ; + put_attr(Y,itf3,itf(n,n,n,n,n,n,nonzero,n,n,n)) + ) + ; + TestY = Y - 0.0, % Y =\= 0 + ( + TestY < -1e-010 -> + + true + ; + TestY > 1e-010 + ) + ). +verify_nonzero(_,_). % X is not nonzero + +% verify_type(type(Type),strictness(Strict),Y,[OL|OLT],OLT) +% +% if possible verifies whether Y satisfies the type and strictness of X +% if not possible to verify, then returns the constraints that follow from +% the type and strictness + +verify_type(type(Type),strictness(Strict),Y) --> + !, + verify_type2(Y,Type,Strict). +verify_type(_,_,_) --> []. + +verify_type2(Y,TypeX,StrictX) --> + {var(Y)}, + !, + verify_type_var(TypeX,Y,StrictX). +verify_type2(Y,TypeX,StrictX) --> + {verify_type_nonvar(TypeX,Y,StrictX)}. + +% verify_type_nonvar(Type,Nonvar,Strictness) +% +% verifies whether the type and strictness are satisfied with the Nonvar + +verify_type_nonvar(t_none,_,_). % no bounds +verify_type_nonvar(t_l(L),Value,S) :- ilb(S,L,Value). % passive lower bound +verify_type_nonvar(t_u(U),Value,S) :- iub(S,U,Value). % passive upper bound +verify_type_nonvar(t_lu(L,U),Value,S) :- % passive lower and upper bound + ilb(S,L,Value), + iub(S,U,Value). +verify_type_nonvar(t_L(L),Value,S) :- ilb(S,L,Value). % active lower bound +verify_type_nonvar(t_U(U),Value,S) :- iub(S,U,Value). % active upper bound +verify_type_nonvar(t_Lu(L,U),Value,S) :- % active lower and passive upper bound + ilb(S,L,Value), + iub(S,U,Value). +verify_type_nonvar(t_lU(L,U),Value,S) :- % passive lower and active upper bound + ilb(S,L,Value), + iub(S,U,Value). + +% ilb(Strict,Lower,Value) & iub(Strict,Upper,Value) +% +% check whether Value is satisfiable with the given lower/upper bound and strictness +% strictness is encoded as follows: +% 2 = strict lower bound +% 1 = strict upper bound +% 3 = strict lower and upper bound +% 0 = no strict bounds + +ilb(S,L,V) :- + S /\ 2 =:= 0, + !, + L - V < 1e-010. +ilb(_,L,V) :- L - V < -1e-010. + +iub(S,U,V) :- + S /\ 1 =:= 0, + !, + V - U < 1e-010. +iub(_,U,V) :- V - U < -1e-010. + + +% +% Running some goals after X=Y simplifies the coding. It should be possible +% to run the goals here and taking care not to put_atts/2 on X ... +% + +% verify_type_var(Type,Var,Strictness,[OutList|OutListTail],OutListTail) +% +% returns the inequalities following from a type and strictness satisfaction test with Var + +verify_type_var(t_none,_,_) --> []. % no bounds, no inequalities +verify_type_var(t_l(L),Y,S) --> llb(S,L,Y). % passive lower bound +verify_type_var(t_u(U),Y,S) --> lub(S,U,Y). % passive upper bound +verify_type_var(t_lu(L,U),Y,S) --> % passive lower and upper bound + llb(S,L,Y), + lub(S,U,Y). +verify_type_var(t_L(L),Y,S) --> llb(S,L,Y). % active lower bound +verify_type_var(t_U(U),Y,S) --> lub(S,U,Y). % active upper bound +verify_type_var(t_Lu(L,U),Y,S) --> % active lower and passive upper bound + llb(S,L,Y), + lub(S,U,Y). +verify_type_var(t_lU(L,U),Y,S) --> % passive lower and active upper bound + llb(S,L,Y), + lub(S,U,Y). + +% llb(Strict,Lower,Value,[OL|OLT],OLT) and lub(Strict,Upper,Value,[OL|OLT],OLT) +% +% returns the inequalities following from the lower and upper bounds and the strictness +% see also lb and ub +llb(S,L,V) --> + {S /\ 2 =:= 0}, + !, + [{L =< V}]. +llb(_,L,V) --> [{L < V}]. + +lub(S,U,V) --> + {S /\ 1 =:= 0}, + !, + [{V =< U}]. +lub(_,U,V) --> [{V < U}]. + + +% +% We used to drop X from the class/basis to avoid trouble with subsequent +% put_atts/2 on X. Now we could let these dead but harmless updates happen. +% In R however, exported bindings might conflict, e.g. 0 \== 0.0 +% +% If X is indep and we do _not_ solve for it, we are in deep shit +% because the ordering is violated. +% +verify_lin(order(OrdX),class(Class),lin(LinX),Y) :- + !, + ( + indep(LinX,OrdX) -> + + detach_bounds_vlv(OrdX,LinX,Class,Y,NewLinX), % if there were bounds, they are requeued already + class_drop(Class,Y), + nf(-Y,NfY), + deref(NfY,LinY), + add_linear_11(NewLinX,LinY,Lind), + ( + nf_coeff_of(Lind,OrdX,_) -> % X is element of Lind + + solve_ord_x(Lind,OrdX,Class) + ; + solve(Lind) % X is gone, can safely solve Lind + ) + ; + class_drop(Class,Y), + nf(-Y,NfY), + deref(NfY,LinY), + add_linear_11(LinX,LinY,Lind), + solve(Lind) + ). +verify_lin(_,_,_,_). + + diff --git a/LGPL/clpr/clpr/nf.pl b/LGPL/clpr/clpr/nf.pl new file mode 100644 index 000000000..0788b8a28 --- /dev/null +++ b/LGPL/clpr/clpr/nf.pl @@ -0,0 +1,1251 @@ +/* $Id: nf.pl,v 1.1 2005-10-28 17:51:01 vsc Exp $ + + Part of CPL(R) (Constraint Logic Programming over Reals) + + Author: Leslie De Koninck + E-mail: Tom.Schrijvers@cs.kuleuven.ac.be + WWW: http://www.swi-prolog.org + http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 + Copyright (C): 2004, K.U. Leuven and + 1992-1995, Austrian Research Institute for + Artificial Intelligence (OFAI), + Vienna, Austria + + This software is part of Leslie De Koninck's master thesis, supervised + by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) + by Christian Holzbaur for SICStus Prolog and distributed under the + license details below with permission from all mentioned authors. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 + of the License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + As a special exception, if you link this library with other files, + compiled with a Free Software compiler, to produce an executable, this + library does not by itself cause the resulting executable to be covered + by the GNU General Public License. This exception does not however + invalidate any other reasons why the executable file might be covered by + the GNU General Public License. +*/ + + +:- module(nf, + [ + {}/1, + nf/2, + entailed/1, + split/3, + repair/2, + nf_constant/2, + wait_linear/3, + nf2term/2 + + ]). + +:- use_module(geler, + [ + geler/2 + ]). +:- use_module(bv, + [ + export_binding/2, + log_deref/4, + solve/1, + 'solve_<'/1, + 'solve_=<'/1, + 'solve_=\\='/1 + ]). +:- use_module(ineq, + [ + ineq_one/4, + ineq_one_s_p_0/1, + ineq_one_s_n_0/1, + ineq_one_n_p_0/1, + ineq_one_n_n_0/1 + ]). +:- use_module(store, + [ + add_linear_11/3, + normalize_scalar/2 + ]). +:- use_module(nfr, + [ + transg/3 + ]). +:- use_module(arith_r, + [ + arith_eps/1, + arith_normalize/2, + integerp/2 + ]). + + +% ------------------------------------------------------------------------- + +% {Constraint} +% +% Adds the constraint Constraint to the constraint store. +% +% First rule is to prevent binding with other rules when a variable is input +% Constraints are converted to normal form and if necessary, submitted to the linear +% equality/inequality solver (bv + ineq) or to the non-linear store (geler) + +{Rel} :- + var(Rel), + !, + raise_exception(instantiation_error({Rel},1)). +{R,Rs} :- + !, + {R}, + !, + {Rs}. +{R;Rs} :- + !, + ({R};{Rs}). % for entailment checking +{L < R} :- + !, + nf(L-R,Nf), + submit_lt(Nf). +{L > R} :- + !, + nf(R-L,Nf), + submit_lt(Nf). +{L =< R} :- + !, + nf(L-R,Nf), + submit_le( Nf). +{<=(L,R)} :- + !, + nf(L-R,Nf), + submit_le(Nf). +{L >= R} :- + !, + nf(R-L,Nf), + submit_le(Nf). +{L =\= R} :- + !, + nf(L-R,Nf), + submit_ne(Nf). +{L =:= R} :- + !, + nf(L-R,Nf), + submit_eq(Nf). +{L = R} :- + !, + nf(L-R,Nf), + submit_eq(Nf). +{Rel} :- raise_exception(type_error({Rel},1,'a constraint',Rel)). + +% entailed(C) +% +% s -> c = ~s v c = ~(s /\ ~c) +% where s is the store and c is the constraint for which +% we want to know whether it is entailed. +% C is negated and added to the store. If this fails, then c is entailed by s + +entailed(C) :- + negate(C,Cn), + \+ {Cn}. + +% negate(C,Res). +% +% Res is the negation of constraint C +% first rule is to prevent binding with other rules when a variable is input + +negate(Rel,_) :- + var(Rel), + !, + raise_exception(instantiation_error(entailed(Rel),1)). +negate((A,B),(Na;Nb)) :- + !, + negate(A,Na), + negate(B,Nb). +negate((A;B),(Na,Nb)) :- + !, + negate(A,Na), + negate(B,Nb). +negate(A=B) :- !. +negate(A>B,A=B) :- !. +negate(A>=B,A A = 0 +% b4) nonlinear -> geler +% c) Nf=[A,B|Rest] +% c1) A=k +% c11) (B=c*X^+1 or B=c*X^-1), Rest=[] -> B=-k/c or B=-c/k +% c12) invertible(A,B) +% c13) linear(B|Rest) +% c14) geler +% c2) linear(Nf) +% c3) nonlinear -> geler + +submit_eq([]). % trivial success: case a +submit_eq([T|Ts]) :- + submit_eq(Ts,T). +submit_eq([],A) :- submit_eq_b(A). % case b +submit_eq([B|Bs],A) :- submit_eq_c(A,B,Bs). % case c + +% submit_eq_b(A) +% +% Handles case b of submit_eq/1 + +% case b1: A is a constant (non-zero) +submit_eq_b(v(_,[])) :- + !, + fail. +% case b2/b3: A is n*X^P => X = 0 +submit_eq_b(v(_,[X^P])) :- + var(X), + P > 0, + !, + Z = 0.0, + export_binding(X,Z). +% case b2: non-linear is invertible: NL(X) = 0 => X - inv(NL)(0) = 0 +submit_eq_b(v(_,[NL^1])) :- + nonvar(NL), + Z = 0.0, + nl_invertible(NL,X,Z,Inv), + !, + nf(-Inv,S), + nf_add(X,S,New), + submit_eq(New). +% case b4: A is non-linear and not invertible => submit equality to geler +submit_eq_b(Term) :- + term_variables(Term,Vs), + geler(Vs,nf:resubmit_eq([Term])). + +% submit_eq_c(A,B,Rest) +% +% Handles case c of submit_eq/1 + +% case c1: A is a constant +submit_eq_c(v(I,[]),B,Rest) :- + !, + submit_eq_c1(Rest,B,I). +% case c2: A,B and Rest are linear +submit_eq_c(A,B,Rest) :- % c2 + A = v(_,[X^1]), + var(X), + B = v(_,[Y^1]), + var(Y), + linear(Rest), + !, + Hom = [A,B|Rest], + % 'solve_='(Hom). + nf_length(Hom,0,Len), + log_deref(Len,Hom,[],HomD), + solve(HomD). +% case c3: A, B or Rest is non-linear => geler +submit_eq_c(A,B,Rest) :- + Norm = [A,B|Rest], + term_variables(Norm,Vs), + geler(Vs,nf:resubmit_eq(Norm)). + +% submit_eq_c1(Rest,B,K) +% +% Handles case c1 of submit_eq/1 + +% case c11: k+cX^1=0 or k+cX^-1=0 +submit_eq_c1([],v(K,[X^P]),I) :- + var(X), + ( + P = 1, + !, + Val is -I/K, + export_binding(X,Val) + ; + P = -1, + !, + Val is -K/I, + export_binding(X,Val) + ). +% case c12: non-linear, invertible: cNL(X)^1+k=0 => inv(NL)(-k/c) = 0 ; +% cNL(X)^-1+k=0 => inv(NL)(-c/k) = 0 +submit_eq_c1([],v(K,[NL^P]),I) :- + nonvar(NL), + ( + P = 1, + Y is -I/K + ; + P = -1, + Y is -K/I + ), + nl_invertible(NL,X,Y,Inv), + !, + nf(-Inv,S), + nf_add(X,S,New), + submit_eq(New). +% case c13: linear: X + Y + Z + c = 0 => +submit_eq_c1(Rest,B,I) :- + B = v(_,[Y^1]), + var(Y), + linear(Rest), + !, + % 'solve_='( [v(I,[]),B|Rest]). + Hom = [B|Rest], + nf_length(Hom,0,Len), + normalize_scalar(I,Nonvar), + log_deref(Len,Hom,[],HomD), + add_linear_11(Nonvar,HomD,LinD), + solve(LinD). +% case c14: other cases => geler +submit_eq_c1(Rest,B,I) :- + Norm = [v(I,[]),B|Rest], + term_variables(Norm,Vs), + geler(Vs,nf:resubmit_eq(Norm)). + +% ----------------------------------------------------------------------- + +% submit_lt(Nf) +% +% Submits the inequality Nf<0 to the constraint store, where Nf is in normal form. + +% 0 < 0 => fail +submit_lt([]) :- fail. +% A + B < 0 +submit_lt([A|As]) :- submit_lt(As,A). + +% submit_lt(As,A) +% +% Does what submit_lt/1 does where Nf = [A|As] + +% v(K,P) < 0 +submit_lt([],v(K,P)) :- submit_lt_b(P,K). +% A + B + Bs < 0 +submit_lt([B|Bs],A) :- submit_lt_c(Bs,A,B). + +% submit_lt_b(P,K) +% +% Does what submit_lt/2 does where A = [v(K,P)] and As = [] + +% c < 0 +submit_lt_b([],I) :- + !, + I - 0.0 < -1e-010. +% cX^1 < 0 : if c < 0 then X > 0, else X < 0 +submit_lt_b([X^1],K) :- + var(X), + !, + ( + 0.0 - K < -1e-010 -> + + ineq_one_s_p_0(X) % X is strictly positive + ; + ineq_one_s_n_0(X) % X is strictly negative + ). +% non-linear => geler +submit_lt_b(P,K) :- + term_variables(P,Vs), + geler(Vs,nf:resubmit_lt([v(K,P)])). + +% submit_lt_c(Bs,A,B) +% +% Does what submit_lt/2 does where As = [B|Bs]. + +% c + kX < 0 => kX < c +submit_lt_c([],A,B) :- + A = v(I,[]), + B = v(K,[Y^1]), + var(Y), + !, + ineq_one(strict,Y,K,I). +% linear < 0 => solve, non-linear < 0 => geler +submit_lt_c(Rest,A,B) :- + Norm = [A,B|Rest], + ( + linear(Norm) -> + + 'solve_<'(Norm) + ; + term_variables(Norm,Vs), + geler(Vs,nf:resubmit_lt(Norm)) + ). + +% submit_le(Nf) +% +% Submits the inequality Nf =< 0 to the constraint store, where Nf is in normal form. +% See also submit_lt/1 + +% 0 =< 0 => success +submit_le([]). +% A + B =< 0 +submit_le([A|As]) :- submit_le(As,A). + +% submit_le(As,A) +% +% See submit_lt/2. This handles less or equal. + +% v(K,P) =< 0 +submit_le([],v(K,P)) :- submit_le_b(P,K). +% A + B + Bs =< 0 +submit_le([B|Bs],A) :- submit_le_c(Bs,A,B). + +% submit_le_b(P,K) +% +% See submit_lt_b/2. This handles less or equal. + +% c =< 0 +submit_le_b([],I) :- + !, + I - 0.0 < 1e-010. +% cX^1 =< 0: if c < 0 then X >= 0, else X =< 0 +submit_le_b([X^1],K) :- + var(X), + !, + ( + 0.0 - K < -1e-010 -> + + ineq_one_n_p_0(X) % X is non-strictly positive + ; + ineq_one_n_n_0(X) % X is non-strictly negative + ). +% cX^P =< 0 => geler +submit_le_b(P,K) :- + term_variables(P,Vs), + geler(Vs,nf:resubmit_le([v(K,P)])). + +% submit_le_c(Bs,A,B) +% +% See submit_lt_c/3. This handles less or equal. + +% c + kX^1 =< 0 => kX =< 0 +submit_le_c([],A,B) :- + A = v(I,[]), + B = v(K,[Y^1]), + var(Y), + !, + ineq_one(nonstrict,Y,K,I). +% A, B & Rest are linear => solve, otherwise => geler +submit_le_c(Rest,A,B) :- + Norm = [A,B|Rest], + ( + linear(Norm) -> + + 'solve_=<'(Norm) + ; + term_variables(Norm,Vs), + geler(Vs,nf:resubmit_le(Norm)) + ). + +% submit_ne(Nf) +% +% Submits the inequality Nf =\= 0 to the constraint store, where Nf is in normal form. +% if Nf is a constant => check constant = 0, else if Nf is linear => solve else => geler + +submit_ne(Norm1) :- + ( + nf_constant(Norm1,K) -> + + TestK is K - 0.0, % K =\= 0 + ( + TestK < -1e-010 -> + + true + ; + TestK > 1e-010 + ) + ; + linear(Norm1) -> + + 'solve_=\\='(Norm1) + ; + term_variables(Norm1,Vs), + geler(Vs,nf:resubmit_ne(Norm1)) + ). + +% linear(A) +% +% succeeds when A is linear: all elements are of the form v(_,[]) or v(_,[X^1]) + +linear([]). +linear(v(_,Ps)) :- linear_ps(Ps). +linear([A|As]) :- + linear(A), + linear(As). + +% linear_ps(A) +% +% Succeeds when A = V^1 with V a variable. This reflects the linearity of v(_,A). + +linear_ps([]). +linear_ps([V^1]) :- var(V). % excludes sin(_), ... + +% +% Goal delays until Term gets linear. +% At this time, Var will be bound to the normalform of Term. +% +:- meta_predicate wait_linear( ?, ?, :). +% +wait_linear(Term,Var,Goal) :- + nf(Term,Nf), + ( + linear(Nf) -> + + Var = Nf, + call(Goal) + ; + term_variables(Nf,Vars), + geler(Vars,nf:wait_linear_retry(Nf,Var,Goal)) + ). +% +% geler clients +% +resubmit_eq(N) :- + repair(N,Norm), + submit_eq(Norm). +resubmit_lt(N) :- + repair(N,Norm), + submit_lt(Norm). +resubmit_le(N) :- + repair(N,Norm), + submit_le(Norm). +resubmit_ne(N) :- + repair(N,Norm), + submit_ne(Norm). +wait_linear_retry(Nf0,Var,Goal) :- + repair(Nf0,Nf), + ( + linear(Nf) -> + + Var = Nf, + call(Goal) + ; + term_variables(Nf,Vars), + geler(Vars,nf:wait_linear_retry(Nf,Var,Goal)) + ). +% ----------------------------------------------------------------------- + +% nl_invertible(F,X,Y,Res) +% +% Res is the evaluation of the inverse of nonlinear function F in variable X where X is Y + +nl_invertible(sin(X),X,Y,Res) :- Res is asin(Y). +nl_invertible(cos(X),X,Y,Res) :- Res is acos(Y). +nl_invertible(tan(X),X,Y,Res) :- Res is atan(Y). +nl_invertible(exp(B,C),X,A,Res) :- + ( + nf_constant(B,Kb) -> + + 0.0 - A < -1e-010, + 0.0 - Kb < -1e-010, + TestKb is Kb - 1.0, + ( + TestKb < -1e-010 -> + + true + ; + TestKb > 1e-010 + ), + X = C, % note delayed unification + Res is log(A)/log(Kb) + ; + nf_constant(C,Kc), + \+ ( + TestA is A - 0.0, % A =:= 0 + TestA =< 1e-010, + TestA >= -1e-010, + Kc - 0.0 < 1e-010 % Kc =< 0 + ), + X = B, % note delayed unification + Res is A**(1.0/Kc) + ). + +% ----------------------------------------------------------------------- + +% nf(Exp,Nf) +% +% Returns in Nf, the normal form of expression Exp +% +% v(A,[B^C,D^E|...]) means A*B^C*D^E*... where A is a scalar (number) +% v(A,[]) means scalar A + +% variable X => 1*X^1 +nf(X,Norm) :- + var(X), + !, + Norm = [v(One,[X^1])], + One = 1.0. +nf(X,Norm) :- + number(X), + !, + nf_number(X,Norm). +% +nf(rat(N,D),Norm) :- + !, + nf_number(rat(N,D),Norm). +% +nf(#(Const),Norm) :- + monash_constant(Const,Value), + !, + ( + rat(1,1) = 1.0 -> + + nf_number(Value,Norm) % swallows #(zero) ... ok in Q + ; + arith_normalize(Value,N), % in R we want it + Norm = [v(N,[])] + ). +% +nf(-A,Norm) :- + !, + nf(A,An), + K = -1.0, + nf_mul_factor(v(K,[]),An,Norm). +nf(+A,Norm) :- + !, + nf(A,Norm). +% +nf(A+B,Norm) :- + !, + nf(A,An), + nf(B,Bn), + nf_add(An,Bn,Norm). +nf(A-B,Norm) :- + !, + nf(A,An), + nf(-B,Bn), + nf_add(An,Bn,Norm). +% +nf(A*B,Norm) :- + !, + nf(A,An), + nf(B,Bn), + nf_mul(An,Bn,Norm). +nf(A/B,Norm) :- + !, + nf(A,An), + nf(B,Bn), + nf_div(Bn,An,Norm). +% non-linear function, one argument: Term = f(Arg) equals f'(Sa1) = Skel +nf(Term,Norm) :- + nonlin_1(Term,Arg,Skel,Sa1), + !, + nf(Arg,An), + nf_nonlin_1(Skel,An,Sa1,Norm). +% non-linear function, two arguments: Term = f(A1,A2) equals f'(Sa1,Sa2) = Skel +nf(Term,Norm) :- + nonlin_2(Term,A1,A2,Skel,Sa1,Sa2), + !, + nf(A1,A1n), + nf(A2,A2n), + nf_nonlin_2(Skel,A1n,A2n,Sa1,Sa2,Norm). +% +nf(Term,_) :- + raise_exception(type_error(nf(Term,_),1,'a numeric expression',Term)). + +% nf_number(N,Res) +% +% If N is a number, N is normalized (in R: n is converted to a float; in Q: n is converted to a normalized rational) + +nf_number(N,Res) :- + nf_number(N), + arith_normalize(N,Normal), + ( + TestNormal is Normal - 0.0, % Normal =:= 0 + TestNormal =< 1e-010, + TestNormal >= -1e-010 -> + + Res = [] + ; + Res = [v(Normal,[])] + ). + +% checks whether N is a number (standard Prolog number or rational) +nf_number(N) :- + number(N), + !. /* MC 980507 */ +nf_number(N) :- + compound(N), + N=rat(_,_). % sicstus + +nonlin_1(abs(X),X,abs(Y),Y). +nonlin_1(sin(X),X,sin(Y),Y). +nonlin_1(cos(X),X,cos(Y),Y). +nonlin_1(tan(X),X,tan(Y),Y). +nonlin_2(min(A,B),A,B,min(X,Y),X,Y). +nonlin_2(max(A,B),A,B,max(X,Y),X,Y). +nonlin_2(exp(A,B),A,B,exp(X,Y),X,Y). +nonlin_2(pow(A,B),A,B,exp(X,Y),X,Y). % pow->exp +nonlin_2(A^B,A,B,exp(X,Y),X,Y). + +nf_nonlin_1(Skel,An,S1,Norm) :- + ( + nf_constant(An,S1) -> + + nl_eval(Skel,Res), + nf_number(Res,Norm) + ; + S1 = An, + One = 1.0, + Norm = [v(One,[Skel^1])]). +nf_nonlin_2(Skel,A1n,A2n,S1,S2,Norm) :- + ( + nf_constant(A1n,S1), + nf_constant(A2n,S2) -> + + nl_eval(Skel,Res), + nf_number(Res,Norm) + ; + Skel=exp(_,_), + nf_constant(A2n,Exp), + integerp(Exp,I) -> + + nf_power(I,A1n,Norm) + ; + S1 = A1n, + S2 = A2n, + One = 1.0, + Norm = [v(One,[Skel^1])] + ). + +% evaluates non-linear functions in one variable where the variable is bound +nl_eval(abs(X),R) :- R is abs(X). +nl_eval(sin(X),R) :- R is sin(X). +nl_eval(cos(X),R) :- R is cos(X). +nl_eval(tan(X),R) :- R is tan(X). +% evaluates non-linear functions in two variables where both variables are bound +nl_eval(min(X,Y),R) :- R is min(X,Y). +nl_eval(max(X,Y),R) :- R is max(X,Y). +nl_eval(exp(X,Y),R) :- R is X**Y. + +monash_constant(X,_) :- + var(X), + !, + fail. +monash_constant(p,3.14259265). +monash_constant(pi,3.14259265). +monash_constant(e,2.71828182). +monash_constant(zero,Eps) :- arith_eps(Eps). + +% +% check if a Nf consists of just a constant +% + +nf_constant([],Z) :- Z = 0.0. +nf_constant([v(K,[])],K). + +% split(NF,SNF,C) +% +% splits a normalform expression NF into two parts: +% - a constant term C (which might be 0) +% - the homogene part of the expression +% +% this method depends on the polynf ordering, i.e. [] < [X^1] ... + +split([],[],Z) :- Z = 0.0. +split([First|T],H,I) :- + ( + First = v(I,[]) -> + + H = T + ; + I = 0.0, + H = [First|T] + ). + +% nf_add(A,B,C): merges two normalized additions into a new normalized addition +% +% a normalized addition is one where the terms are ordered, e.g. X^1 < Y^1, X^1 < X^2 etc. +% terms in the same variable with the same exponent are added, +% e.g. when A contains v(5,[X^1]) and B contains v(4,[X^1]) then C contains v(9,[X^1]). + +nf_add([],Bs,Bs). +nf_add([A|As],Bs,Cs) :- nf_add(Bs,A,As,Cs). + +nf_add([],A,As,Cs) :- Cs = [A|As]. +nf_add([B|Bs],A,As,Cs) :- + A = v(Ka,Pa), + B = v(Kb,Pb), + compare(Rel,Pa,Pb), + nf_add_case(Rel,A,As,Cs,B,Bs,Ka,Kb,Pa). + +% nf_add_case(Rel,A,As,Cs,B,Bs,Ka,Kb,Pa) +% +% merges sorted lists [A|As] and [B|Bs] into new sorted list Cs +% A = v(Ka,Pa) and B = v(Kb,_) +% Rel is the ordering relation (<, > or =) between A and B. +% when Rel is =, Ka and Kb are added to form a new scalar for Pa +nf_add_case(<,A,As,Cs,B,Bs,_,_,_) :- + Cs = [A|Rest], + nf_add(As,B,Bs,Rest). +nf_add_case(>,A,As,Cs,B,Bs,_,_,_) :- + Cs = [B|Rest], + nf_add(Bs,A,As,Rest). +nf_add_case(=,_,As,Cs,_,Bs,Ka,Kb,Pa) :- + Kc is Ka + Kb, + ( + TestKc = Kc - 0.0, % Kc =:= 0 + TestKc =< 1e-010, + TestKc >= -1e-010 -> + + nf_add(As,Bs,Cs) + ; + Cs = [v(Kc,Pa)|Rest], + nf_add(As,Bs,Rest) + ). + +nf_mul(A,B,Res) :- + nf_length(A,0,LenA), + nf_length(B,0,LenB), + nf_mul_log(LenA,A,[],LenB,B,Res). + +nf_mul_log(0,As,As,_,_,[]) :- !. +nf_mul_log(1,[A|As],As,Lb,B,R) :- + !, + nf_mul_factor_log(Lb,B,[],A,R). +nf_mul_log(2,[A1,A2|As],As,Lb,B,R) :- + !, + nf_mul_factor_log(Lb,B,[],A1,A1b), + nf_mul_factor_log(Lb,B,[],A2,A2b), + nf_add(A1b,A2b,R). +nf_mul_log(N,A0,A2,Lb,B,R) :- + P is N>>1, + Q is N-P, + nf_mul_log(P,A0,A1,Lb,B,Rp), + nf_mul_log(Q,A1,A2,Lb,B,Rq), + nf_add(Rp,Rq,R). + + +% nf_add_2: does the same thing as nf_add, but only has 2 elements to combine. +nf_add_2(Af,Bf,Res) :- % unfold: nf_add([Af],[Bf],Res). + Af = v(Ka,Pa), + Bf = v(Kb,Pb), + compare(Rel,Pa,Pb), + nf_add_2_case(Rel,Af,Bf,Res,Ka,Kb,Pa). + +nf_add_2_case(<,Af,Bf,[Af,Bf],_,_,_). +nf_add_2_case(>,Af,Bf,[Bf,Af],_,_,_). +nf_add_2_case(=,_, _,Res,Ka,Kb,Pa) :- + Kc is Ka + Kb, + ( + TestKc is Kc - 0.0, % Kc =:= 0 + TestKc =< 1e-010, + TestKc >= -1e-010 -> + + Res = [] + ; + Res=[v(Kc,Pa)] + ). + +% nf_mul_k(A,B,C) +% +% C is the result of the multiplication of each element of A (of the form v(_,_)) with scalar B (which shouldn't be 0) +nf_mul_k([],_,[]). +nf_mul_k([v(I,P)|Vs],K,[v(Ki,P)|Vks]) :- + Ki is K*I, + nf_mul_k(Vs,K,Vks). + +% nf_mul_factor(A,Sum,Res) +% +% multiplies each element of the list Sum with factor A which is of the form v(_,_) +% and puts the result in the sorted list Res. +nf_mul_factor(v(K,[]),Sum,Res) :- + !, + nf_mul_k(Sum,K,Res). +nf_mul_factor(F,Sum,Res) :- + nf_length(Sum,0,Len), + nf_mul_factor_log(Len,Sum,[],F,Res). + +% nf_mul_factor_log(Len,[Sum|SumTail],SumTail,F,Res) +% +% multiplies each element of Sum with F and puts the result in the sorted list Res +% Len is the length of Sum +% Sum is split logarithmically each step + +nf_mul_factor_log(0,As,As,_,[]) :- !. +nf_mul_factor_log(1,[A|As],As,F,[R]) :- + !, + mult(A,F,R). +nf_mul_factor_log(2,[A,B|As],As,F,Res) :- + !, + mult(A,F,Af), + mult(B,F,Bf), + nf_add_2(Af,Bf,Res). +nf_mul_factor_log(N,A0,A2,F,R) :- + P is N>>1, % P is rounded(N/2) + Q is N-P, + nf_mul_factor_log(P,A0,A1,F,Rp), + nf_mul_factor_log(Q,A1,A2,F,Rq), + nf_add(Rp,Rq,R). + +% mult(A,B,C) +% +% multiplies A and B into C each of the form v(_,_) + +mult(v(Ka,La),v(Kb,Lb),v(Kc,Lc)) :- + Kc is Ka*Kb, + pmerge(La,Lb,Lc). + +% pmerge(A,B,C) +% +% multiplies A and B into sorted C, where each is of the form of the second argument of v(_,_) + +pmerge([],Bs,Bs). +pmerge([A|As],Bs,Cs) :- pmerge(Bs,A,As,Cs). + +pmerge([],A,As,Res) :- Res = [A|As]. +pmerge([B|Bs],A,As,Res) :- + A = Xa^Ka, + B = Xb^Kb, + compare(R,Xa,Xb), + pmerge_case(R,A,As,Res,B,Bs,Ka,Kb,Xa). + +% pmerge_case(Rel,A,As,Res,B,Bs,Ka,Kb,Xa) +% +% multiplies and sorts [A|As] with [B|Bs] into Res where each is of the form of +% the second argument of v(_,_) +% +% A is Xa^Ka and B is Xb^Kb, Rel is ordening relation between Xa and Xb + +pmerge_case(<,A,As,Res,B,Bs,_,_,_) :- + Res = [A|Tail], + pmerge(As,B,Bs,Tail). +pmerge_case(>,A,As,Res,B,Bs,_,_,_) :- + Res = [B|Tail], + pmerge(Bs,A,As,Tail). +pmerge_case(=,_,As,Res,_,Bs,Ka,Kb,Xa) :- + Kc is Ka + Kb, + ( + Kc=:=0 -> + + pmerge(As,Bs,Res) + ; + Res = [Xa^Kc|Tail], + pmerge(As,Bs,Tail) + ). + +% nf_div(Factor,In,Out) +% +% Out is the result of the division of each element in In (which is of the form v(_,_)) by Factor. + +% division by zero +nf_div([],_,_) :- + !, + zero_division. +% division by v(K,P) => multiplication by v(1/K,P^-1) +nf_div([v(K,P)],Sum,Res) :- + !, + Ki is 1.0/K, + mult_exp(P,-1,Pi), + nf_mul_factor(v(Ki,Pi),Sum,Res). +nf_div(D,A,[v(One,[(A/D)^1])]) :- One = 1.0. + +% zero_division +% +% called when a division by zero is performed +zero_division :- fail. % raise_exception(_) ? + +% mult_exp(In,Factor,Out) +% +% Out is the result of the multiplication of the exponents of the elements in In +% (which are of the form X^Exp by Factor. +mult_exp([],_,[]). +mult_exp([X^P|Xs],K,[X^I|Tail]) :- + I is K*P, + mult_exp(Xs,K,Tail). +% +% raise to integer powers +% +% | ?- time({(1+X+Y+Z)^15=0}). (sicstus, try with SWI) +% Timing 00:00:02.610 2.610 iterative +% Timing 00:00:00.660 0.660 binomial +nf_power(N,Sum,Norm) :- + integer(N), + compare(Rel,N,0), + ( + Rel = (<) -> + + Pn is -N, + % nf_power_pos(Pn,Sum,Inorm), + binom(Sum,Pn,Inorm), + One = 1.0, + nf_div(Inorm,[v(One,[])],Norm) + ; + Rel = (>) -> + + % nf_power_pos(N,Sum,Norm) + binom(Sum,N,Norm) + ; + Rel = (=) -> % 0^0 is indeterminate but we say 1 + + One = 1.0, + Norm = [v(One,[])] + ). +% +% N>0 +% +% iterative method: X^N = X*(X^N-1) +nf_power_pos(1,Sum,Norm) :- + !, + Sum = Norm. +nf_power_pos(N,Sum,Norm) :- + N1 is N-1, + nf_power_pos(N1,Sum,Pn1), + nf_mul(Sum,Pn1,Norm). +% +% N>0 +% +% binomial method +binom(Sum,1,Power) :- + !, + Power = Sum. +binom([],_,[]). +binom([A|Bs],N,Power) :- + ( + Bs = [] -> + + + nf_power_factor(A,N,Ap), + Power = [Ap] + ; + Bs = [_|_] -> + + One = 1.0, + factor_powers(N,A,v(One,[]),Pas), + sum_powers(N,Bs,[v(One,[])],Pbs,[]), + combine_powers(Pas,Pbs,0,N,1,[],Power) + ). + +combine_powers([],[],_,_,_,Pi,Pi). +combine_powers([A|As],[B|Bs],L,R,C,Pi,Po) :- + nf_mul(A,B,Ab), + arith_normalize(C,Cn), + nf_mul_k(Ab,Cn,Abc), + nf_add(Abc,Pi,Pii), + L1 is L+1, + R1 is R-1, + C1 is C*R//L1, + combine_powers(As,Bs,L1,R1,C1,Pii,Po). + +nf_power_factor(v(K,P),N,v(Kn,Pn)) :- + arith_normalize(N,Nn), + Kn is K**Nn, + mult_exp(P,N,Pn). + +factor_powers(0,_,Prev,[[Prev]]) :- !. +factor_powers(N,F,Prev,[[Prev]|Ps]) :- + N1 is N-1, + mult(Prev,F,Next), + factor_powers(N1,F,Next,Ps). +sum_powers(0,_,Prev,[Prev|Lt],Lt) :- !. +sum_powers(N,S,Prev,L0,Lt) :- + N1 is N-1, + nf_mul(S,Prev,Next), + sum_powers(N1,S,Next,L0,[Prev|Lt]). + +% ------------------------------------------------------------------------------ +repair(Sum,Norm) :- + nf_length(Sum,0,Len), + repair_log(Len,Sum,[],Norm). +repair_log(0,As,As,[]) :- !. +repair_log(1,[v(Ka,Pa)|As],As,R) :- + !, + repair_term(Ka,Pa,R). +repair_log(2,[v(Ka,Pa),v(Kb,Pb)|As],As,R) :- + !, + repair_term(Ka,Pa,Ar), + repair_term(Kb,Pb,Br), + nf_add(Ar,Br,R). +repair_log(N,A0,A2,R) :- + P is N>>1, + Q is N-P, + repair_log(P,A0,A1,Rp), + repair_log(Q,A1,A2,Rq), + nf_add(Rp,Rq,R). + +repair_term(K,P,Norm) :- + length(P,Len), + One = 1.0, + repair_p_log(Len,P,[],Pr,[v(One,[])],Sum), + nf_mul_factor(v(K,Pr),Sum,Norm). + +repair_p_log(0,Ps,Ps,[],L0,L0) :- !. +repair_p_log(1,[X^P|Ps],Ps,R,L0,L1) :- + !, + repair_p(X,P,R,L0,L1). +repair_p_log(2,[X^Px,Y^Py|Ps],Ps,R,L0,L2) :- + !, + repair_p(X,Px,Rx,L0,L1), + repair_p(Y,Py,Ry,L1,L2), + pmerge(Rx,Ry,R). +repair_p_log(N,P0,P2,R,L0,L2) :- + P is N>>1, + Q is N-P, + repair_p_log(P,P0,P1,Rp,L0,L1), + repair_p_log(Q,P1,P2,Rq,L1,L2), + pmerge(Rp,Rq,R). + +repair_p(Term,P,[Term^P],L0,L0) :- var(Term). +repair_p(Term,P,[],L0,L1) :- + nonvar(Term), + repair_p_one(Term,TermN), + nf_power(P,TermN,TermNP), + nf_mul(TermNP,L0,L1). +% +% An undigested term a/b is distinguished from an +% digested one by the fact that its arguments are +% digested -> cuts after repair of args! +% +repair_p_one(Term,TermN) :- + nf_number(Term,TermN), % freq. shortcut for nf/2 case below + !. +repair_p_one(A1/A2,TermN) :- + repair(A1,A1n), + repair(A2,A2n), + !, + nf_div(A2n,A1n,TermN). +repair_p_one(Term,TermN) :- + nonlin_1(Term,Arg,Skel,Sa), + repair(Arg,An), + !, + nf_nonlin_1(Skel,An,Sa,TermN). +repair_p_one(Term,TermN) :- + nonlin_2(Term,A1,A2,Skel,Sa1,Sa2), + repair(A1,A1n), + repair(A2,A2n), + !, + nf_nonlin_2(Skel,A1n,A2n,Sa1,Sa2,TermN). +repair_p_one(Term,TermN) :- + nf(Term,TermN). + +nf_length([],Li,Li). +nf_length([_|R],Li,Lo) :- + Lii is Li+1, + nf_length(R,Lii,Lo). +% ------------------------------------------------------------------------------ +% nf2term(NF,Term) +% +% transforms a normal form into a readable term + +% empty normal form = 0 +nf2term([],Z) :- Z = 0.0. +% term is first element (+ next elements) +nf2term([F|Fs],T) :- + f02t(F,T0), % first element + yfx(Fs,T0,T). % next elements + +yfx([],T0,T0). +yfx([F|Fs],T0,TN) :- + fn2t(F,Ft,Op), + T1 =.. [Op,T0,Ft], + yfx(Fs,T1,TN). + +% f02t(v(K,P),T) +% +% transforms the first element of the normal form (something of the form v(K,P)) +% into a readable term +f02t(v(K,P),T) :- + ( + P = [] -> % just a constant + + T = K + ; + TestK is K - 0.0, % K =:= 0 + TestK =< 1e-010, + TestK >= -1e-010 -> % constant = 1,P =\= [], don't show constant + + p2term(P,T) + ; + TestK is K - -1.0, % K =:= -1 + TestK =< 1e-010, + TestK >= -1e-010 -> + + T = -Pt, % constant = -1, P =\= [], only show sign (-) p2term(P,Pt) + p2term(P,Pt) + ; + T = K*Pt, + p2term(P,Pt) + ). + +% f02t(v(K,P),T,Op) +% +% transforms a next element of the normal form (something of the form v(K,P)) into a readable term +fn2t(v(K,P),Term,Op) :- + ( + TestK is K - 1.0, % K =:= 1 + TestK =< 1e-010, + TestK >= -1e-010 -> % no constant, Op = + + + Term = Pt, + Op = + + ; + TestK is K - -1.0, % K =:= -1 + TestK =< 1e-010, + TestK >= -1e-010 -> % no constant, Op = - + + Term = Pt, + Op = - + ; + K - 0.0 < -1e-010 -> % constant, Op = - + + Kf is -K, + Term = Kf*Pt, + Op = - + ; % constant, Op = + + Term = K*Pt, + Op = + + ), + p2term(P,Pt). + +% transforms the P part in v(_,P) into a readable term +p2term([X^P|Xs],Term) :- + ( + Xs = [] -> + + pe2term(X,Xt), + exp2term(P,Xt,Term) + ; + Xs = [_|_] -> + + Term = Xst*Xtp, + pe2term(X,Xt), + exp2term(P,Xt,Xtp), + p2term(Xs,Xst) + ). + +% +exp2term(1,X,X) :- !. +exp2term(-1,X,One/X) :- + !, + One = 1.0. +exp2term(P,X,Term) :- + arith_normalize(P,Pn), + % Term = exp(X,Pn). + Term = X^Pn. + +pe2term(X,Term) :- + var(X), + Term = X. +pe2term(X,Term) :- + nonvar(X), + X =.. [F|Args], + pe2term_args(Args,Argst), + Term =.. [F|Argst]. + +pe2term_args([],[]). +pe2term_args([A|As],[T|Ts]) :- + nf2term(A,T), + pe2term_args(As,Ts). diff --git a/LGPL/clpr/clpr/nfr.pl b/LGPL/clpr/clpr/nfr.pl new file mode 100644 index 000000000..7d8bc392c --- /dev/null +++ b/LGPL/clpr/clpr/nfr.pl @@ -0,0 +1,87 @@ +/* $Id: nfr.pl,v 1.1 2005-10-28 17:51:01 vsc Exp $ + + Part of CPL(R) (Constraint Logic Programming over Reals) + + Author: Leslie De Koninck + E-mail: Tom.Schrijvers@cs.kuleuven.ac.be + WWW: http://www.swi-prolog.org + http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 + Copyright (C): 2004, K.U. Leuven and + 1992-1995, Austrian Research Institute for + Artificial Intelligence (OFAI), + Vienna, Austria + + This software is part of Leslie De Koninck's master thesis, supervised + by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) + by Christian Holzbaur for SICStus Prolog and distributed under the + license details below with permission from all mentioned authors. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 + of the License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + As a special exception, if you link this library with other files, + compiled with a Free Software compiler, to produce an executable, this + library does not by itself cause the resulting executable to be covered + by the GNU General Public License. This exception does not however + invalidate any other reasons why the executable file might be covered by + the GNU General Public License. +*/ + + +:- module(nfr, + [ + transg/3 + ]). + +:- use_module(nf, + [ + nf2term/2 + ]). + +% transg(Goal,[OutList|OutListTail],OutListTail) +% +% puts the equalities and inequalities that are implied by the elements in Goal +% in the difference list OutList +% +% called by geler.pl for project.pl + +transg(resubmit_eq(Nf)) --> + { + nf2term([],Z), + nf2term(Nf,Term) + }, + [clpr:{Term=Z}]. +transg(resubmit_lt(Nf)) --> + { + nf2term([],Z), + nf2term(Nf,Term) + }, + [clpr:{Term + { + nf2term([],Z), + nf2term(Nf,Term) + }, + [clpr:{Term= + { + nf2term([],Z), + nf2term(Nf,Term) + }, + [clpr:{Term=\=Z}]. +transg(wait_linear_retry(Nf,Res,Goal)) --> + { + nf2term(Nf,Term) + }, + [clpr:{Term=Res},Goal]. diff --git a/LGPL/clpr/clpr/ordering.pl b/LGPL/clpr/clpr/ordering.pl new file mode 100644 index 000000000..639e55ded --- /dev/null +++ b/LGPL/clpr/clpr/ordering.pl @@ -0,0 +1,210 @@ +/* $Id: ordering.pl,v 1.1 2005-10-28 17:51:01 vsc Exp $ + + Part of CPL(R) (Constraint Logic Programming over Reals) + + Author: Leslie De Koninck + E-mail: Tom.Schrijvers@cs.kuleuven.ac.be + WWW: http://www.swi-prolog.org + http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 + Copyright (C): 2004, K.U. Leuven and + 1992-1995, Austrian Research Institute for + Artificial Intelligence (OFAI), + Vienna, Austria + + This software is part of Leslie De Koninck's master thesis, supervised + by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) + by Christian Holzbaur for SICStus Prolog and distributed under the + license details below with permission from all mentioned authors. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 + of the License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + As a special exception, if you link this library with other files, + compiled with a Free Software compiler, to produce an executable, this + library does not by itself cause the resulting executable to be covered + by the GNU General Public License. This exception does not however + invalidate any other reasons why the executable file might be covered by + the GNU General Public License. +*/ + + +:- module(ordering, + [ + combine/3, + ordering/1, + arrangement/2 + ]). +:- use_module(class, + [ + class_get_prio/2, + class_put_prio/2 + ]). +:- use_module(bv, + [ + var_intern/2 + ]). + +:- use_module(ugraphs, + [ + add_edges/3, + add_vertices/3, + top_sort/2 + ]). + +ordering(X) :- + var(X), + !, + fail. +ordering(A>B) :- + !, + ordering(B + + var_intern(X,Class) + ; + true + ), + join_class(Xs,Class). + +% combine(Ga,Gb,Gc) +% +% Combines the vertices of Ga and Gb into Gc. + +combine(Ga,Gb,Gc) :- + normalize(Ga,Gan), + normalize(Gb,Gbn), + ugraphs:graph_union(Gan,Gbn,Gc). + +% +% both Ga and Gb might have their internal ordering invalidated +% because of bindings and aliasings +% + +normalize([],[]) :- !. +normalize(G,Gsgn) :- + G = [_|_], + keysort(G,Gs), % sort vertices on key + group(Gs,Gsg), % concatenate vertices with the same key + normalize_vertices(Gsg,Gsgn). % normalize + +normalize_vertices([],[]). +normalize_vertices([X-Xnb|Xs],Res) :- + ( + normalize_vertex(X,Xnb,Xnorm) -> + + Res = [Xnorm|Xsn], + normalize_vertices(Xs,Xsn) + ; + normalize_vertices(Xs,Res) + ). + +% normalize_vertex(X,Nbs,X-Nbss) +% +% Normalizes a vertex X-Nbs into X-Nbss by sorting Nbs, removing duplicates (also of X) +% and removing non-vars. + +normalize_vertex(X,Nbs,X-Nbsss) :- + var(X), + sort(Nbs,Nbss), + strip_nonvar(Nbss,X,Nbsss). + +% strip_nonvar(Nbs,X,Res) +% +% Turns vertext X-Nbs into X-Res by removing occurrences of X from Nbs and removing +% non-vars. This to normalize after bindings have occurred. See also normalize_vertex/3. + +strip_nonvar([],_,[]). +strip_nonvar([X|Xs],Y,Res) :- + ( + X==Y -> % duplicate of Y + + strip_nonvar(Xs,Y,Res) + ; + var(X) -> % var: keep + + Res = [X|Stripped], + strip_nonvar(Xs,Y,Stripped) + ; % nonvar: remove + nonvar(X), + Res = [] % because Vars []. +gen_edges([X|Xs]) --> + gen_edges(Xs,X), + gen_edges(Xs). + +gen_edges([],_) --> []. +gen_edges([Y|Ys],X) --> + [X-Y], + gen_edges(Ys,X). + +% group(Vert,Res) +% +% Concatenates vertices with the same key. + +group([],[]). +group([K-Kl|Ks],Res) :- + group(Ks,K,Kl,Res). + +group([],K,Kl,[K-Kl]). +group([L-Ll|Ls],K,Kl,Res) :- + ( + K==L -> + + append(Kl,Ll,KLl), + group(Ls,K,KLl,Res) + ; + Res = [K-Kl|Tail], + group(Ls,L,Ll,Tail) + ). + + + + + diff --git a/LGPL/clpr/clpr/project.pl b/LGPL/clpr/clpr/project.pl new file mode 100644 index 000000000..b5650f0ab --- /dev/null +++ b/LGPL/clpr/clpr/project.pl @@ -0,0 +1,286 @@ +/* $Id: project.pl,v 1.1 2005-10-28 17:51:01 vsc Exp $ + + Part of CPL(R) (Constraint Logic Programming over Reals) + + Author: Leslie De Koninck + E-mail: Tom.Schrijvers@cs.kuleuven.ac.be + WWW: http://www.swi-prolog.org + http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 + Copyright (C): 2004, K.U. Leuven and + 1992-1995, Austrian Research Institute for + Artificial Intelligence (OFAI), + Vienna, Austria + + This software is part of Leslie De Koninck's master thesis, supervised + by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) + by Christian Holzbaur for SICStus Prolog and distributed under the + license details below with permission from all mentioned authors. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 + of the License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + As a special exception, if you link this library with other files, + compiled with a Free Software compiler, to produce an executable, this + library does not by itself cause the resulting executable to be covered + by the GNU General Public License. This exception does not however + invalidate any other reasons why the executable file might be covered by + the GNU General Public License. +*/ + +% +% Answer constraint projection +% + +%:- public project_attributes/2. % xref.pl + +:- module(project, + [ + drop_dep/1, + drop_dep_one/1, + make_target_indep/2, + project_attributes/2 + ]). +:- use_module(class, + [ + class_allvars/2 + ]). +:- use_module(geler, + [ + project_nonlin/3 + + ]). +:- use_module(fourmotz, + [ + fm_elim/3 + ]). +:- use_module(redund, + [ + redundancy_vars/1, + systems/3 + ]). +:- use_module(bv, + [ + pivot/4 + ]). +:- use_module(ordering, + [ + arrangement/2 + ]). +:- use_module(store, + [ + indep/2, + renormalize/2 + ]). + +% +% interface predicate +% +% May be destructive (either acts on a copy or in a failure loop) +% +project_attributes(TargetVars,Cvas) :- + sort(TargetVars,Tvs), % duplicates ? + sort(Cvas,Avs), % duplicates ? + mark_target(Tvs), + project_nonlin(Tvs,Avs,NlReachable), + ( + Tvs == [] -> + + drop_lin_atts(Avs) + ; + redundancy_vars(Avs), % removes redundant bounds (redund.pl) + make_target_indep(Tvs,Pivots), % pivot partners are marked to be kept during elim. + mark_target(NlReachable), % after make_indep to express priority + drop_dep(Avs), + fm_elim(Avs,Tvs,Pivots), + impose_ordering(Avs) + ). + +% mark_target(Vars) +% +% Marks the variables in Vars as target variables. + +mark_target([]). +mark_target([V|Vs]) :- + get_attr(V,itf3,(Ty,St,Li,Or,Cl,Fo,No,_,RAtt)), + put_attr(V,itf3,(Ty,St,Li,Or,Cl,Fo,No,target,RAtt)), + mark_target(Vs). + +% mark_keep(Vars) +% +% Mark the variables in Vars to be kept during elimination. + +mark_keep([]). +mark_keep([V|Vs]) :- + get_attr(V,itf3,(Ty,St,Li,Or,Cl,Fo,No,Ta,KI,_)), + put_attr(V,itf3,(Ty,St,Li,Or,Cl,Fo,No,Ta,KI,keep)), + mark_keep(Vs). + +% +% Collect the pivots in reverse order +% We have to protect the target variables pivot partners +% from redundancy eliminations triggered by fm_elim, +% in order to allow for reverse pivoting. +% +make_target_indep(Ts,Ps) :- make_target_indep(Ts,[],Ps). + +% make_target_indep(Targets,Pivots,PivotsTail) +% +% Tries to make as many targetvariables independent by pivoting them with a non-target +% variable. The pivots are stored as T:NT where T is a target variable and NT a non-target +% variable. The non-target variables are marked to be kept during redundancy eliminations. + +make_target_indep([],Ps,Ps). +make_target_indep([T|Ts],Ps0,Pst) :- + ( + get_attr(T,itf3,(type(Type),_,lin(Lin),_)), + Lin = [_,_|H], + nontarget(H,Nt) -> + + Ps1 = [T:Nt|Ps0], + get_attr(Nt,itf3,(Ty,St,Li,order(Ord),class(Class),Fo,No,Ta,KI,_)), + put_attr(Nt,itf3,(Ty,St,Li,order(Ord),class(Class),Fo,No,Ta,KI,keep)), + pivot(T,Class,Ord,Type) + ; + Ps1 = Ps0 + ), + make_target_indep(Ts,Ps1,Pst). + +% nontarget(Hom,Nt) +% +% Finds a nontarget variable in homogene part Hom. +% Hom contains elements of the form l(V*K,OrdV). +% A nontarget variable has no target attribute and no keep_indep attribute. + +nontarget([l(V*_,_)|Vs],Nt) :- + ( + get_attr(V,itf3,(_,_,_,_,_,_,_,n,n,_)) -> + + Nt = V + ; + nontarget(Vs,Nt) + ). + +% drop_dep(Vars) +% +% Does drop_dep_one/1 on each variable in Vars. + +drop_dep(Vs) :- + var(Vs), + !. +drop_dep([]). +drop_dep([V|Vs]) :- + drop_dep_one(V), + drop_dep(Vs). + +% drop_dep_one(V) +% +% If V is an unbounded dependent variable that isn't a target variable, shouldn't be kept +% and is not nonzero, drops all linear attributes of V. +% The linear attributes are: type, strictness, linear equation (lin), class and order. + +drop_dep_one(V) :- + get_attr(V,itf3,(type(t_none),_,lin(Lin),order(OrdV),_,Fo,n,n,KI,n)), + \+ indep(Lin,OrdV), + !, + put_attr(V,itf3,(n,n,n,n,n,Fo,n,n,KI,n)). +drop_dep_one(_). + +% drop_lin_atts(Vs) +% +% Removes the linear attributes of the variables in Vs. +% The linear attributes are type, strictness, linear equation (lin), order and class. + +drop_lin_atts([]). +drop_lin_atts([V|Vs]) :- + get_attr(V,itf3,(_,_,_,_,_,RAtt)), + put_attr(V,itf3,(n,n,n,n,n,RAtt)), + drop_lin_atts(Vs). + +impose_ordering(Cvas) :- + systems(Cvas,[],Sys), + impose_ordering_sys(Sys). + + + +impose_ordering_sys([]). +impose_ordering_sys([S|Ss]) :- + arrangement(S,Arr), % ordering.pl + arrange(Arr,S), + impose_ordering_sys(Ss). + +arrange([],_). +arrange(Arr,S) :- + Arr = [_|_], + class_allvars(S,All), + order(Arr,1,N), + order(All,N,_), + renorm_all(All), + arrange_pivot(All). + +order(Xs,N,M) :- + var(Xs), + !, + N = M. +order([],N,N). +order([X|Xs],N,M) :- + ( + get_attr(X,itf3,(_,_,_,order(O),_)), + var(O) -> + + O = N, + N1 is N+1, + order(Xs,N1,M) + ; + order(Xs,N,M) + ). + +% renorm_all(Vars) +% +% Renormalizes all linear equations of the variables in difference list Vars to reflect +% their new ordering. + +renorm_all(Xs) :- + var(Xs), + !. +renorm_all([X|Xs]) :- + ( + get_attr(X,itf3,(Ty,St,lin(Lin),RAtt)) -> + + renormalize(Lin,New), + put_attr(X,itf3,(Ty,St,lin(New),RAtt)), + renorm_all(Xs) + ; + renorm_all(Xs) + ). + +% arrange_pivot(Vars) +% +% If variable X of Vars has type t_none and has a higher order than the first element of +% its linear equation, then it is pivoted with that element. + +arrange_pivot(Xs) :- + var(Xs), + !. +arrange_pivot([X|Xs]) :- + ( + get_attr(X,itf3,(type(t_none),_,lin(Lin),order(OrdX),_)), + Lin = [_,_|[l(Y*_,_)|_]], + get_attr(Y,itf3,(_,_,_,order(OrdY),class(Class),_)), + compare(<,OrdY,OrdX) -> + + pivot(X,Class,OrdY,t_none), + arrange_pivot(Xs) + ; + arrange_pivot(Xs) + ). diff --git a/LGPL/clpr/clpr/redund.pl b/LGPL/clpr/clpr/redund.pl new file mode 100644 index 000000000..331d7faaa --- /dev/null +++ b/LGPL/clpr/clpr/redund.pl @@ -0,0 +1,304 @@ +/* $Id: redund.pl,v 1.1 2005-10-28 17:51:01 vsc Exp $ + + Part of CPL(R) (Constraint Logic Programming over Reals) + + Author: Leslie De Koninck + E-mail: Tom.Schrijvers@cs.kuleuven.ac.be + WWW: http://www.swi-prolog.org + http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 + Copyright (C): 2004, K.U. Leuven and + 1992-1995, Austrian Research Institute for + Artificial Intelligence (OFAI), + Vienna, Austria + + This software is part of Leslie De Koninck's master thesis, supervised + by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) + by Christian Holzbaur for SICStus Prolog and distributed under the + license details below with permission from all mentioned authors. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 + of the License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + As a special exception, if you link this library with other files, + compiled with a Free Software compiler, to produce an executable, this + library does not by itself cause the resulting executable to be covered + by the GNU General Public License. This exception does not however + invalidate any other reasons why the executable file might be covered by + the GNU General Public License. +*/ + + +:- module(redund, + [ + redundancy_vars/1, + systems/3 + ]). +:- use_module(bv, + [ + detach_bounds/1, + intro_at/3 + ]). +:- use_module(class, + [ + class_allvars/2 + ]). +:- use_module(nf, [{}/1]). + +% +% redundancy removal (semantic definition) +% +% done: +% +) deal with active bounds +% +) indep t_[lu] -> t_none invalidates invariants (fixed) +% + +% systems(Vars,SystemsIn,SystemsOut) +% +% Returns in SystemsOut the different classes to which variables in Vars +% belong. Every class only appears once in SystemsOut. + +systems([],Si,Si). +systems([V|Vs],Si,So) :- + ( + var(V), + get_attr(V,itf3,(_,_,_,_,class(C),_)), + not_memq(Si,C) -> + + systems(Vs,[C|Si],So) + ; + systems(Vs,Si,So) + ). + +% not_memq(Lst,El) +% +% Succeeds if El is not a member of Lst (doesn't use unification). + +not_memq([],_). +not_memq([Y|Ys],X) :- + X \== Y, + not_memq(Ys,X). + +% redundancy_systems(Classes) +% +% Does redundancy removal via redundancy_vs/1 on all variables in the classes Classes. + +redundancy_systems([]). +redundancy_systems([S|Sys]) :- + class_allvars(S,All), + redundancy_vs(All), + redundancy_systems(Sys). + +% redundancy_vars(Vs) +% +% Does the same thing as redundancy_vs/1 but has some extra timing facilities that +% may be used. + +redundancy_vars(Vs) :- + !, + redundancy_vs(Vs). +redundancy_vars(Vs) :- + statistics(runtime,[Start|_]), + redundancy_vs(Vs), + statistics(runtime,[End|_]), + Duration is End-Start, + format(user_error,"% Redundancy elimination took ~d msec~n",Duration). + + +% redundancy_vs(Vs) +% +% Removes redundant bounds from the variables in Vs via redundant/3 + +redundancy_vs(Vs) :- + var(Vs), + !. +redundancy_vs([]). +redundancy_vs([V|Vs]) :- + ( + get_attr(V,itf3,(type(Type),strictness(Strict),_)), + redundant(Type,V,Strict) -> + + redundancy_vs(Vs) + ; + redundancy_vs(Vs) + ). + +% redundant(Type,Var,Strict) +% +% Removes redundant bounds from variable Var with type Type and strictness Strict. +% A redundant bound is one that is satisfied anyway (so adding the inverse of the bound +% makes the system infeasible. This predicate can either fail or succeed but a success +% doesn't necessarily mean a redundant bound. + +redundant(t_l(L),X,Strict) :- + detach_bounds(X), % drop temporarily + % if not redundant, backtracking will restore bound + negate_l(Strict,L,X), + red_t_l. % negate_l didn't fail, redundant bound +redundant(t_u(U),X,Strict) :- + detach_bounds(X), + negate_u(Strict,U,X), + red_t_u. +redundant(t_lu(L,U),X,Strict) :- + strictness_parts(Strict,Sl,Su), + ( + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_u(U)),strictness(Su),RAtt)), + negate_l(Strict,L,X) -> + red_t_l, + ( + redundant(t_u(U),X,Strict) -> + + true + ; + true + ) + ; + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_l(L)),strictness(Sl),RAtt)), + negate_u(Strict,U,X) -> + + red_t_u + ; + true + ). +redundant(t_L(L),X,Strict) :- + Bound is -L, + intro_at(X,Bound,t_none), % drop temporarily + detach_bounds(X), + negate_l(Strict,L,X), + red_t_L. +redundant(t_U(U),X,Strict) :- + Bound is -U, + intro_at(X,Bound,t_none), % drop temporarily + detach_bounds(X), + negate_u(Strict,U,X), + red_t_U. +redundant(t_Lu(L,U),X,Strict) :- + strictness_parts(Strict,Sl,Su), + ( + Bound is -L, + intro_at(X,Bound,t_u(U)), + get_attr(X,itf3,(Ty,_,RAtt)), + put_attr(X,itf3,(Ty,strictness(Su),RAtt)), + negate_l(Strict,L,X) -> + + red_t_l, + ( + redundant(t_u(U),X,Strict) -> + + true + ; + true + ) + ; + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_L(L)),strictness(Sl),RAtt)), + negate_u(Strict,U,X) -> + + red_t_u + ; + true + ). +redundant(t_lU(L,U),X,Strict) :- + strictness_parts(Strict,Sl,Su), + ( + get_attr(X,itf3,(_,_,RAtt)), + put_attr(X,itf3,(type(t_U(U)),strictness(Su),RAtt)), + negate_l(Strict,L,X) -> + + red_t_l, + ( + redundant(t_U(U),X,Strict) -> + + true + ; + true + ) + ; + Bound is -U, + intro_at(X,Bound,t_l(L)), + get_attr(X,itf3,(Ty,_,RAtt)), + put_attr(X,itf3,(Ty,strictness(Sl),RAtt)), + negate_u(Strict,U,X) -> + + red_t_u + ; + true + ). + +% strictness_parts(Strict,Lower,Upper) +% +% Splits strictness Strict into two parts: one related to the lowerbound and +% one related to the upperbound. + +strictness_parts(Strict,Lower,Upper) :- + Lower is Strict /\ 2, + Upper is Strict /\ 1. + +% negate_l(Strict,Lowerbound,X) +% +% Fails if X does not necessarily satisfy the lowerbound and strictness +% In other words: if adding the inverse of the lowerbound (X < L or X =< L) +% does not result in a failure, this predicate fails. + +negate_l(0,L,X) :- + { L > X }, + !, + fail. +negate_l(1,L,X) :- + { L > X }, + !, + fail. +negate_l(2,L,X) :- + { L >= X }, + !, + fail. +negate_l(3,L,X) :- + { L >= X }, + !, + fail. +negate_l(_,_,_). + +% negate_u(Strict,Upperbound,X) +% +% Fails if X does not necessarily satisfy the upperbound and strictness +% In other words: if adding the inverse of the upperbound (X > U or X >= U) +% does not result in a failure, this predicate fails. + +negate_u(0,U,X) :- + { U < X }, + !, + fail. +negate_u(1,U,X) :- + { U =< X }, + !, + fail. +negate_u(2,U,X) :- + { U < X }, + !, + fail. +negate_u(3,U,X) :- + { U =< X }, + !, + fail. +negate_u(_,_,_). + +% Profiling: these predicates are called during redundant and can be used +% to count the number of redundant bounds. + +red_t_l. +red_t_u. +red_t_L. +red_t_U. + diff --git a/LGPL/clpr/clpr/store.pl b/LGPL/clpr/clpr/store.pl new file mode 100644 index 000000000..eeab1a114 --- /dev/null +++ b/LGPL/clpr/clpr/store.pl @@ -0,0 +1,485 @@ +/* $Id: store.pl,v 1.1 2005-10-28 17:51:01 vsc Exp $ + + Part of CPL(R) (Constraint Logic Programming over Reals) + + Author: Leslie De Koninck + E-mail: Tom.Schrijvers@cs.kuleuven.ac.be + WWW: http://www.swi-prolog.org + http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 + Copyright (C): 2004, K.U. Leuven and + 1992-1995, Austrian Research Institute for + Artificial Intelligence (OFAI), + Vienna, Austria + + This software is part of Leslie De Koninck's master thesis, supervised + by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) + by Christian Holzbaur for SICStus Prolog and distributed under the + license details below with permission from all mentioned authors. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 + of the License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + As a special exception, if you link this library with other files, + compiled with a Free Software compiler, to produce an executable, this + library does not by itself cause the resulting executable to be covered + by the GNU General Public License. This exception does not however + invalidate any other reasons why the executable file might be covered by + the GNU General Public License. +*/ + +:- module(store, + [ + add_linear_11/3, + add_linear_f1/4, + add_linear_ff/5, + normalize_scalar/2, + delete_factor/4, + mult_linear_factor/3, + nf_rhs_x/4, + indep/2, + isolate/3, + nf_substitute/4, + mult_hom/3, + nf2sum/3, + nf_coeff_of/3, + renormalize/2 + ]). +:- use_module(arith_r, + [ + arith_normalize/2 + ]). + + +% normalize_scalar(S,[N,Z]) +% +% Transforms a scalar S into a linear expression [S,0] + +normalize_scalar(S,[N,Z]) :- + arith_normalize(S,N), + Z = 0.0. + +% renormalize(List,Lin) +% +% Renormalizes the not normalized linear expression in List into +% a normalized one. It does so to take care of unifications. +% (e.g. when a variable X is bound to a constant, the constant is added to +% the constant part of the linear expression; when a variable X is bound to +% another variable Y, the scalars of both are added) + +renormalize(List,Lin) :- + List = [I,R|Hom], + length(Hom,Len), + renormalize_log(Len,Hom,[],Lin0), + add_linear_11([I,R],Lin0,Lin). + +% renormalize_log(Len,Hom,HomTail,Lin) +% +% Logarithmically renormalizes the homogene part of a not normalized +% linear expression. See also renormalize/2. + +renormalize_log(1,[Term|Xs],Xs,Lin) :- + !, + Term = l(X*_,_), + renormalize_log_one(X,Term,Lin). +renormalize_log(2,[A,B|Xs],Xs,Lin) :- + !, + A = l(X*_,_), + B = l(Y*_,_), + renormalize_log_one(X,A,LinA), + renormalize_log_one(Y,B,LinB), + add_linear_11(LinA,LinB,Lin). +renormalize_log(N,L0,L2,Lin) :- + P is N>>1, + Q is N-P, + renormalize_log(P,L0,L1,Lp), + renormalize_log(Q,L1,L2,Lq), + add_linear_11(Lp,Lq,Lin). + +% renormalize_log_one(X,Term,Res) +% +% Renormalizes a term in X: if X is a nonvar, the term becomes a scalar. + +renormalize_log_one(X,Term,Res) :- + var(X), + Term = l(X*K,_), + get_attr(X,itf3,(_,_,_,order(OrdX),_)), % Order might have changed + Z = 0.0, + Res = [Z,Z,l(X*K,OrdX)]. +renormalize_log_one(X,Term,Res) :- + nonvar(X), + Term = l(X*K,_), + Xk is X*K, + normalize_scalar(Xk,Res). + +% ----------------------------- sparse vector stuff ---------------------------- % + +% add_linear_ff(LinA,Ka,LinB,Kb,LinC) +% +% Linear expression LinC is the result of the addition of the 2 linear expressions +% LinA and LinB, each one multiplied by a scalar (Ka for LinA and Kb for LinB). + +add_linear_ff(LinA,Ka,LinB,Kb,LinC) :- + LinA = [Ia,Ra|Ha], + LinB = [Ib,Rb|Hb], + LinC = [Ic,Rc|Hc], + Ic is Ia*Ka+Ib*Kb, + Rc is Ra*Ka+Rb*Kb, + add_linear_ffh(Ha,Ka,Hb,Kb,Hc). + +% add_linear_ffh(Ha,Ka,Hb,Kb,Hc) +% +% Homogene part Hc is the result of the addition of the 2 homogene parts Ha and Hb, +% each one multiplied by a scalar (Ka for Ha and Kb for Hb) + +add_linear_ffh([],_,Ys,Kb,Zs) :- mult_hom(Ys,Kb,Zs). +add_linear_ffh([l(X*Kx,OrdX)|Xs],Ka,Ys,Kb,Zs) :- + add_linear_ffh(Ys,X,Kx,OrdX,Xs,Zs,Ka,Kb). + +% add_linear_ffh(Ys,X,Kx,OrdX,Xs,Zs,Ka,Kb) +% +% Homogene part Zs is the result of the addition of the 2 homogene parts Ys and +% [l(X*Kx,OrdX)|Xs], each one multiplied by a scalar (Ka for [l(X*Kx,OrdX)|Xs] and Kb for Ys) + +add_linear_ffh([],X,Kx,OrdX,Xs,Zs,Ka,_) :- mult_hom([l(X*Kx,OrdX)|Xs],Ka,Zs). +add_linear_ffh([l(Y*Ky,OrdY)|Ys],X,Kx,OrdX,Xs,Zs,Ka,Kb) :- + compare(Rel,OrdX,OrdY), + ( + Rel = (=), + !, + Kz is Kx*Ka+Ky*Kb, + ( + TestKz is Kz - 0.0, % Kz =:= 0 + TestKz =< 1e-010, + TestKz >= -1e-010 -> + + !, + add_linear_ffh(Xs,Ka,Ys,Kb,Zs) + ; + Zs = [l(X*Kz,OrdX)|Ztail], + add_linear_ffh(Xs,Ka,Ys,Kb,Ztail) + ) + ; + Rel = (<), + !, + Zs = [l(X*Kz,OrdX)|Ztail], + Kz is Kx*Ka, + add_linear_ffh(Xs,Y,Ky,OrdY,Ys,Ztail,Kb,Ka) + ; + Rel = (>), + !, + Zs = [l(Y*Kz,OrdY)|Ztail], + Kz is Ky*Kb, + add_linear_ffh(Ys,X,Kx,OrdX,Xs,Ztail,Ka,Kb) + ). + +% add_linear_f1(LinA,Ka,LinB,LinC) +% +% special case of add_linear_ff with Kb = 1 + +add_linear_f1(LinA,Ka,LinB,LinC) :- + LinA = [Ia,Ra|Ha], + LinB = [Ib,Rb|Hb], + LinC = [Ic,Rc|Hc], + Ic is Ia*Ka+Ib, + Rc is Ra*Ka+Rb, + add_linear_f1h(Ha,Ka,Hb,Hc). + +% add_linear_f1h(Ha,Ka,Hb,Hc) +% +% special case of add_linear_ffh/5 with Kb = 1 + +add_linear_f1h([],_,Ys,Ys). +add_linear_f1h([l(X*Kx,OrdX)|Xs],Ka,Ys,Zs) :- + add_linear_f1h(Ys,X,Kx,OrdX,Xs,Zs,Ka). + +% add_linear_f1h(Ys,X,Kx,OrdX,Xs,Zs,Ka) +% +% special case of add_linear_ffh/8 with Kb = 1 + +add_linear_f1h([],X,Kx,OrdX,Xs,Zs,Ka) :- mult_hom([l(X*Kx,OrdX)|Xs],Ka,Zs). +add_linear_f1h([l(Y*Ky,OrdY)|Ys],X,Kx,OrdX,Xs,Zs,Ka) :- + compare(Rel,OrdX,OrdY), + ( + Rel = (=), + !, + Kz is Kx*Ka+Ky, + ( + TestKz is Kz - 0.0, % Kz =:= 0 + TestKz =< 1e-010, + TestKz >= -1e-010 -> + + !, + add_linear_f1h(Xs,Ka,Ys,Zs) + ; + Zs = [l(X*Kz,OrdX)|Ztail], + add_linear_f1h(Xs,Ka,Ys,Ztail) + ) + ; + Rel = (<), + !, + Zs = [l(X*Kz,OrdX)|Ztail], + Kz is Kx*Ka, + add_linear_f1h(Xs,Ka,[l(Y*Ky,OrdY)|Ys],Ztail) + ; + Rel = (>), + !, + Zs = [l(Y*Ky,OrdY)|Ztail], + add_linear_f1h(Ys,X,Kx,OrdX,Xs,Ztail,Ka) + ). + +% add_linear_11(LinA,LinB,LinC) +% +% special case of add_linear_ff with Ka = 1 and Kb = 1 + +add_linear_11(LinA,LinB,LinC) :- + LinA = [Ia,Ra|Ha], + LinB = [Ib,Rb|Hb], + LinC = [Ic,Rc|Hc], + Ic is Ia+Ib, + Rc is Ra+Rb, + add_linear_11h(Ha,Hb,Hc). + +% add_linear_11h(Ha,Hb,Hc) +% +% special case of add_linear_ffh/5 with Ka = 1 and Kb = 1 + +add_linear_11h([],Ys,Ys). +add_linear_11h([l(X*Kx,OrdX)|Xs],Ys,Zs) :- + add_linear_11h(Ys,X,Kx,OrdX,Xs,Zs). + +% add_linear_11h(Ys,X,Kx,OrdX,Xs,Zs) +% +% special case of add_linear_ffh/8 with Ka = 1 and Kb = 1 + +add_linear_11h([],X,Kx,OrdX,Xs,[l(X*Kx,OrdX)|Xs]). +add_linear_11h([l(Y*Ky,OrdY)|Ys],X,Kx,OrdX,Xs,Zs) :- + compare(Rel,OrdX,OrdY), + ( + Rel = (=), + !, + Kz is Kx+Ky, + ( + TestKz is Kz - 0.0, % Kz =:= 0 + TestKz =< 1e-010, + TestKz >= -1e-010 -> + + !, + add_linear_11h(Xs,Ys,Zs) + ; + Zs = [l(X*Kz,OrdX)|Ztail], + add_linear_11h(Xs,Ys,Ztail) + ) + ; + Rel = (<), + !, + Zs = [l(X*Kx,OrdX)|Ztail], + add_linear_11h(Xs,Y,Ky,OrdY,Ys,Ztail) + ; + Rel = (>), + !, + Zs = [l(Y*Ky,OrdY)|Ztail], + add_linear_11h(Ys,X,Kx,OrdX,Xs,Ztail) + ). + +% mult_linear_factor(Lin,K,Res) +% +% Linear expression Res is the result of multiplication of linear +% expression Lin by scalar K + +mult_linear_factor(Lin,K,Mult) :- + TestK is K - 1.0, % K =:= 1 + TestK =< 1e-010, + TestK >= -1e-010, % avoid copy + !, + Mult = Lin. +mult_linear_factor(Lin,K,Res) :- + Lin = [I,R|Hom], + Res = [Ik,Rk|Mult], + Ik is I*K, + Rk is R*K, + mult_hom(Hom,K,Mult). + +% mult_hom(Hom,K,Res) +% +% Homogene part Res is the result of multiplication of homogene part +% Hom by scalar K + +mult_hom([],_,[]). +mult_hom([l(A*Fa,OrdA)|As],F,[l(A*Fan,OrdA)|Afs]) :- + Fan is F*Fa, + mult_hom(As,F,Afs). + +% nf_substitute(Ord,Def,Lin,Res) +% +% Linear expression Res is the result of substitution of Var in +% linear expression Lin, by its definition in the form of linear +% expression Def + +nf_substitute(OrdV,LinV,LinX,LinX1) :- + delete_factor(OrdV,LinX,LinW,K), + add_linear_f1(LinV,K,LinW,LinX1). + +% delete_factor(Ord,Lin,Res,Coeff) +% +% Linear expression Res is the result of the deletion of the term +% Var*Coeff where Var has ordering Ord from linear expression Lin + +delete_factor(OrdV,Lin,Res,Coeff) :- + Lin = [I,R|Hom], + Res = [I,R|Hdel], + delete_factor_hom(OrdV,Hom,Hdel,Coeff). + +% delete_factor_hom(Ord,Hom,Res,Coeff) +% +% Homogene part Res is the result of the deletion of the term +% Var*Coeff from homogene part Hom + +delete_factor_hom(VOrd,[Car|Cdr],RCdr,RKoeff) :- + Car = l(_*Koeff,Ord), + compare(Rel,VOrd,Ord), + ( + Rel= (=), + !, + RCdr = Cdr, + RKoeff=Koeff + ; + Rel= (>), + !, + RCdr = [Car|RCdr1], + delete_factor_hom(VOrd,Cdr,RCdr1,RKoeff) + ). + + +% nf_coeff_of(Lin,OrdX,Coeff) +% +% Linear expression Lin contains the term l(X*Coeff,OrdX) + +nf_coeff_of(Lin,VOrd,Coeff) :- + Lin = [_,_|Hom], + nf_coeff_hom(Hom,VOrd,Coeff), + !. + +% nf_coeff_hom(Lin,OrdX,Coeff) +% +% Linear expression Lin contains the term l(X*Coeff,OrdX) where the +% order attribute of X = OrdX + +nf_coeff_hom([l(_*K,OVar)|Vs],OVid,Coeff) :- + compare(Rel,OVid,OVar), + ( + Rel = (=), + !, + Coeff = K + ; + Rel = (>), + !, + nf_coeff_hom(Vs,OVid,Coeff) + ). + +% nf_rhs_x(Lin,OrdX,Rhs,K) +% +% Rhs = R + I where Lin = [I,R|Hom] and l(X*K,OrdX) is a term of Hom + +nf_rhs_x(Lin,OrdX,Rhs,K) :- + Lin = [I,R|Tail], + nf_coeff_hom(Tail,OrdX,K), + Rhs is R+I. % late because X may not occur in H + +% isolate(OrdN,Lin,Lin1) +% +% Linear expression Lin1 is the result of the transformation of linear expression +% Lin = 0 which contains the term l(New*K,OrdN) into an equivalent expression Lin1 = New. + +isolate(OrdN,Lin,Lin1) :- + delete_factor(OrdN,Lin,Lin0,Coeff), + K is -1.0/Coeff, + mult_linear_factor(Lin0,K,Lin1). + +% indep(Lin,OrdX) +% +% succeeds if Lin = [0,_|[l(X*1,OrdX)]] + +indep(Lin,OrdX) :- + Lin = [I,_|[l(_*K,OrdY)]], + OrdX == OrdY, + TestK is K - 1.0, % K =:= 1 + TestK =< 1e-010, + TestK >= -1e-010, + TestI is I - 0.0, % I =:= 0 + TestI =< 1e-010, + TestI >= -1e-010. + +% nf2sum(Lin,Sofar,Term) +% +% Transforms a linear expression into a sum +% (e.g. the expression [5,_,[l(X*2,OrdX),l(Y*-1,OrdY)]] gets transformed into 5 + 2*X - Y) + +nf2sum([],I,I). +nf2sum([X|Xs],I,Sum) :- + ( + TestI is I - 0.0, % I =:= 0 + TestI =< 1e-010, + TestI >= -1e-010 -> + + X = l(Var*K,_), + ( + TestK is K - 1.0, % K =:= 1 + TestK =< 1e-010, + TestK >= -1e-010 -> + + hom2sum(Xs,Var,Sum) + ; + TestK is K - -1.0, % K =:= -1 + TestK =< 1e-010, + TestK >= -1e-010 -> + + hom2sum(Xs,-Var,Sum) + ; + hom2sum(Xs,K*Var,Sum) + ) + ; + hom2sum([X|Xs],I,Sum) + ). + +% hom2sum(Hom,Sofar,Term) +% +% Transforms a linear expression into a sum +% this predicate handles all but the first term +% (the first term does not need a concatenation symbol + or -) +% see also nf2sum/3 + +hom2sum([],Term,Term). +hom2sum([l(Var*K,_)|Cs],Sofar,Term) :- + ( + TestK is K - 1.0, % K =:= 1 + TestK =< 1e-010, + TestK >= -1e-010 -> + + Next = Sofar + Var + ; + TestK is K - -1.0, % K =:= -1 + TestK =< 1e-010, + TestK >= -1e-010 -> + + Next = Sofar - Var + + ; + K - 0.0 < -1e-010 -> + + Ka is -K, + Next = Sofar - Ka*Var + ; + Next = Sofar + K*Var + ), + hom2sum(Cs,Next,Term). diff --git a/LGPL/clpr/clpr/ugraphs.pl b/LGPL/clpr/clpr/ugraphs.pl new file mode 100644 index 000000000..85f82d23a --- /dev/null +++ b/LGPL/clpr/clpr/ugraphs.pl @@ -0,0 +1,828 @@ +/* Copyright(C) 1992, Swedish Institute of Computer Science */ + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% File : UGRAPHS.PL % +% Maintainer : Mats Carlsson % +% New versions of transpose/2, reduce/2, top_sort/2 by Dan Sahlin % +% Updated: 3 September 1999 % +% Purpose: Unweighted graph-processing utilities % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +/* Adapted from shared code written by Richard A O'Keefe */ + +/* An unweighted directed graph (ugraph) is represented as a list of + (vertex-neighbors) pairs, where the pairs are in standard order + (as produced by keysort with unique keys) and the neighbors of + each vertex are also in standard order (as produced by sort), and + every neighbor appears as a vertex even if it has no neighbors + itself. + + An undirected graph is represented as a directed graph where for + each edge (U,V) there is a symmetric edge (V,U). + + An edge (U,V) is represented as the term U-V. + + A vertex can be any term. Two vertices are distinct iff they are + not identical (==). + + A path is represented as a list of vertices. + No vertex can appear twice in a path. +*/ + +:- module(ugraphs, [ + vertices_edges_to_ugraph/3, + vertices/2, + edges/2, + add_vertices/3, + del_vertices/3, + add_edges/3, + del_edges/3, + transpose/2, + neighbors/3, + neighbours/3, + complement/2, + compose/3, + transitive_closure/2, + symmetric_closure/2, + top_sort/2, + max_path/5, + min_path/5, + min_paths/3, + path/3, + reduce/2, + reachable/3, + random_ugraph/3, + min_tree/3, + clique/3, + independent_set/3, + coloring/3, + colouring/3 + ]). +:- use_module(library(ordsets), + [ ord_del_element/3, + ord_add_element/3, + ord_subset/2, + ord_union/3, + ord_union/4, + ord_disjoint/2, + ord_intersection/3 + ]). + + +:- use_module(library(lists), [ + append/3, + member/2, + reverse/2 + ]). + +:- use_module(library(assoc), [ + list_to_assoc/2, + ord_list_to_assoc/2, + get_assoc/3, + get_assoc/5 + ]). + +/*:- use_module(random, [ + random/1 + ]). +*/ +% vertices_edges_to_ugraph(+Vertices, +Edges, -Graph) +% is true if Vertices is a list of vertices, Edges is a list of edges, +% and Graph is a graph built from Vertices and Edges. Vertices and +% Edges may be in any order. The vertices mentioned in Edges do not +% have to occur explicitly in Vertices. Vertices may be used to +% specify vertices that are not connected to any edges. + +vertices_edges_to_ugraph(Vertices0, Edges, Graph) :- + sort(Vertices0, Vertices1), + sort(Edges, EdgeSet), + edges_vertices(EdgeSet, Bag), + sort(Bag, Vertices2), + ord_union(Vertices1, Vertices2, VertexSet), + group_edges(VertexSet, EdgeSet, Graph). + +edges_vertices([], []). +edges_vertices([From-To|Edges], [From,To|Vertices]) :- + edges_vertices(Edges, Vertices). + +group_edges([], _, []). +group_edges([Vertex|Vertices], Edges, [Vertex-Neibs|G]) :- + group_edges(Edges, Vertex, Neibs, RestEdges), + group_edges(Vertices, RestEdges, G). + +group_edges([V0-X|Edges], V, [X|Neibs], RestEdges) :- V0==V, !, + group_edges(Edges, V, Neibs, RestEdges). +group_edges(Edges, _, [], Edges). + + + +% vertices(+Graph, -Vertices) +% unifies Vertices with the vertices in Graph. + +vertices([], []). +vertices([Vertex-_|Graph], [Vertex|Vertices]) :- vertices(Graph, Vertices). + + + +% edges(+Graph, -Edges) +% unifies Edges with the edges in Graph. + +edges([], []). +edges([Vertex-Neibs|G], Edges) :- + edges(Neibs, Vertex, Edges, MoreEdges), + edges(G, MoreEdges). + +edges([], _, Edges, Edges). +edges([Neib|Neibs], Vertex, [Vertex-Neib|Edges], MoreEdges) :- + edges(Neibs, Vertex, Edges, MoreEdges). + + + +% add_vertices(+Graph1, +Vertices, -Graph2) +% is true if Graph2 is Graph1 with Vertices added to it. + +add_vertices(Graph0, Vs0, Graph) :- + sort(Vs0, Vs), + vertex_units(Vs, Graph1), + graph_union(Graph0, Graph1, Graph). + + + +% del_vertices(+Graph1, +Vertices, -Graph2) +% is true if Graph2 is Graph1 with Vertices and all edges to and from +% Vertices removed from it. + +del_vertices(Graph0, Vs0, Graph) :- + sort(Vs0, Vs), + graph_del_vertices(Graph0, Vs, Vs, Graph). + + + +% add_edges(+Graph1, +Edges, -Graph2) +% is true if Graph2 is Graph1 with Edges and their "to" and "from" +% vertices added to it. + +add_edges(Graph0, Edges0, Graph) :- + sort(Edges0, EdgeSet), + edges_vertices(EdgeSet, Vs0), + sort(Vs0, Vs), + group_edges(Vs, EdgeSet, Graph1), + graph_union(Graph0, Graph1, Graph). + + + +% del_edges(+Graph1, +Edges, -Graph2) +% is true if Graph2 is Graph1 with Edges removed from it. + +del_edges(Graph0, Edges0, Graph) :- + sort(Edges0, EdgeSet), + edges_vertices(EdgeSet, Vs0), + sort(Vs0, Vs), + group_edges(Vs, EdgeSet, Graph1), + graph_difference(Graph0, Graph1, Graph). + + + +vertex_units([], []). +vertex_units([V|Vs], [V-[]|Us]) :- vertex_units(Vs, Us). + + +graph_union(G0, [], G) :- !, G = G0. +graph_union([], G, G). +graph_union([V1-N1|G1], [V2-N2|G2], G) :- + compare(C, V1, V2), + graph_union(C, V1, N1, G1, V2, N2, G2, G). + +graph_union(<, V1, N1, G1, V2, N2, G2, [V1-N1|G]) :- + graph_union(G1, [V2-N2|G2], G). +graph_union(=, V, N1, G1, _, N2, G2, [V-N|G]) :- + ord_union(N1, N2, N), + graph_union(G1, G2, G). +graph_union(>, V1, N1, G1, V2, N2, G2, [V2-N2|G]) :- + graph_union([V1-N1|G1], G2, G). + + + +graph_difference(G0, [], G) :- !, G = G0. +graph_difference([], _, []). +graph_difference([V1-N1|G1], [V2-N2|G2], G) :- + compare(C, V1, V2), + graph_difference(C, V1, N1, G1, V2, N2, G2, G). + +graph_difference(<, V1, N1, G1, V2, N2, G2, [V1-N1|G]) :- + graph_difference(G1, [V2-N2|G2], G). +graph_difference(=, V, N1, G1, _, N2, G2, [V-N|G]) :- + ord_subtract(N1, N2, N), + graph_difference(G1, G2, G). +graph_difference(>, V1, N1, G1, _, _, G2, G) :- + graph_difference([V1-N1|G1], G2, G). + + +graph_del_vertices(G1, [], Set, G) :- !, + graph_del_vertices(G1, Set, G). +graph_del_vertices([], _, _, []). +graph_del_vertices([V1-N1|G1], [V2|Vs], Set, G) :- + compare(C, V1, V2), + graph_del_vertices(C, V1, N1, G1, V2, Vs, Set, G). + +graph_del_vertices(<, V1, N1, G1, V2, Vs, Set, [V1-N|G]) :- + ord_subtract(N1, Set, N), + graph_del_vertices(G1, [V2|Vs], Set, G). +graph_del_vertices(=, _, _, G1, _, Vs, Set, G) :- + graph_del_vertices(G1, Vs, Set, G). +graph_del_vertices(>, V1, N1, G1, _, Vs, Set, G) :- + graph_del_vertices([V1-N1|G1], Vs, Set, G). + +graph_del_vertices([], _, []). +graph_del_vertices([V1-N1|G1], Set, [V1-N|G]) :- + ord_subtract(N1, Set, N), + graph_del_vertices(G1, Set, G). + + + +% transpose(+Graph, -Transpose) +% is true if Transpose is the graph computed by replacing each edge +% (u,v) in Graph by its symmetric edge (v,u). It can only be used +% one way around. The cost is O(N log N). + +transpose(Graph, Transpose) :- + transpose_edges(Graph, TEdges, []), + sort(TEdges, TEdges2), + vertices(Graph, Vertices), + group_edges(Vertices, TEdges2, Transpose). + +transpose_edges([]) --> []. +transpose_edges([Vertex-Neibs|G]) --> + transpose_edges(Neibs, Vertex), + transpose_edges(G). + +transpose_edges([], _) --> []. +transpose_edges([Neib|Neibs], Vertex) --> [Neib-Vertex], + transpose_edges(Neibs, Vertex). + + +% neighbours(+Vertex, +Graph, -Neighbors) +% neighbors(+Vertex, +Graph, -Neighbors) +% is true if Vertex is a vertex in Graph and Neighbors are its neighbors. + +neighbours(Vertex, Graph, Neighbors) :- + neighbors(Vertex, Graph, Neighbors). + +neighbors(V, [V0-Neighbors|_], Neighbors) :- V0==V, !. +neighbors(V, [_|Graph], Neighbors) :- neighbors(V, Graph, Neighbors). + + + +% complement(+Graph, -Complement) +% Complement is the complement graph of Graph, i.e. the graph that has +% the same vertices as Graph but only the edges that are not in Graph. + +complement(Graph, Complement) :- + vertices(Graph, Vertices), + complement(Graph, Vertices, Complement). + +complement([], _, []). +complement([V-Neibs|Graph], Vertices, [V-Others|Complement]) :- + ord_add_element(Neibs, V, Neibs1), + ord_subtract(Vertices, Neibs1, Others), + complement(Graph, Vertices, Complement). + + + +% compose(+G1, +G2, -Composition) +% computes Composition as the composition of two graphs, which need +% not have the same set of vertices. + +compose(G1, G2, Composition) :- + vertices(G1, V1), + vertices(G2, V2), + ord_union(V1, V2, V), + compose(V, G1, G2, Composition). + +compose([], _, _, []). +compose([V0|Vertices], [V-Neibs|G1], G2, [V-Comp|Composition]) :- V==V0, !, + compose1(Neibs, G2, [], Comp), + compose(Vertices, G1, G2, Composition). +compose([V|Vertices], G1, G2, [V-[]|Composition]) :- + compose(Vertices, G1, G2, Composition). + +compose1([V1|Vs1], [V2-N2|G2], SoFar, Comp) :- !, + compare(Rel, V1, V2), + compose1(Rel, V1, Vs1, V2, N2, G2, SoFar, Comp). +compose1(_, _, Comp, Comp). + +compose1(<, _, Vs1, V2, N2, G2, SoFar, Comp) :- + compose1(Vs1, [V2-N2|G2], SoFar, Comp). +compose1(>, V1, Vs1, _, _, G2, SoFar, Comp) :- + compose1([V1|Vs1], G2, SoFar, Comp). +compose1(=, V1, Vs1, V1, N2, G2, SoFar, Comp) :- + ord_union(N2, SoFar, Next), + compose1(Vs1, G2, Next, Comp). + + + +% transitive_closure(+Graph, -Closure) +% computes Closure as the transitive closure of Graph in O(N^3) time. + +transitive_closure(Graph, Closure) :- + warshall(Graph, Graph, Closure). + +warshall([], Closure, Closure). +warshall([V-_|G], E, Closure) :- + neighbors(V, E, Y), + warshall(E, V, Y, NewE), + warshall(G, NewE, Closure). + +warshall([], _, _, []). +warshall([X-Neibs|G], V, Y, [X-NewNeibs|NewG]) :- + ord_subset([V], Neibs), !, + ord_del_element(Y, X, Y1), + ord_union(Neibs, Y1, NewNeibs), + warshall(G, V, Y, NewG). +warshall([X-Neibs|G], V, Y, [X-Neibs|NewG]) :- + warshall(G, V, Y, NewG). + + + +% symmetric_closure(+Graph, -Closure) +% computes Closure as the symmetric closure of Graph, i.e. for each +% edge (u,v) in Graph, add its symmetric edge (v,u). Approx O(N log N) +% time. This is useful for making a directed graph undirected. + +symmetric_closure(Graph, Closure) :- + transpose(Graph, Transpose), + symmetric_closure(Graph, Transpose, Closure). + +symmetric_closure([], [], []). +symmetric_closure([V-Neibs1|Graph], [V-Neibs2|Transpose], [V-Neibs|Closure]) :- + ord_union(Neibs1, Neibs2, Neibs), + symmetric_closure(Graph, Transpose, Closure). + + + +% top_sort(+Graph, -Sorted) +% finds a topological ordering of a Graph and returns the ordering +% as a list of Sorted vertices. Fails iff no ordering exists, i.e. +% iff the graph contains cycles. Approx O(N log N) time. + +top_sort(Graph, Sorted) :- + fanin_counts(Graph, Counts), + get_top_elements(Counts, Top, 0, I), + ord_list_to_assoc(Counts, Map), + top_sort(Top, I, Map, Sorted). + +top_sort([], 0, _, []). +top_sort([V-VN|Top0], I, Map0, [V|Sorted]) :- + dec_counts(VN, I, J, Map0, Map, Top0, Top), + top_sort(Top, J, Map, Sorted). + +dec_counts([], I, I, Map, Map, Top, Top). +dec_counts([N|Ns], I, K, Map0, Map, Top0, Top) :- + get_assoc(N, Map0, NN-C0, Map1, NN-C), + C is C0-1, + (C=:=0 -> J is I-1, Top1 = [N-NN|Top0]; J = I, Top1 = Top0), + dec_counts(Ns, J, K, Map1, Map, Top1, Top). + +get_top_elements([], [], I, I). +get_top_elements([V-(VN-C)|Counts], Top0, I, K) :- + (C=:=0 -> J = I, Top0 = [V-VN|Top1]; J is I+1, Top0 = Top1), + get_top_elements(Counts, Top1, J, K). + +fanin_counts(Graph, Counts) :- + transpose_edges(Graph, Edges0, []), + keysort(Edges0, Edges), + fanin_counts(Graph, Edges, Counts). + +fanin_counts([], [], []). +fanin_counts([V-VN|Graph], Edges, [V-(VN-C)|Counts]) :- + fanin_counts(Edges, V, 0, C, Edges1), + fanin_counts(Graph, Edges1, Counts). + +fanin_counts([V-_|Edges0], V0, C0, C, Edges) :- + V==V0, !, + C1 is C0+1, + fanin_counts(Edges0, V0, C1, C, Edges). +fanin_counts(Edges, _, C, C, Edges). + + +% max_path(+V1, +V2, +Graph, -Path, -Cost) +% is true if Path is a list of vertices constituting a longest path +% of cost Cost from V1 to V2 in Graph, there being no cyclic paths from +% V1 to V2. Takes O(N^2) time. + +max_path(Initial, Final, Graph, Path, Cost) :- + transpose(Graph, TGraph), + max_path_init(Initial, Final, Graph, TGraph, TGraph2, Order), + max_path_init(TGraph2, Val0), + max_path(Order, TGraph2, Val0, Val), + max_path_select(Val, Path, Cost). + +max_path_init(Initial, Final, Graph, TGraph, TGraph2, Order) :- + reachable(Initial, Graph, InitialReachable), + reachable(Final, TGraph, FinalReachable), + ord_intersection(InitialReachable, FinalReachable, Reachable), + subgraph(TGraph, Reachable, TGraph2), + top_sort(TGraph2, Order). + +max_path_init([], []). +max_path_init([V-_|G], [V-([]-0)|Val]) :- max_path_init(G, Val). + +max_path_select([V-(Longest-Max)|Val], Path, Cost) :- + max_path_select(Val, V, Longest, Path, Max, Cost). + +max_path_select([], V, Path, [V|Path], Cost, Cost). +max_path_select([V1-(Path1-Cost1)|Val], V2, Path2, Path, Cost2, Cost) :- + ( Cost1>Cost2 -> + max_path_select(Val, V1, Path1, Path, Cost1, Cost) + ; max_path_select(Val, V2, Path2, Path, Cost2, Cost) + ). + +max_path([], _, Val, Val). +max_path([V|Order], Graph, Val0, Val) :- + neighbors(V, Graph, Neibs), + neighbors(V, Val0, Item), + max_path_update(Neibs, V-Item, Val0, Val1), + max_path(Order, Graph, Val1, Val). + +%% [MC] 3.8.6: made determinate +max_path_update([], _, Val, Val). +max_path_update([N|Neibs], Item, [Item0|Val0], Val) :- + Item0 = V0-(_-Cost0), + N==V0, !, + Item = V-(Path-Cost), + Cost1 is Cost+1, + ( Cost1>Cost0 -> Val = [V0-([V|Path]-Cost1)|Val1] + ; Val = [Item0|Val1] + ), + max_path_update(Neibs, Item, Val0, Val1). +max_path_update(Neibs, Item, [X|Val0], [X|Val]) :- + Neibs = [_|_], + max_path_update(Neibs, Item, Val0, Val). + +subgraph([], _, []). +subgraph([V-Neibs|Graph], Vs, [V-Neibs1|Subgraph]) :- + ord_subset([V], Vs), !, + ord_intersection(Neibs, Vs, Neibs1), + subgraph(Graph, Vs, Subgraph). +subgraph([_|Graph], Vs, Subgraph) :- + subgraph(Graph, Vs, Subgraph). + + + +% min_path(+V1, +V2, +Graph, -Path, -Length) +% is true if Path is a list of vertices constituting a shortest path +% of length Length from V1 to V2 in Graph. Takes O(N^2) time. + +min_path(Initial, Final, Graph, Path, Length) :- + min_path([[Initial]|Q], Q, [Initial], Final, Graph, Rev), + reverse(Rev, Path), + length(Path, N), + Length is N-1. + +min_path(Head0, Tail0, Closed0, Final, Graph, Rev) :- + Head0 \== Tail0, + Head0 = [Sofar|Head], + Sofar = [V|_], + ( V==Final -> Rev = Sofar + ; neighbors(V, Graph, Neibs), + ord_union(Closed0, Neibs, Closed, Neibs1), + enqueue(Neibs1, Sofar, Tail0, Tail), + min_path(Head, Tail, Closed, Final, Graph, Rev) + ). + +enqueue([], _) --> []. +enqueue([V|Vs], Sofar) --> [[V|Sofar]], enqueue(Vs, Sofar). + + + +% min_paths(+Vertex, +Graph, -Tree) +% is true if Tree is a tree of all the shortest paths from Vertex to +% every other vertex in Graph. This is the single-source shortest +% paths problem. The algorithm is straightforward. + +min_paths(Vertex, Graph, Tree) :- + min_paths([Vertex], Graph, [Vertex], List), + keysort(List, Tree). + +min_paths([], _, _, []). +min_paths([Q|R], Graph, Reach0, [Q-New|List]) :- + neighbors(Q, Graph, Neibs), + ord_union(Reach0, Neibs, Reach, New), + append(R, New, S), + min_paths(S, Graph, Reach, List). + + + +% path(+Vertex, +Graph, -Path) +% is given a Graph and a Vertex of that Graph, and returns a maximal +% Path rooted at Vertex, enumerating more Paths on backtracking. + +path(Initial, Graph, Path) :- + path([Initial], [], Graph, Path). + +path(Q, Not, Graph, Path) :- + Q = [Qhead|_], + neighbors(Qhead, Graph, Neibs), + ord_subtract(Neibs, Not, Neibs1), + ( Neibs1 = [] -> reverse(Q, Path) + ; ord_add_element(Not, Qhead, Not1), + member(N, Neibs1), + path([N|Q], Not1, Graph, Path) + ). + + + +% reduce(+Graph, -Reduced) +% is true if Reduced is the reduced graph for Graph. The vertices of +% the reduced graph are the strongly connected components of Graph. +% There is an edge in Reduced from u to v iff there is an edge in +% Graph from one of the vertices in u to one of the vertices in v. A +% strongly connected component is a maximal set of vertices where +% each vertex has a path to every other vertex. +% Algorithm from "Algorithms" by Sedgewick, page 482, Tarjan's algorithm. +% Approximately linear in the maximum of arcs and nodes (O(N log N)). + +reduce(Graph, Reduced) :- + strong_components(Graph, SCCS, Map), + reduced_vertices_edges(Graph, Vertices, Map, Edges, []), + sort(Vertices, Vertices1), + sort(Edges, Edges1), + group_edges(Vertices1, Edges1, Reduced), + sort(SCCS, Vertices1). + +strong_components(Graph, SCCS, A) :- + nodeinfo(Graph, Nodeinfo, Vertices), + ord_list_to_assoc(Nodeinfo, A0), + visit(Vertices, 0, _, A0, A, 0, _, [], _, SCCS, []). + +visit([], Min, Min, A, A, I, I, Stk, Stk) --> []. +visit([V|Vs], Min0, Min, A0, A, I, M, Stk0, Stk) --> + {get_assoc(V, A0, node(Ns,J,Eq), A1, node(Ns,K,Eq))}, + ( {J>0} -> + {J=K, J=Min1, A1=A3, I=L, Stk0=Stk2} + ; {K is I+1}, + visit(Ns, K, Min1, A1, A2, K, L, [V|Stk0], Stk1), + ( {K>Min1} -> {A2=A3, Stk1=Stk2} + ; pop(V, Eq, A2, A3, Stk1, Stk2, []) + ) + ), + {Min2 is min(Min0,Min1)}, + visit(Vs, Min2, Min, A3, A, L, M, Stk2, Stk). + +pop(V, Eq, A0, A, [V1|Stk0], Stk, SCC0) --> + {get_assoc(V1, A0, node(Ns,_,Eq), A1, node(Ns,16'100000,Eq))}, + ( {V==V1} -> [SCC], {A1=A, Stk0=Stk, sort([V1|SCC0], SCC)} + ; pop(V, Eq, A1, A, Stk0, Stk, [V1|SCC0]) + ). + +nodeinfo([], [], []). +nodeinfo([V-Ns|G], [V-node(Ns,0,_)|Nodeinfo], [V|Vs]) :- + nodeinfo(G, Nodeinfo, Vs). + +reduced_vertices_edges([], [], _) --> []. +reduced_vertices_edges([V-Neibs|Graph], [V1|Vs], Map) --> + {get_assoc(V, Map, N), N=node(_,_,V1)}, + reduced_edges(Neibs, V1, Map), + reduced_vertices_edges(Graph, Vs, Map). + +reduced_edges([], _, _) --> []. +reduced_edges([V|Vs], V1, Map) --> + {get_assoc(V, Map, N), N=node(_,_,V2)}, + ({V1==V2} -> []; [V1-V2]), + reduced_edges(Vs, V1, Map). + + +% reachable(+Vertex, +Graph, -Reachable) +% is given a Graph and a Vertex of that Graph, and returns the set +% of vertices that are Reachable from that Vertex. Takes O(N^2) +% time. + +reachable(Initial, Graph, Reachable) :- + reachable([Initial], Graph, [Initial], Reachable). + +reachable([], _, Reachable, Reachable). +reachable([Q|R], Graph, Reach0, Reachable) :- + neighbors(Q, Graph, Neighbors), + ord_union(Reach0, Neighbors, Reach1, New), + append(R, New, S), + reachable(S, Graph, Reach1, Reachable). + + + +% random_ugraph(+P, +N, -Graph) +% where P is a probability, unifies Graph with a random graph of N +% vertices where each possible edge is included with probability P. + +random_ugraph(P, N, Graph) :- + ( float(P), P >= 0.0, P =< 1.0 -> true + ; prolog:illarg(domain(float,between(0.0,1.0)), + random_ugraph(P,N,Graph), 1) + ), + ( integer(N), N >= 0 -> true + ; prolog:illarg(domain(integer,>=(0)), + random_ugraph(P,N,Graph), 2) + ), + random_ugraph(0, N, P, Graph). + +random_ugraph(N, N, _, Graph) :- !, Graph = []. +random_ugraph(I, N, P, [J-List|Graph]) :- + J is I+1, + random_neighbors(N, J, P, List, []), + random_ugraph(J, N, P, Graph). + +random_neighbors(0, _, _, S0, S) :- !, S = S0. +random_neighbors(N, J, P, S0, S) :- + ( N==J -> S1 = S + ; random(X), X > P -> S1 = S + ; S1 = [N|S] + ), + M is N-1, + random_neighbors(M, J, P, S0, S1). + +random(X) :- % JW: was undefined. Assume + X is random(10000)/10000. % we need 0<=X<=1 + + +% min_tree(+Graph, -Tree, -Cost) +% is true if Tree is a spanning tree of an *undirected* Graph with +% cost Cost, if it exists. Using a version of Prim's algorithm. + +min_tree([V-Neibs|Graph], Tree, Cost) :- + length(Graph, Cost), + prim(Cost, Neibs, Graph, [V], Edges), + vertices_edges_to_ugraph([], Edges, Tree). + +%% [MC] 3.8.6: made determinate +prim(0, [], _, _, []) :- !. +prim(I, [V|Vs], Graph, Reach, [V-W,W-V|Edges]) :- + neighbors(V, Graph, Neibs), + ord_subtract(Neibs, Reach, Neibs1), + ord_subtract(Neibs, Neibs1, [W|_]), + ord_add_element(Reach, V, Reach1), + ord_union(Vs, Neibs1, Vs1), + J is I-1, + prim(J, Vs1, Graph, Reach1, Edges). + + + +% clique(+Graph, +K, -Clique) +% is true if Clique is a maximal clique (complete subgraph) of N +% vertices of an *undirected* Graph, where N>=K. Adapted from +% Algorithm 457, "Finding All Cliques of an Undirected Graph [H]", +% Version 1, by Coen Bron and Joep Kerbosch, CACM vol. 6 no. 9 pp. +% 575-577, Sept. 1973. + +clique(Graph, K, Clique) :- + ( integer(K), K >= 0 -> true + ; prolog:illarg(domain(integer,>=(0)), + clique(Graph,K,Clique), 2) + ), + J is K-1, + prune(Graph, [], J, Graph1), + clique1(Graph1, J, Clique). + +clique1([], J, []) :- J < 0. +clique1([C-Neibs|Graph], J, [C|Clique]) :- + neighbor_graph(Graph, Neibs, C, Vs, Graph1), + J1 is J-1, + prune(Graph1, Vs, J1, Graph2), + clique1(Graph2, J1, Clique). +clique1([C-Neibs|Graph], J, Clique) :- + prune(Graph, [C], J, Graph2), + clique1(Graph2, J, Clique), + \+ ord_subset(Clique, Neibs). + +neighbor_graph([], _, _, [], []). +neighbor_graph([V0-Neibs0|Graph0], [V|Vs], W, Del, [V-Neibs|Graph]) :- + V0==V, !, + ord_del_element(Neibs0, W, Neibs), + neighbor_graph(Graph0, Vs, W, Del, Graph). +neighbor_graph([V-_|Graph0], Vs, W, [V|Del], Graph) :- + neighbor_graph(Graph0, Vs, W, Del, Graph). + +prune(Graph0, [], K, Graph) :- K =< 0, !, Graph = Graph0. +prune(Graph0, Vs0, K, Graph) :- + prune(Graph0, Vs0, K, Graph1, Vs1), + ( Vs1==[] -> Graph = Graph1 + ; prune(Graph1, Vs1, K, Graph) + ). + +prune([], _, _, [], []). +prune([V-Ns0|Graph0], Vs1, K, [V-Ns|Graph], Vs2) :- + ord_disjoint([V], Vs1), + ord_subtract(Ns0, Vs1, Ns), + length(Ns, I), + I >= K, !, + prune(Graph0, Vs1, K, Graph, Vs2). +prune([V-_|Graph0], Vs1, K, Graph, Vs2) :- + ( ord_disjoint([V], Vs1) -> Vs2 = [V|Vs3] + ; Vs2 = Vs3 + ), + prune(Graph0, Vs1, K, Graph, Vs3). + + + +% independent_set(+Graph, +K, -Set) +% is true if Set is a maximal independent (unconnected) set of N +% vertices of an *undirected* Graph, where N>=K. + +independent_set(Graph, K, Set) :- + complement(Graph, Complement), + clique(Complement, K, Set). + + + +% colouring(+Graph, +K, -Coloring) +% coloring(+Graph, +K, -Coloring) +% is true if Coloring is a mapping from vertices to colors 1..N of +% an *undirected* Graph such that all edges have distinct end colors, +% where N== 0 -> true + ; prolog:illarg(domain(integer,>=(0)), + coloring(Graph,K,Coloring), 2) + ), + color_map(Graph, Coloring), + color_map(Graph, Graph1, Coloring, Coloring), + coloring(Graph1, K, 0, [], Stack), + color_stack(Stack). + +coloring([], _, _, Stk0, Stk) :- !, Stk0 = Stk. +coloring(Graph, K, InUse, Stk0, Stk) :- + select_vertex(Graph, K, Compare, -, 0+0, V-Ns), + graph_del_vertices(Graph, [V], [V], Graph1), + ( Compare = < -> + coloring(Graph1, K, InUse, [V-Ns|Stk0], Stk) + ; M is min(K,InUse+1), + vertex_color(Ns, 1, M, V), + add_color(Graph1, Ns, V, Graph2), + InUse1 is max(V,InUse), + coloring(Graph2, K, InUse1, Stk0, Stk) + ). + +% select_vertex(+Graph, +K, -Comp, +Pair0, +Rank, -Pair) +% return any vertex with degree=), break ties using +% max. degree + +select_vertex([], _, >=, Pair, _, Pair). +select_vertex([V-Neibs|Graph], K, Comp, Pair0, Rank0, Pair) :- + evaluate_vertex(Neibs, 0, Rank), + ( Rank < K -> Comp = <, Pair = V-Neibs + ; Rank @> Rank0 -> + select_vertex(Graph, K, Comp, V-Neibs, Rank, Pair) + ; select_vertex(Graph, K, Comp, Pair0, Rank0, Pair) + ). + +evaluate_vertex([V|Neibs], Deg, Rank) :- var(V), !, + Deg1 is Deg+1, + evaluate_vertex(Neibs, Deg1, Rank). +evaluate_vertex(Neibs, Deg, Sat+Deg) :- + prolog:length(Neibs, 0, Sat). + +add_color([], _, _, []). +add_color([V-Neibs|Graph], [W|Vs], C, [V-Neibs1|Graph1]) :- V==W, !, + ord_add_element(Neibs, C, Neibs1), + add_color(Graph, Vs, C, Graph1). +add_color([Pair|Graph], Vs, C, [Pair|Graph1]) :- + add_color(Graph, Vs, C, Graph1). + +vertex_color([V|Vs], I, M, Color) :- V@, Ns, [_|Coloring], Ns1) :- + map_colors(Ns, Coloring, Ns1).