1/* $Id$ 2 3 Part of CLP(R) (Constraint Logic Programming over Reals) 4 5 Author: Leslie De Koninck 6 E-mail: Leslie.DeKoninck@cs.kuleuven.be 7 WWW: http://www.swi-prolog.org 8 http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 9 Copyright (C): 2004, K.U. Leuven and 10 1992-1995, Austrian Research Institute for 11 Artificial Intelligence (OFAI), 12 Vienna, Austria 13 14 This software is part of Leslie De Koninck's master thesis, supervised 15 by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) 16 by Christian Holzbaur for SICStus Prolog and distributed under the 17 license details below with permission from all mentioned authors. 18 19 This program is free software; you can redistribute it and/or 20 modify it under the terms of the GNU General Public License 21 as published by the Free Software Foundation; either version 2 22 of the License, or (at your option) any later version. 23 24 This program is distributed in the hope that it will be useful, 25 but WITHOUT ANY WARRANTY; without even the implied warranty of 26 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 27 GNU General Public License for more details. 28 29 You should have received a copy of the GNU Lesser General Public 30 License along with this library; if not, write to the Free Software 31 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 32 33 As a special exception, if you link this library with other files, 34 compiled with a Free Software compiler, to produce an executable, this 35 library does not by itself cause the resulting executable to be covered 36 by the GNU General Public License. This exception does not however 37 invalidate any other reasons why the executable file might be covered by 38 the GNU General Public License. 39*/ 40 41 42:- module(fourmotz_r, 43 [ 44 fm_elim/3 45 ]). 46:- use_module(bv_r, 47 [ 48 allvars/2, 49 basis_add/2, 50 detach_bounds/1, 51 pivot/5, 52 var_with_def_intern/4 53 ]). 54:- use_module('../clpqr/class', 55 [ 56 class_allvars/2 57 ]). 58:- use_module('../clpqr/project', 59 [ 60 drop_dep/1, 61 drop_dep_one/1, 62 make_target_indep/2 63 ]). 64:- use_module('../clpqr/redund', 65 [ 66 redundancy_vars/1 67 ]). 68:- use_module(store_r, 69 [ 70 add_linear_11/3, 71 add_linear_f1/4, 72 indep/2, 73 nf_coeff_of/3, 74 normalize_scalar/2 75 ]). 76 77 78 79fm_elim(Vs,Target,Pivots) :- 80 prefilter(Vs,Vsf), 81 fm_elim_int(Vsf,Target,Pivots). 82 83% prefilter(Vars,Res) 84% 85% filters out target variables and variables that do not occur in bounded linear equations. 86% Stores that the variables in Res are to be kept independent. 87 88prefilter([],[]). 89prefilter([V|Vs],Res) :- 90 ( get_attr(V,clpqr_itf,Att), 91 arg(9,Att,n), 92 occurs(V) % V is a nontarget variable that occurs in a bounded linear equation 93 -> Res = [V|Tail], 94 setarg(10,Att,keep_indep), 95 prefilter(Vs,Tail) 96 ; prefilter(Vs,Res) 97 ). 98 99% 100% the target variables are marked with an attribute, and we get a list 101% of them as an argument too 102% 103fm_elim_int([],_,Pivots) :- % done 104 unkeep(Pivots). 105fm_elim_int(Vs,Target,Pivots) :- 106 Vs = [_|_], 107 ( best(Vs,Best,Rest) 108 -> occurences(Best,Occ), 109 elim_min(Best,Occ,Target,Pivots,NewPivots) 110 ; % give up 111 NewPivots = Pivots, 112 Rest = [] 113 ), 114 fm_elim_int(Rest,Target,NewPivots). 115 116% best(Vs,Best,Rest) 117% 118% Finds the variable with the best result (lowest Delta) in fm_cp_filter 119% and returns the other variables in Rest. 120 121best(Vs,Best,Rest) :- 122 findall(Delta-N,fm_cp_filter(Vs,Delta,N),Deltas), 123 keysort(Deltas,[_-N|_]), 124 select_nth(Vs,N,Best,Rest). 125 126% fm_cp_filter(Vs,Delta,N) 127% 128% For an indepenent variable V in Vs, which is the N'th element in Vs, 129% find how many inequalities are generated when this variable is eliminated. 130% Note that target variables and variables that only occur in unbounded equations 131% should have been removed from Vs via prefilter/2 132 133fm_cp_filter(Vs,Delta,N) :- 134 length(Vs,Len), % Len = number of variables in Vs 135 mem(Vs,X,Vst), % Selects a variable X in Vs, Vst is the list of elements after X in Vs 136 get_attr(X,clpqr_itf,Att), 137 arg(4,Att,lin(Lin)), 138 arg(5,Att,order(OrdX)), 139 arg(9,Att,n), % no target variable 140 indep(Lin,OrdX), % X is an independent variable 141 occurences(X,Occ), 142 Occ = [_|_], 143 cp_card(Occ,0,Lnew), 144 length(Occ,Locc), 145 Delta is Lnew-Locc, 146 length(Vst,Vstl), 147 N is Len-Vstl. % X is the Nth element in Vs 148 149% mem(Xs,X,XsT) 150% 151% If X is a member of Xs, XsT is the list of elements after X in Xs. 152 153mem([X|Xs],X,Xs). 154mem([_|Ys],X,Xs) :- mem(Ys,X,Xs). 155 156% select_nth(List,N,Nth,Others) 157% 158% Selects the N th element of List, stores it in Nth and returns the rest of the list in Others. 159 160select_nth(List,N,Nth,Others) :- 161 select_nth(List,1,N,Nth,Others). 162 163select_nth([X|Xs],N,N,X,Xs) :- !. 164select_nth([Y|Ys],M,N,X,[Y|Xs]) :- 165 M1 is M+1, 166 select_nth(Ys,M1,N,X,Xs). 167 168% 169% fm_detach + reverse_pivot introduce indep t_none, which 170% invalidates the invariants 171% 172elim_min(V,Occ,Target,Pivots,NewPivots) :- 173 crossproduct(Occ,New,[]), 174 activate_crossproduct(New), 175 reverse_pivot(Pivots), 176 fm_detach(Occ), 177 allvars(V,All), 178 redundancy_vars(All), % only for New \== [] 179 make_target_indep(Target,NewPivots), 180 drop_dep(All). 181 182% 183% restore NF by reverse pivoting 184% 185reverse_pivot([]). 186reverse_pivot([I:D|Ps]) :- 187 get_attr(D,clpqr_itf,AttD), 188 arg(2,AttD,type(Dt)), 189 setarg(11,AttD,n), % no longer 190 get_attr(I,clpqr_itf,AttI), 191 arg(2,AttI,type(It)), 192 arg(5,AttI,order(OrdI)), 193 arg(6,AttI,class(ClI)), 194 pivot(D,ClI,OrdI,Dt,It), 195 reverse_pivot(Ps). 196 197% unkeep(Pivots) 198% 199% 200 201unkeep([]). 202unkeep([_:D|Ps]) :- 203 get_attr(D,clpqr_itf,Att), 204 setarg(11,Att,n), 205 drop_dep_one(D), 206 unkeep(Ps). 207 208 209% 210% All we drop are bounds 211% 212fm_detach( []). 213fm_detach([V:_|Vs]) :- 214 detach_bounds(V), 215 fm_detach(Vs). 216 217% activate_crossproduct(Lst) 218% 219% For each inequality Lin =< 0 (or Lin < 0) in Lst, a new variable is created: 220% Var = Lin and Var =< 0 (or Var < 0). Var is added to the basis. 221 222activate_crossproduct([]). 223activate_crossproduct([lez(Strict,Lin)|News]) :- 224 var_with_def_intern(t_u(0.0),Var,Lin,Strict), 225 % Var belongs to same class as elements in Lin 226 basis_add(Var,_), 227 activate_crossproduct(News). 228 229% ------------------------------------------------------------------------------ 230 231% crossproduct(Lst,Res,ResTail) 232% 233% See crossproduct/4 234% This predicate each time puts the next element of Lst as First in crossproduct/4 235% and lets the rest be Next. 236 237crossproduct([]) --> []. 238crossproduct([A|As]) --> 239 crossproduct(As,A), 240 crossproduct(As). 241 242% crossproduct(Next,First,Res,ResTail) 243% 244% Eliminates a variable in linear equations First + Next and stores the generated 245% inequalities in Res. 246% Let's say A:K1 = First and B:K2 = first equation in Next. 247% A = ... + K1*V + ... 248% B = ... + K2*V + ... 249% Let K = -K2/K1 250% then K*A + B = ... + 0*V + ... 251% from the bounds of A and B, via cross_lower/7 and cross_upper/7, new inequalities 252% are generated. Then the same is done for B:K2 = next element in Next. 253 254crossproduct([],_) --> []. 255crossproduct([B:Kb|Bs],A:Ka) --> 256 { 257 get_attr(A,clpqr_itf,AttA), 258 arg(2,AttA,type(Ta)), 259 arg(3,AttA,strictness(Sa)), 260 arg(4,AttA,lin(LinA)), 261 get_attr(B,clpqr_itf,AttB), 262 arg(2,AttB,type(Tb)), 263 arg(3,AttB,strictness(Sb)), 264 arg(4,AttB,lin(LinB)), 265 K is -Kb/Ka, 266 add_linear_f1(LinA,K,LinB,Lin) % Lin doesn't contain the target variable anymore 267 }, 268 ( { K > 1.0e-10 } % K > 0: signs were opposite 269 -> { Strict is Sa \/ Sb }, 270 cross_lower(Ta,Tb,K,Lin,Strict), 271 cross_upper(Ta,Tb,K,Lin,Strict) 272 ; % La =< A =< Ua -> -Ua =< -A =< -La 273 { 274 flip(Ta,Taf), 275 flip_strict(Sa,Saf), 276 Strict is Saf \/ Sb 277 }, 278 cross_lower(Taf,Tb,K,Lin,Strict), 279 cross_upper(Taf,Tb,K,Lin,Strict) 280 ), 281 crossproduct(Bs,A:Ka). 282 283% cross_lower(Ta,Tb,K,Lin,Strict,Res,ResTail) 284% 285% Generates a constraint following from the bounds of A and B. 286% When A = LinA and B = LinB then Lin = K*LinA + LinB. Ta is the type 287% of A and Tb is the type of B. Strict is the union of the strictness 288% of A and B. If K is negative, then Ta should have been flipped (flip/2). 289% The idea is that if La =< A =< Ua and Lb =< B =< Ub (=< can also be <) 290% then if K is positive, K*La + Lb =< K*A + B =< K*Ua + Ub. 291% if K is negative, K*Ua + Lb =< K*A + B =< K*La + Ub. 292% This predicate handles the first inequality and adds it to Res in the form 293% lez(Sl,Lhs) meaning K*La + Lb - (K*A + B) =< 0 or K*Ua + Lb - (K*A + B) =< 0 294% with Sl being the strictness and Lhs the lefthandside of the equation. 295% See also cross_upper/7 296 297cross_lower(Ta,Tb,K,Lin,Strict) --> 298 { 299 lower(Ta,La), 300 lower(Tb,Lb), 301 !, 302 L is K*La+Lb, 303 normalize_scalar(L,Ln), 304 add_linear_f1(Lin,-1.0,Ln,Lhs), 305 Sl is Strict >> 1 % normalize to upper bound 306 }, 307 [ lez(Sl,Lhs) ]. 308cross_lower(_,_,_,_,_) --> []. 309 310% cross_upper(Ta,Tb,K,Lin,Strict,Res,ResTail) 311% 312% See cross_lower/7 313% This predicate handles the second inequality: 314% -(K*Ua + Ub) + K*A + B =< 0 or -(K*La + Ub) + K*A + B =< 0 315 316cross_upper(Ta,Tb,K,Lin,Strict) --> 317 { 318 upper(Ta,Ua), 319 upper(Tb,Ub), 320 !, 321 U is -(K*Ua+Ub), 322 normalize_scalar(U,Un), 323 add_linear_11(Un,Lin,Lhs), 324 Su is Strict /\ 1 % normalize to upper bound 325 }, 326 [ lez(Su,Lhs) ]. 327cross_upper(_,_,_,_,_) --> []. 328 329% lower(Type,Lowerbound) 330% 331% Returns the lowerbound of type Type if it has one. 332% E.g. if type = t_l(L) then Lowerbound is L, 333% if type = t_lU(L,U) then Lowerbound is L, 334% if type = t_u(U) then fails 335 336lower(t_l(L),L). 337lower(t_lu(L,_),L). 338lower(t_L(L),L). 339lower(t_Lu(L,_),L). 340lower(t_lU(L,_),L). 341 342% upper(Type,Upperbound) 343% 344% Returns the upperbound of type Type if it has one. 345% See lower/2 346 347upper(t_u(U),U). 348upper(t_lu(_,U),U). 349upper(t_U(U),U). 350upper(t_Lu(_,U),U). 351upper(t_lU(_,U),U). 352 353% flip(Type,FlippedType) 354% 355% Flips the lower and upperbound, so the old lowerbound becomes the new upperbound and 356% vice versa. 357 358flip(t_l(X),t_u(X)). 359flip(t_u(X),t_l(X)). 360flip(t_lu(X,Y),t_lu(Y,X)). 361flip(t_L(X),t_u(X)). 362flip(t_U(X),t_l(X)). 363flip(t_lU(X,Y),t_lu(Y,X)). 364flip(t_Lu(X,Y),t_lu(Y,X)). 365 366% flip_strict(Strict,FlippedStrict) 367% 368% Does what flip/2 does, but for the strictness. 369 370flip_strict(0,0). 371flip_strict(1,2). 372flip_strict(2,1). 373flip_strict(3,3). 374 375% cp_card(Lst,CountIn,CountOut) 376% 377% Counts the number of bounds that may generate an inequality in 378% crossproduct/3 379 380cp_card([],Ci,Ci). 381cp_card([A|As],Ci,Co) :- 382 cp_card(As,A,Ci,Cii), 383 cp_card(As,Cii,Co). 384 385% cp_card(Next,First,CountIn,CountOut) 386% 387% Counts the number of bounds that may generate an inequality in 388% crossproduct/4. 389 390cp_card([],_,Ci,Ci). 391cp_card([B:Kb|Bs],A:Ka,Ci,Co) :- 392 get_attr(A,clpqr_itf,AttA), 393 arg(2,AttA,type(Ta)), 394 get_attr(B,clpqr_itf,AttB), 395 arg(2,AttB,type(Tb)), 396 K is -Kb/Ka, 397 ( K > 1.0e-10 % K > 0: signs were opposite 398 -> cp_card_lower(Ta,Tb,Ci,Cii), 399 cp_card_upper(Ta,Tb,Cii,Ciii) 400 ; flip(Ta,Taf), 401 cp_card_lower(Taf,Tb,Ci,Cii), 402 cp_card_upper(Taf,Tb,Cii,Ciii) 403 ), 404 cp_card(Bs,A:Ka,Ciii,Co). 405 406% cp_card_lower(TypeA,TypeB,SIn,SOut) 407% 408% SOut = SIn + 1 if both TypeA and TypeB have a lowerbound. 409 410cp_card_lower(Ta,Tb,Si,So) :- 411 lower(Ta,_), 412 lower(Tb,_), 413 !, 414 So is Si+1. 415cp_card_lower(_,_,Si,Si). 416 417% cp_card_upper(TypeA,TypeB,SIn,SOut) 418% 419% SOut = SIn + 1 if both TypeA and TypeB have an upperbound. 420 421cp_card_upper(Ta,Tb,Si,So) :- 422 upper(Ta,_), 423 upper(Tb,_), 424 !, 425 So is Si+1. 426cp_card_upper(_,_,Si,Si). 427 428% ------------------------------------------------------------------------------ 429 430% occurences(V,Occ) 431% 432% Returns in Occ the occurrences of variable V in the linear equations of dependent variables 433% with bound =\= t_none in the form of D:K where D is a dependent variable and K is the scalar 434% of V in the linear equation of D. 435 436occurences(V,Occ) :- 437 get_attr(V,clpqr_itf,Att), 438 arg(5,Att,order(OrdV)), 439 arg(6,Att,class(C)), 440 class_allvars(C,All), 441 occurences(All,OrdV,Occ). 442 443% occurences(De,OrdV,Occ) 444% 445% Returns in Occ the occurrences of variable V with order OrdV in the linear equations of 446% dependent variables De with bound =\= t_none in the form of D:K where D is a dependent 447% variable and K is the scalar of V in the linear equation of D. 448 449occurences(De,_,[]) :- 450 var(De), 451 !. 452occurences([D|De],OrdV,Occ) :- 453 ( get_attr(D,clpqr_itf,Att), 454 arg(2,Att,type(Type)), 455 arg(4,Att,lin(Lin)), 456 occ_type_filter(Type), 457 nf_coeff_of(Lin,OrdV,K) 458 -> Occ = [D:K|Occt], 459 occurences(De,OrdV,Occt) 460 ; occurences(De,OrdV,Occ) 461 ). 462 463% occ_type_filter(Type) 464% 465% Succeeds when Type is any other type than t_none. Is used in occurences/3 and occurs/2 466 467occ_type_filter(t_l(_)). 468occ_type_filter(t_u(_)). 469occ_type_filter(t_lu(_,_)). 470occ_type_filter(t_L(_)). 471occ_type_filter(t_U(_)). 472occ_type_filter(t_lU(_,_)). 473occ_type_filter(t_Lu(_,_)). 474 475% occurs(V) 476% 477% Checks whether variable V occurs in a linear equation of a dependent variable with a bound 478% =\= t_none. 479 480occurs(V) :- 481 get_attr(V,clpqr_itf,Att), 482 arg(5,Att,order(OrdV)), 483 arg(6,Att,class(C)), 484 class_allvars(C,All), 485 occurs(All,OrdV). 486 487% occurs(De,OrdV) 488% 489% Checks whether variable V with order OrdV occurs in a linear equation of any dependent variable 490% in De with a bound =\= t_none. 491 492occurs(De,_) :- 493 var(De), 494 !, 495 fail. 496occurs([D|De],OrdV) :- 497 ( get_attr(D,clpqr_itf,Att), 498 arg(2,Att,type(Type)), 499 arg(4,Att,lin(Lin)), 500 occ_type_filter(Type), 501 nf_coeff_of(Lin,OrdV,_) 502 -> true 503 ; occurs(De,OrdV) 504 )