1/* 2 3 Part of CLP(Q) (Constraint Logic Programming over Rationals) 4 5 Author: Leslie De Koninck 6 E-mail: Leslie.DeKoninck@cs.kuleuven.be 7 WWW: http://www.swi-prolog.org 8 http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 9 Copyright (C): 2006, K.U. Leuven and 10 1992-1995, Austrian Research Institute for 11 Artificial Intelligence (OFAI), 12 Vienna, Austria 13 14 This software is based on CLP(Q,R) by Christian Holzbaur for SICStus 15 Prolog and distributed under the license details below with permission from 16 all mentioned authors. 17 18 This program is free software; you can redistribute it and/or 19 modify it under the terms of the GNU General Public License 20 as published by the Free Software Foundation; either version 2 21 of the License, or (at your option) any later version. 22 23 This program is distributed in the hope that it will be useful, 24 but WITHOUT ANY WARRANTY; without even the implied warranty of 25 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 26 GNU General Public License for more details. 27 28 You should have received a copy of the GNU Lesser General Public 29 License along with this library; if not, write to the Free Software 30 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 31 32 As a special exception, if you link this library with other files, 33 compiled with a Free Software compiler, to produce an executable, this 34 library does not by itself cause the resulting executable to be covered 35 by the GNU General Public License. This exception does not however 36 invalidate any other reasons why the executable file might be covered by 37 the GNU General Public License. 38*/ 39 40:- module(bv_q, 41 [ 42 allvars/2, 43 backsubst/3, 44 backsubst_delta/4, 45 basis_add/2, 46 dec_step/2, 47 deref/2, 48 deref_var/2, 49 detach_bounds/1, 50 detach_bounds_vlv/5, 51 determine_active_dec/1, 52 determine_active_inc/1, 53 dump_var/6, 54 dump_nz/5, 55 export_binding/1, 56 get_or_add_class/2, 57 inc_step/2, 58 intro_at/3, 59 iterate_dec/2, 60 lb/3, 61 pivot_a/4, 62 pivot/5, 63 rcbl_status/6, 64 reconsider/1, 65 same_class/2, 66 solve/1, 67 solve_ord_x/3, 68 ub/3, 69 unconstrained/4, 70 var_intern/2, 71 var_intern/3, 72 var_with_def_assign/2, 73 var_with_def_intern/4, 74 maximize/1, 75 minimize/1, 76 sup/2, 77 sup/4, 78 inf/2, 79 inf/4, 80 'solve_<'/1, 81 'solve_=<'/1, 82 'solve_=\\='/1, 83 log_deref/4 84 ]). 85:- use_module(store_q, 86 [ 87 add_linear_11/3, 88 add_linear_f1/4, 89 add_linear_ff/5, 90 delete_factor/4, 91 indep/2, 92 isolate/3, 93 nf2sum/3, 94 nf_rhs_x/4, 95 nf_substitute/4, 96 normalize_scalar/2, 97 mult_hom/3, 98 mult_linear_factor/3 99 ]). 100:- use_module('../clpqr/class', 101 [ 102 class_allvars/2, 103 class_basis/2, 104 class_basis_add/3, 105 class_basis_drop/2, 106 class_basis_pivot/3, 107 class_new/5 108 ]). 109:- use_module(ineq_q, 110 [ 111 ineq/4 112 ]). 113:- use_module(nf_q, 114 [ 115 {}/1, 116 split/3, 117 wait_linear/3 118 ]). 119:- use_module(bb_q, 120 [ 121 vertex_value/2 122 ]). 123:- use_module(library(ordsets), 124 [ 125 ord_add_element/3 126 ]). 127 128% For the rhs maint. the following events are important: 129% 130% -) introduction of an indep var at active bound B 131% -) narrowing of active bound 132% -) swap active bound 133% -) pivot 134% 135 136% a variables bound (L/U) can have the states: 137% 138% -) t_none no bounds 139% -) t_l inactive lower bound 140% -) t_u inactive upper bound 141% -) t_L active lower bound 142% -) t_U active upper bound 143% -) t_lu inactive lower and upper bound 144% -) t_Lu active lower bound and inactive upper bound 145% -) t_lU inactive lower bound and active upper bound 146 147% ----------------------------------- deref ----------------------------------- 148% 149 150% deref(Lin,Lind) 151% 152% Makes a linear equation of the form [v(I,[])|H] into a solvable linear 153% equation. 154% If the variables are new, they are initialized with the linear equation X=X. 155 156deref(Lin,Lind) :- 157 split(Lin,H,I), 158 normalize_scalar(I,Nonvar), 159 length(H,Len), 160 log_deref(Len,H,[],Restd), 161 add_linear_11(Nonvar,Restd,Lind). 162 163% log_deref(Len,[Vs|VsTail],VsTail,Res) 164% 165% Logarithmically converts a linear equation in normal form ([v(_,_)|_]) into a 166% linear equation in solver form ([I,R,K*X|_]). Res contains the result, Len is 167% the length of the part to convert and [Vs|VsTail] is a difference list 168% containing the equation in normal form. 169 170log_deref(0,Vs,Vs,Lin) :- 171 !, 172 Lin = [0,0]. 173log_deref(1,[v(K,[X^1])|Vs],Vs,Lin) :- 174 !, 175 deref_var(X,Lx), 176 mult_linear_factor(Lx,K,Lin). 177log_deref(2,[v(Kx,[X^1]),v(Ky,[Y^1])|Vs],Vs,Lin) :- 178 !, 179 deref_var(X,Lx), 180 deref_var(Y,Ly), 181 add_linear_ff(Lx,Kx,Ly,Ky,Lin). 182log_deref(N,V0,V2,Lin) :- 183 P is N >> 1, 184 Q is N - P, 185 log_deref(P,V0,V1,Lp), 186 log_deref(Q,V1,V2,Lq), 187 add_linear_11(Lp,Lq,Lin). 188 189% deref_var(X,Lin) 190% 191% Returns the equation of variable X. If X is a new variable, a new equation 192% X = X is made. 193 194deref_var(X,Lin) :- 195 ( get_attr(X,itf,Att) 196 -> ( \+ arg(1,Att,clpq) 197 -> throw(error(permission_error('mix CLP(Q) variables with', 198 'CLP(R) variables:',X),context(_))) 199 ; arg(4,Att,lin(Lin)) 200 -> true 201 ; setarg(2,Att,type(t_none)), 202 setarg(3,Att,strictness(0)), 203 Lin = [0,0,l(X*1,Ord)], 204 setarg(4,Att,lin(Lin)), 205 setarg(5,Att,order(Ord)) 206 ) 207 ; Lin = [0,0,l(X*1,Ord)], 208 put_attr(X,itf,t(clpq,type(t_none),strictness(0), 209 lin(Lin),order(Ord),n,n,n,n,n,n)) 210 ). 211 212% TODO 213% 214% 215 216var_with_def_assign(Var,Lin) :- 217 Lin = [I,_|Hom], 218 ( Hom = [] 219 -> % X=k 220 Var = I 221 ; Hom = [l(V*K,_)|Cs] 222 -> ( Cs = [], 223 K =:= 1, 224 I =:= 0 225 -> % X=Y 226 Var = V 227 ; % general case 228 var_with_def_intern(t_none,Var,Lin,0) 229 ) 230 ). 231 232% var_with_def_intern(Type,Var,Lin,Strictness) 233% 234% Makes Lin the linear equation of new variable Var, makes all variables of 235% Lin, and Var of the same class and bounds Var by type(Type) and 236% strictness(Strictness) 237 238var_with_def_intern(Type,Var,Lin,Strict) :- 239 put_attr(Var,itf,t(clpq,type(Type),strictness(Strict),lin(Lin), 240 order(_),n,n,n,n,n,n)), % check uses 241 Lin = [_,_|Hom], 242 get_or_add_class(Var,Class), 243 same_class(Hom,Class). 244 245% TODO 246% 247% 248 249var_intern(Type,Var,Strict) :- 250 put_attr(Var,itf,t(clpq,type(Type),strictness(Strict), 251 lin([0,0,l(Var*1,Ord)]),order(Ord),n,n,n,n,n,n)), 252 get_or_add_class(Var,_Class). 253 254% TODO 255% 256% 257 258var_intern(Var,Class) :- % for ordered/1 but otherwise free vars 259 get_attr(Var,itf,Att), 260 arg(2,Att,type(_)), 261 arg(4,Att,lin(_)), 262 !, 263 get_or_add_class(Var,Class). 264var_intern(Var,Class) :- 265 put_attr(Var,itf,t(clpq,type(t_none),strictness(0), 266 lin([0,0,l(Var*1,Ord)]),order(Ord),n,n,n,n,n,n)), 267 get_or_add_class(Var,Class). 268 269% ----------------------------------------------------------------------------- 270 271% export_binding(Lst) 272% 273% Binds variables X to Y where Lst contains elements of the form [X-Y]. 274 275export_binding([]). 276export_binding([X-Y|Gs]) :- 277 Y = X, 278 export_binding(Gs). 279 280% 'solve_='(Nf) 281% 282% Solves linear equation Nf = 0 where Nf is in normal form. 283 284'solve_='(Nf) :- 285 deref(Nf,Nfd), % dereferences and turns Nf into solvable form Nfd 286 solve(Nfd). 287 288% 'solve_=\\='(Nf) 289% 290% Solves linear inequality Nf =\= 0 where Nf is in normal form. 291 292'solve_=\\='(Nf) :- 293 deref(Nf,Lind), % dereferences and turns Nf into solvable form Lind 294 Lind = [Inhom,_|Hom], 295 ( Hom = [] 296 -> Inhom =\= 0 297 ; % make new variable Nz = Lind 298 var_with_def_intern(t_none,Nz,Lind,0), 299 % make Nz nonzero 300 get_attr(Nz,itf,Att), 301 setarg(8,Att,nonzero) 302 ). 303 304% 'solve_<'(Nf) 305% 306% Solves linear inequality Nf < 0 where Nf is in normal form. 307 308'solve_<'(Nf) :- 309 split(Nf,H,I), 310 ineq(H,I,Nf,strict). 311 312% 'solve_=<'(Nf) 313% 314% Solves linear inequality Nf =< 0 where Nf is in normal form. 315 316'solve_=<'(Nf) :- 317 split(Nf,H,I), 318 ineq(H,I,Nf,nonstrict). 319 320maximize(Term) :- 321 minimize(-Term). 322 323% 324% This is NOT coded as minimize(Expr) :- inf(Expr,Expr). 325% 326% because the new version of inf/2 only visits 327% the vertex where the infimum is assumed and returns 328% to the 'current' vertex via backtracking. 329% The rationale behind this construction is to eliminate 330% all garbage in the solver data structures produced by 331% the pivots on the way to the extremal point caused by 332% {inf,sup}/{2,4}. 333% 334% If we are after the infimum/supremum for minimizing/maximizing, 335% this strategy may have adverse effects on performance because 336% the simplex algorithm is forced to re-discover the 337% extremal vertex through the equation {Inf =:= Expr}. 338% 339% Thus the extra code for {minimize,maximize}/1. 340% 341% In case someone comes up with an example where 342% 343% inf(Expr,Expr) 344% 345% outperforms the provided formulation for minimize - so be it. 346% Both forms are available to the user. 347% 348minimize(Term) :- 349 wait_linear(Term,Nf,minimize_lin(Nf)). 350 351% minimize_lin(Lin) 352% 353% Minimizes the linear expression Lin. It does so by making a new 354% variable Dep and minimizes its value. 355 356minimize_lin(Lin) :- 357 deref(Lin,Lind), 358 var_with_def_intern(t_none,Dep,Lind,0), 359 determine_active_dec(Lind), 360 iterate_dec(Dep,Inf), 361 { Dep =:= Inf }. 362 363sup(Expression,Sup) :- 364 sup(Expression,Sup,[],[]). 365 366sup(Expression,Sup,Vector,Vertex) :- 367 inf(-Expression,-Sup,Vector,Vertex). 368 369inf(Expression,Inf) :- 370 inf(Expression,Inf,[],[]). 371 372inf(Expression,Inf,Vector,Vertex) :- 373 % wait until Expression becomes linear, Nf contains linear Expression 374 % in normal form 375 wait_linear(Expression,Nf,inf_lin(Nf,Inf,Vector,Vertex)). 376 377inf_lin(Lin,_,Vector,_) :- 378 deref(Lin,Lind), 379 var_with_def_intern(t_none,Dep,Lind,0), % make new variable Dep = Lind 380 determine_active_dec(Lind), % minimizes Lind 381 iterate_dec(Dep,Inf), 382 vertex_value(Vector,Values), 383 nb_setval(inf,[Inf|Values]), 384 fail. 385inf_lin(_,Infimum,_,Vertex) :- 386 catch(nb_getval(inf,L),_,fail), 387 nb_delete(inf), 388 assign([Infimum|Vertex],L). 389 390% assign(L1,L2) 391% 392% The elements of L1 are pairwise assigned to the elements of L2 393% by means of asserting {X =:= Y} where X is an element of L1 and Y 394% is the corresponding element of L2. 395 396assign([],[]). 397assign([X|Xs],[Y|Ys]) :- 398 {X =:= Y}, % more defensive/expressive than X=Y 399 assign(Xs,Ys). 400 401% --------------------------------- optimization ------------------------------ 402% 403% The _sn(S) =< 0 row might be temporarily infeasible. 404% We use reconsider/1 to fix this. 405% 406% s(S) e [_,0] = d +xi ... -xj, Rhs > 0 so we want to decrease s(S) 407% 408% positive xi would have to be moved towards their lower bound, 409% negative xj would have to be moved towards their upper bound, 410% 411% the row s(S) does not limit the lower bound of xi 412% the row s(S) does not limit the upper bound of xj 413% 414% a) if some other row R is limiting xk, we pivot(R,xk), 415% s(S) will decrease and get more feasible until (b) 416% b) if there is no limiting row for some xi: we pivot(s(S),xi) 417% xj: we pivot(s(S),xj) 418% which cures the infeasibility in one step 419% 420 421 422% iterate_dec(OptVar,Opt) 423% 424% Decreases the bound on the variables of the linear equation of OptVar as much 425% as possible and returns the resulting optimal bound in Opt. Fails if for some 426% variable, a status of unlimited is found. 427 428iterate_dec(OptVar,Opt) :- 429 get_attr(OptVar,itf,Att), 430 arg(4,Att,lin([I,R|H])), 431 dec_step(H,Status), 432 ( Status = applied 433 -> iterate_dec(OptVar,Opt) 434 ; Status = optimum, 435 Opt is R + I 436 ). 437 438% iterate_inc(OptVar,Opt) 439% 440% Increases the bound on the variables of the linear equation of OptVar as much 441% as possible and returns the resulting optimal bound in Opt. Fails if for some 442% variable, a status of unlimited is found. 443 444iterate_inc(OptVar,Opt) :- 445 get_attr(OptVar,itf,Att), 446 arg(4,Att,lin([I,R|H])), 447 inc_step(H,Status), 448 ( Status = applied 449 -> iterate_inc(OptVar,Opt) 450 ; Status = optimum, 451 Opt is R + I 452 ). 453 454% 455% Status = {optimum,unlimited(Indep,DepT),applied} 456% If Status = optimum, the tables have not been changed at all. 457% Searches left to right, does not try to find the 'best' pivot 458% Therefore we might discover unboundedness only after a few pivots 459% 460 461 462dec_step_cont([],optimum,Cont,Cont). 463dec_step_cont([l(V*K,OrdV)|Vs],Status,ContIn,ContOut) :- 464 get_attr(V,itf,Att), 465 arg(2,Att,type(W)), 466 arg(6,Att,class(Class)), 467 ( dec_step_2_cont(W,l(V*K,OrdV),Class,Status,ContIn,ContOut) 468 -> true 469 ; dec_step_cont(Vs,Status,ContIn,ContOut) 470 ). 471 472inc_step_cont([],optimum,Cont,Cont). 473inc_step_cont([l(V*K,OrdV)|Vs],Status,ContIn,ContOut) :- 474 get_attr(V,itf,Att), 475 arg(2,Att,type(W)), 476 arg(6,Att,class(Class)), 477 ( inc_step_2_cont(W,l(V*K,OrdV),Class,Status,ContIn,ContOut) 478 -> true 479 ; inc_step_cont(Vs,Status,ContIn,ContOut) 480 ). 481 482dec_step_2_cont(t_U(U),l(V*K,OrdV),Class,Status,ContIn,ContOut) :- 483 K > 0, 484 ( lb(Class,OrdV,Vub-Vb-_) 485 -> % found a lower bound 486 Status = applied, 487 pivot_a(Vub,V,Vb,t_u(U)), 488 replace_in_cont(ContIn,Vub,V,ContOut) 489 ; Status = unlimited(V,t_u(U)), 490 ContIn = ContOut 491 ). 492dec_step_2_cont(t_lU(L,U),l(V*K,OrdV),Class,applied,ContIn,ContOut) :- 493 K > 0, 494 Init is L - U, 495 class_basis(Class,Deps), 496 lb(Deps,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_), 497 pivot_b(Vub,V,Vb,t_lu(L,U)), 498 replace_in_cont(ContIn,Vub,V,ContOut). 499dec_step_2_cont(t_L(L),l(V*K,OrdV),Class,Status,ContIn,ContOut) :- 500 K < 0, 501 ( ub(Class,OrdV,Vub-Vb-_) 502 -> Status = applied, 503 pivot_a(Vub,V,Vb,t_l(L)), 504 replace_in_cont(ContIn,Vub,V,ContOut) 505 ; Status = unlimited(V,t_l(L)), 506 ContIn = ContOut 507 ). 508dec_step_2_cont(t_Lu(L,U),l(V*K,OrdV),Class,applied,ContIn,ContOut) :- 509 K < 0, 510 Init is U - L, 511 class_basis(Class,Deps), 512 ub(Deps,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_), 513 pivot_b(Vub,V,Vb,t_lu(L,U)), 514 replace_in_cont(ContIn,Vub,V,ContOut). 515dec_step_2_cont(t_none,l(V*_,_),_,unlimited(V,t_none),Cont,Cont). 516 517 518 519inc_step_2_cont(t_U(U),l(V*K,OrdV),Class,Status,ContIn,ContOut) :- 520 K < 0, 521 ( lb(Class,OrdV,Vub-Vb-_) 522 -> Status = applied, 523 pivot_a(Vub,V,Vb,t_u(U)), 524 replace_in_cont(ContIn,Vub,V,ContOut) 525 ; Status = unlimited(V,t_u(U)), 526 ContIn = ContOut 527 ). 528inc_step_2_cont(t_lU(L,U),l(V*K,OrdV),Class,applied,ContIn,ContOut) :- 529 K < 0, 530 Init is L - U, 531 class_basis(Class,Deps), 532 lb(Deps,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_), 533 pivot_b(Vub,V,Vb,t_lu(L,U)), 534 replace_in_cont(ContIn,Vub,V,ContOut). 535inc_step_2_cont(t_L(L),l(V*K,OrdV),Class,Status,ContIn,ContOut) :- 536 K > 0, 537 ( ub(Class,OrdV,Vub-Vb-_) 538 -> Status = applied, 539 pivot_a(Vub,V,Vb,t_l(L)), 540 replace_in_cont(ContIn,Vub,V,ContOut) 541 ; Status = unlimited(V,t_l(L)), 542 ContIn = ContOut 543 ). 544inc_step_2_cont(t_Lu(L,U),l(V*K,OrdV),Class,applied,ContIn,ContOut) :- 545 K > 0, 546 Init is U - L, 547 class_basis(Class,Deps), 548 ub(Deps,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_), 549 pivot_b(Vub,V,Vb,t_lu(L,U)), 550 replace_in_cont(ContIn,Vub,V,ContOut). 551inc_step_2_cont(t_none,l(V*_,_),_,unlimited(V,t_none),Cont,Cont). 552 553replace_in_cont([],_,_,[]). 554replace_in_cont([H1|T1],X,Y,[H2|T2]) :- 555 ( H1 == X 556 -> H2 = Y, 557 T1 = T2 558 ; H2 = H1, 559 replace_in_cont(T1,X,Y,T2) 560 ). 561 562dec_step([],optimum). 563dec_step([l(V*K,OrdV)|Vs],Status) :- 564 get_attr(V,itf,Att), 565 arg(2,Att,type(W)), 566 arg(6,Att,class(Class)), 567 ( dec_step_2(W,l(V*K,OrdV),Class,Status) 568 -> true 569 ; dec_step(Vs,Status) 570 ). 571 572dec_step_2(t_U(U),l(V*K,OrdV),Class,Status) :- 573 K > 0, 574 ( lb(Class,OrdV,Vub-Vb-_) 575 -> % found a lower bound 576 Status = applied, 577 pivot_a(Vub,V,Vb,t_u(U)) 578 ; Status = unlimited(V,t_u(U)) 579 ). 580dec_step_2(t_lU(L,U),l(V*K,OrdV),Class,applied) :- 581 K > 0, 582 Init is L - U, 583 class_basis(Class,Deps), 584 lb(Deps,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_), 585 pivot_b(Vub,V,Vb,t_lu(L,U)). 586dec_step_2(t_L(L),l(V*K,OrdV),Class,Status) :- 587 K < 0, 588 ( ub(Class,OrdV,Vub-Vb-_) 589 -> Status = applied, 590 pivot_a(Vub,V,Vb,t_l(L)) 591 ; Status = unlimited(V,t_l(L)) 592 ). 593dec_step_2(t_Lu(L,U),l(V*K,OrdV),Class,applied) :- 594 K < 0, 595 Init is U - L, 596 class_basis(Class,Deps), 597 ub(Deps,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_), 598 pivot_b(Vub,V,Vb,t_lu(L,U)). 599dec_step_2(t_none,l(V*_,_),_,unlimited(V,t_none)). 600 601inc_step([],optimum). % if status has not been set yet: no changes 602inc_step([l(V*K,OrdV)|Vs],Status) :- 603 get_attr(V,itf,Att), 604 arg(2,Att,type(W)), 605 arg(6,Att,class(Class)), 606 ( inc_step_2(W,l(V*K,OrdV),Class,Status) 607 -> true 608 ; inc_step(Vs,Status) 609 ). 610 611inc_step_2(t_U(U),l(V*K,OrdV),Class,Status) :- 612 K < 0, 613 ( lb(Class,OrdV,Vub-Vb-_) 614 -> Status = applied, 615 pivot_a(Vub,V,Vb,t_u(U)) 616 ; Status = unlimited(V,t_u(U)) 617 ). 618inc_step_2(t_lU(L,U),l(V*K,OrdV),Class,applied) :- 619 K < 0, 620 Init is L - U, 621 class_basis(Class,Deps), 622 lb(Deps,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_), 623 pivot_b(Vub,V,Vb,t_lu(L,U)). 624inc_step_2(t_L(L),l(V*K,OrdV),Class,Status) :- 625 K > 0, 626 ( ub(Class,OrdV,Vub-Vb-_) 627 -> Status = applied, 628 pivot_a(Vub,V,Vb,t_l(L)) 629 ; Status = unlimited(V,t_l(L)) 630 ). 631inc_step_2(t_Lu(L,U),l(V*K,OrdV),Class,applied) :- 632 K > 0, 633 Init is U - L, 634 class_basis(Class,Deps), 635 ub(Deps,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_), 636 pivot_b(Vub,V,Vb,t_lu(L,U)). 637inc_step_2(t_none,l(V*_,_),_,unlimited(V,t_none)). 638 639% ------------------------- find the most constraining row -------------------- 640% 641% The code for the lower and the upper bound are dual versions of each other. 642% The only difference is in the orientation of the comparisons. 643% Indeps are ruled out by their types. 644% If there is no bound, this fails. 645% 646% *** The actual lb and ub on an indep variable X are [lu]b + b(X), where b(X) 647% is the value of the active bound. 648% 649% Nota bene: We must NOT consider infeasible rows as candidates to 650% leave the basis! 651% 652% ub(Class,OrdX,Ub) 653% 654% See lb/3: this is similar 655 656ub(Class,OrdX,Ub) :- 657 class_basis(Class,Deps), 658 ub_first(Deps,OrdX,Ub). 659 660% ub_first(Deps,X,Dep-W-Ub) 661% 662% Finds the tightest upperbound for variable X from the linear equations of 663% basis variables Deps, and puts the resulting bound in Ub. Dep is the basis 664% variable that generates the bound, and W is bound of that variable that has 665% to be activated to achieve this. 666 667ub_first([Dep|Deps],OrdX,Tightest) :- 668 ( get_attr(Dep,itf,Att), 669 arg(2,Att,type(Type)), 670 arg(4,Att,lin(Lin)), 671 ub_inner(Type,OrdX,Lin,W,Ub), 672 Ub >= 0 673 -> ub(Deps,OrdX,Dep-W-Ub,Tightest) 674 ; ub_first(Deps,OrdX,Tightest) 675 ). 676 677% ub(Deps,OrdX,TightestIn,TightestOut) 678% 679% See lb/4: this is similar 680 681ub([],_,T0,T0). 682ub([Dep|Deps],OrdX,T0,T1) :- 683 ( get_attr(Dep,itf,Att), 684 arg(2,Att,type(Type)), 685 arg(4,Att,lin(Lin)), 686 ub_inner(Type,OrdX,Lin,W,Ub), 687 T0 = _-Ubb, 688 Ub < Ubb, 689 Ub >= 0 690 -> ub(Deps,OrdX,Dep-W-Ub,T1) % tighter bound, use new bound 691 ; ub(Deps,OrdX,T0,T1) % no tighter bound, keep current one 692 ). 693 694% ub_inner(Type,OrdX,Lin,W,Ub) 695% 696% See lb_inner/5: this is similar 697 698ub_inner(t_l(L),OrdX,Lin,t_L(L),Ub) :- 699 nf_rhs_x(Lin,OrdX,Rhs,K), 700 K < 0, 701 Ub is (L - Rhs) rdiv K. 702ub_inner(t_u(U),OrdX,Lin,t_U(U),Ub) :- 703 nf_rhs_x(Lin,OrdX,Rhs,K), 704 K > 0, 705 Ub is (U - Rhs) rdiv K. 706ub_inner(t_lu(L,U),OrdX,Lin,W,Ub) :- 707 nf_rhs_x(Lin,OrdX,Rhs,K), 708 ( K < 0 % use lowerbound 709 -> W = t_Lu(L,U), 710 Ub = (L - Rhs) rdiv K 711 ; K > 0 % use upperbound 712 -> W = t_lU(L,U), 713 Ub = (U - Rhs) rdiv K 714 ). 715 716% lb(Class,OrdX,Lb) 717% 718% Returns in Lb how much we can lower the upperbound of X without violating 719% a bound of the basisvariables. 720% Lb has the form Dep-W-Lb with Dep the variable whose bound is violated when 721% lowering the bound for X more, W the actual bound that has to be activated 722% and Lb the amount that the upperbound can be lowered. 723% X has ordering OrdX and class Class. 724 725lb(Class,OrdX,Lb) :- 726 class_basis(Class,Deps), 727 lb_first(Deps,OrdX,Lb). 728 729% lb_first(Deps,OrdX,Tightest) 730% 731% Returns in Tightest how much we can lower the upperbound of X without 732% violating a bound of Deps. 733% Tightest has the form Dep-W-Lb with Dep the variable whose bound is violated 734% when lowering the bound for X more, W the actual bound that has to be 735% activated and Lb the amount that the upperbound can be lowered. X has 736% ordering attribute OrdX. 737 738lb_first([Dep|Deps],OrdX,Tightest) :- 739 ( get_attr(Dep,itf,Att), 740 arg(2,Att,type(Type)), 741 arg(4,Att,lin(Lin)), 742 lb_inner(Type,OrdX,Lin,W,Lb), 743 Lb =< 0 % Lb > 0 means a violated bound 744 -> lb(Deps,OrdX,Dep-W-Lb,Tightest) 745 ; lb_first(Deps,OrdX,Tightest) 746 ). 747 748% lb(Deps,OrdX,TightestIn,TightestOut) 749% 750% See lb_first/3: this one does the same thing, but is used for the steps after 751% the first one and remembers the tightest bound so far. 752 753lb([],_,T0,T0). 754lb([Dep|Deps],OrdX,T0,T1) :- 755 ( get_attr(Dep,itf,Att), 756 arg(2,Att,type(Type)), 757 arg(4,Att,lin(Lin)), 758 lb_inner(Type,OrdX,Lin,W,Lb), 759 T0 = _-Lbb, 760 Lb > Lbb, % choose the least lowering, others might violate 761 % bounds 762 Lb =< 0 % violation of a bound (without lowering) 763 -> lb(Deps,OrdX,Dep-W-Lb,T1) 764 ; lb(Deps,OrdX,T0,T1) 765 ). 766 767% lb_inner(Type,X,Lin,W,Lb) 768% 769% Returns in Lb how much lower we can make X without violating a bound 770% by using the linear equation Lin of basis variable B which has type 771% Type and which has to activate a bound (type W) to do so. 772% 773% E.g. when B has a lowerbound L, then L should always be smaller than I + R. 774% So a lowerbound of X (which has scalar K in Lin), could be at most 775% (L-(I+R))/K lower than its upperbound (if K is positive). 776% Also note that Lb should always be smaller than 0, otherwise the row is 777% not feasible. 778% X has ordering attribute OrdX. 779 780lb_inner(t_l(L),OrdX,Lin,t_L(L),Lb) :- 781 nf_rhs_x(Lin,OrdX,Rhs,K), % if linear equation Lin contains the term 782 % X*K, Rhs is the right hand side of that 783 % equation 784 K > 0, 785 Lb is (L - Rhs) rdiv K. 786lb_inner(t_u(U),OrdX,Lin,t_U(U),Lb) :- 787 nf_rhs_x(Lin,OrdX,Rhs,K), 788 K < 0, % K < 0 789 Lb is (U - Rhs) rdiv K. 790lb_inner(t_lu(L,U),OrdX,Lin,W,Lb) :- 791 nf_rhs_x(Lin,OrdX,Rhs,K), 792 ( K < 0 793 -> W = t_lU(L,U), 794 Lb is (U - Rhs) rdiv K 795 ; K > 0 796 -> W = t_Lu(L,U), 797 Lb is (L - Rhs) rdiv K 798 ). 799 800% ---------------------------------- equations -------------------------------- 801% 802% backsubstitution will not make the system infeasible, if the bounds on the 803% indep vars are obeyed, but some implied values might pop up in rows where X 804% occurs 805% -) special case X=Y during bs -> get rid of dependend var(s), alias 806% 807 808solve(Lin) :- 809 Lin = [I,_|H], 810 solve(H,Lin,I,Bindings,[]), 811 export_binding(Bindings). 812 813% solve(Hom,Lin,I,Bind,BindT) 814% 815% Solves a linear equation Lin = [I,_|H] = 0 and exports the generated bindings 816 817solve([],_,I,Bind0,Bind0) :- 818 !, 819 I =:= 0. 820solve(H,Lin,_,Bind0,BindT) :- 821 sd(H,[],ClassesUniq,9-9-0,Category-Selected-_,NV,NVT), 822 get_attr(Selected,itf,Att), 823 arg(5,Att,order(Ord)), 824 isolate(Ord,Lin,Lin1), % Lin = 0 => Selected = Lin1 825 ( Category = 1 % classless variable, no bounds 826 -> setarg(4,Att,lin(Lin1)), 827 Lin1 = [Inhom,_|Hom], 828 bs_collect_binding(Hom,Selected,Inhom,Bind0,BindT), 829 eq_classes(NV,NVT,ClassesUniq) 830 ; Category = 2 % class variable, no bounds 831 -> arg(6,Att,class(NewC)), 832 class_allvars(NewC,Deps), 833 ( ClassesUniq = [_] % rank increasing 834 -> bs_collect_bindings(Deps,Ord,Lin1,Bind0,BindT) 835 ; Bind0 = BindT, 836 bs(Deps,Ord,Lin1) 837 ), 838 eq_classes(NV,NVT,ClassesUniq) 839 ; Category = 3 % classless variable, all variables in Lin and 840 % Selected are bounded 841 -> arg(2,Att,type(Type)), 842 setarg(4,Att,lin(Lin1)), 843 deactivate_bound(Type,Selected), 844 eq_classes(NV,NVT,ClassesUniq), 845 basis_add(Selected,Basis), 846 undet_active(Lin1), % we can't tell which bound will likely be a 847 % problem at this point 848 Lin1 = [Inhom,_|Hom], 849 bs_collect_binding(Hom,Selected,Inhom,Bind0,Bind1), % only if 850 % Hom = [] 851 rcbl(Basis,Bind1,BindT) % reconsider entire basis 852 ; Category = 4 % class variable, all variables in Lin and Selected 853 % are bounded 854 -> arg(2,Att,type(Type)), 855 arg(6,Att,class(NewC)), 856 class_allvars(NewC,Deps), 857 ( ClassesUniq = [_] % rank increasing 858 -> bs_collect_bindings(Deps,Ord,Lin1,Bind0,Bind1) 859 ; Bind0 = Bind1, 860 bs(Deps,Ord,Lin1) 861 ), 862 deactivate_bound(Type,Selected), 863 basis_add(Selected,Basis), 864 % eq_classes( NV, NVT, ClassesUniq), 865 % 4 -> var(NV) 866 equate(ClassesUniq,_), 867 undet_active(Lin1), 868 rcbl(Basis,Bind1,BindT) 869 ). 870 871% 872% Much like solve, but we solve for a particular variable of type t_none 873% 874 875% solve_x(H,Lin,I,X,[Bind|BindT],BindT) 876% 877% 878 879solve_x(Lin,X) :- 880 Lin = [I,_|H], 881 solve_x(H,Lin,I,X,Bindings,[]), 882 export_binding(Bindings). 883 884solve_x([],_,I,_,Bind0,Bind0) :- 885 !, 886 I =:= 0. 887solve_x(H,Lin,_,X,Bind0,BindT) :- 888 sd(H,[],ClassesUniq,9-9-0,_,NV,NVT), 889 get_attr(X,itf,Att), 890 arg(5,Att,order(OrdX)), 891 isolate(OrdX,Lin,Lin1), 892 ( arg(6,Att,class(NewC)) 893 -> class_allvars(NewC,Deps), 894 ( ClassesUniq = [_] % rank increasing 895 -> bs_collect_bindings(Deps,OrdX,Lin1,Bind0,BindT) 896 ; Bind0 = BindT, 897 bs(Deps,OrdX,Lin1) 898 ), 899 eq_classes(NV,NVT,ClassesUniq) 900 ; setarg(4,Att,lin(Lin1)), 901 Lin1 = [Inhom,_|Hom], 902 bs_collect_binding(Hom,X,Inhom,Bind0,BindT), 903 eq_classes(NV,NVT,ClassesUniq) 904 ). 905 906% solve_ord_x(Lin,OrdX,ClassX) 907% 908% Does the same thing as solve_x/2, but has the ordering of X and its class as 909% input. This also means that X has a class which is not sure in solve_x/2. 910 911solve_ord_x(Lin,OrdX,ClassX) :- 912 Lin = [I,_|H], 913 solve_ord_x(H,Lin,I,OrdX,ClassX,Bindings,[]), 914 export_binding(Bindings). 915 916solve_ord_x([],_,I,_,_,Bind0,Bind0) :- 917 I =:= 0. 918solve_ord_x([_|_],Lin,_,OrdX,ClassX,Bind0,BindT) :- 919 isolate(OrdX,Lin,Lin1), 920 Lin1 = [_,_|H1], 921 sd(H1,[],ClassesUniq1,9-9-0,_,NV,NVT), % do sd on Lin without X, then 922 % add class of X 923 ord_add_element(ClassesUniq1,ClassX,ClassesUniq), 924 class_allvars(ClassX,Deps), 925 ( ClassesUniq = [_] % rank increasing 926 -> bs_collect_bindings(Deps,OrdX,Lin1,Bind0,BindT) 927 ; Bind0 = BindT, 928 bs(Deps,OrdX,Lin1) 929 ), 930 eq_classes(NV,NVT,ClassesUniq). 931 932% sd(H,[],ClassesUniq,9-9-0,Category-Selected-_,NV,NVT) 933 934% sd(Hom,ClassesIn,ClassesOut,PreferenceIn,PreferenceOut,[NV|NVTail],NVTail) 935% 936% ClassesOut is a sorted list of the different classes that are either in 937% ClassesIn or that are the classes of the variables in Hom. Variables that do 938% not belong to a class yet, are put in the difference list NV. 939 940sd([],Class0,Class0,Preference0,Preference0,NV0,NV0). 941sd([l(X*K,_)|Xs],Class0,ClassN,Preference0,PreferenceN,NV0,NVt) :- 942 get_attr(X,itf,Att), 943 ( arg(6,Att,class(Xc)) % old: has class 944 -> NV0 = NV1, 945 ord_add_element(Class0,Xc,Class1), 946 ( arg(2,Att,type(t_none)) 947 -> preference(Preference0,2-X-K,Preference1) 948 % has class, no bounds => category 2 949 ; preference(Preference0,4-X-K,Preference1) 950 % has class, is bounded => category 4 951 ) 952 ; % new: has no class 953 Class1 = Class0, 954 NV0 = [X|NV1], % X has no class yet, add to list of new variables 955 ( arg(2,Att,type(t_none)) 956 -> preference(Preference0,1-X-K,Preference1) 957 % no class, no bounds => category 1 958 ; preference(Preference0,3-X-K,Preference1) 959 % no class, is bounded => category 3 960 ) 961 ), 962 sd(Xs,Class1,ClassN,Preference1,PreferenceN,NV1,NVt). 963 964% 965% A is best sofar, B is current 966% smallest prefered 967preference(A,B,Pref) :- 968 A = Px-_-_, 969 B = Py-_-_, 970 ( Px < Py 971 -> Pref = A 972 ; Pref = B 973 ). 974 975% eq_classes(NV,NVTail,Cs) 976% 977% Attaches all classless variables NV to a new class and equates all other 978% classes with this class. The equate operation only happens after attach_class 979% because the unification of classes can bind the tail of the AllVars attribute 980% to a nonvar and then the attach_class operation wouldn't work. 981 982eq_classes(NV,_,Cs) :- 983 var(NV), 984 !, 985 equate(Cs,_). 986eq_classes(NV,NVT,Cs) :- 987 class_new(Su,clpq,NV,NVT,[]), % make a new class Su with NV as the variables 988 attach_class(NV,Su), % attach the variables NV to Su 989 equate(Cs,Su). 990 991equate([],_). 992equate([X|Xs],X) :- equate(Xs,X). 993 994% 995% assert: none of the Vars has a class attribute yet 996% 997attach_class(Xs,_) :- 998 var(Xs), % Tail 999 !. 1000attach_class([X|Xs],Class) :- 1001 get_attr(X,itf,Att), 1002 setarg(6,Att,class(Class)), 1003 attach_class(Xs,Class). 1004 1005% unconstrained(Lin,Uc,Kuc,Rest) 1006% 1007% Finds an unconstrained variable Uc (type(t_none)) in Lin with scalar Kuc and 1008% removes it from Lin to return Rest. 1009 1010unconstrained(Lin,Uc,Kuc,Rest) :- 1011 Lin = [_,_|H], 1012 sd(H,[],_,9-9-0,Category-Uc-_,_,_), 1013 Category =< 2, 1014 get_attr(Uc,itf,Att), 1015 arg(5,Att,order(OrdUc)), 1016 delete_factor(OrdUc,Lin,Rest,Kuc). 1017 1018% 1019% point the vars in Lin into the same equivalence class 1020% maybe join some global data 1021% 1022same_class([],_). 1023same_class([l(X*_,_)|Xs],Class) :- 1024 get_or_add_class(X,Class), 1025 same_class(Xs,Class). 1026 1027% get_or_add_class(X,Class) 1028% 1029% Returns in Class the class of X if X has one, or a new class where X now 1030% belongs to if X didn't have one. 1031 1032get_or_add_class(X,Class) :- 1033 get_attr(X,itf,Att), 1034 arg(1,Att,CLP), 1035 ( arg(6,Att,class(ClassX)) 1036 -> ClassX = Class 1037 ; setarg(6,Att,class(Class)), 1038 class_new(Class,CLP,[X|Tail],Tail,[]) 1039 ). 1040 1041% allvars(X,Allvars) 1042% 1043% Allvars is a list of all variables in the class to which X belongs. 1044 1045allvars(X,Allvars) :- 1046 get_attr(X,itf,Att), 1047 arg(6,Att,class(C)), 1048 class_allvars(C,Allvars). 1049 1050% deactivate_bound(Type,Variable) 1051% 1052% The Type of the variable is changed to reflect the deactivation of its 1053% bounds. 1054% t_L(_) becomes t_l(_), t_lU(_,_) becomes t_lu(_,_) and so on. 1055 1056deactivate_bound(t_l(_),_). 1057deactivate_bound(t_u(_),_). 1058deactivate_bound(t_lu(_,_),_). 1059deactivate_bound(t_L(L),X) :- 1060 get_attr(X,itf,Att), 1061 setarg(2,Att,type(t_l(L))). 1062deactivate_bound(t_Lu(L,U),X) :- 1063 get_attr(X,itf,Att), 1064 setarg(2,Att,type(t_lu(L,U))). 1065deactivate_bound(t_U(U),X) :- 1066 get_attr(X,itf,Att), 1067 setarg(2,Att,type(t_u(U))). 1068deactivate_bound(t_lU(L,U),X) :- 1069 get_attr(X,itf,Att), 1070 setarg(2,Att,type(t_lu(L,U))). 1071 1072% intro_at(X,Value,Type) 1073% 1074% Variable X gets new type Type which reflects the activation of a bound with 1075% value Value. In the linear equations of all the variables belonging to the 1076% same class as X, X is substituted by [0,Value,X] to reflect the new active 1077% bound. 1078 1079intro_at(X,Value,Type) :- 1080 get_attr(X,itf,Att), 1081 arg(5,Att,order(Ord)), 1082 arg(6,Att,class(Class)), 1083 setarg(2,Att,type(Type)), 1084 ( Value =:= 0 1085 -> true 1086 ; backsubst_delta(Class,Ord,X,Value) 1087 ). 1088 1089% undet_active(Lin) 1090% 1091% For each variable in the homogene part of Lin, a bound is activated 1092% if an inactive bound exists. (t_l(L) becomes t_L(L) and so on) 1093 1094undet_active([_,_|H]) :- 1095 undet_active_h(H). 1096 1097% undet_active_h(Hom) 1098% 1099% For each variable in homogene part Hom, a bound is activated if an 1100% inactive bound exists (t_l(L) becomes t_L(L) and so on) 1101 1102undet_active_h([]). 1103undet_active_h([l(X*_,_)|Xs]) :- 1104 get_attr(X,itf,Att), 1105 arg(2,Att,type(Type)), 1106 undet_active(Type,X), 1107 undet_active_h(Xs). 1108 1109% undet_active(Type,Var) 1110% 1111% An inactive bound of Var is activated if such exists 1112% t_lu(L,U) is arbitrarily chosen to become t_Lu(L,U) 1113 1114undet_active(t_none,_). % type_activity 1115undet_active(t_L(_),_). 1116undet_active(t_Lu(_,_),_). 1117undet_active(t_U(_),_). 1118undet_active(t_lU(_,_),_). 1119undet_active(t_l(L),X) :- intro_at(X,L,t_L(L)). 1120undet_active(t_u(U),X) :- intro_at(X,U,t_U(U)). 1121undet_active(t_lu(L,U),X) :- intro_at(X,L,t_Lu(L,U)). 1122 1123% determine_active_dec(Lin) 1124% 1125% Activates inactive bounds on the variables of Lin if such bounds exist. 1126% If the type of a variable is t_none, this fails. This version is aimed 1127% to make the R component of Lin as small as possible in order not to violate 1128% an upperbound (see reconsider/1) 1129 1130determine_active_dec([_,_|H]) :- 1131 determine_active(H,-1). 1132 1133% determine_active_inc(Lin) 1134% 1135% Activates inactive bounds on the variables of Lin if such bounds exist. 1136% If the type of a variable is t_none, this fails. This version is aimed 1137% to make the R component of Lin as large as possible in order not to violate 1138% a lowerbound (see reconsider/1) 1139 1140determine_active_inc([_,_|H]) :- 1141 determine_active(H,1). 1142 1143% determine_active(Hom,S) 1144% 1145% For each variable in Hom, activates its bound if it is not yet activated. 1146% For the case of t_lu(_,_) the lower or upper bound is activated depending on 1147% K and S: 1148% If sign of K*S is negative, then lowerbound, otherwise upperbound. 1149 1150determine_active([],_). 1151determine_active([l(X*K,_)|Xs],S) :- 1152 get_attr(X,itf,Att), 1153 arg(2,Att,type(Type)), 1154 determine_active(Type,X,K,S), 1155 determine_active(Xs,S). 1156 1157determine_active(t_L(_),_,_,_). 1158determine_active(t_Lu(_,_),_,_,_). 1159determine_active(t_U(_),_,_,_). 1160determine_active(t_lU(_,_),_,_,_). 1161determine_active(t_l(L),X,_,_) :- intro_at(X,L,t_L(L)). 1162determine_active(t_u(U),X,_,_) :- intro_at(X,U,t_U(U)). 1163determine_active(t_lu(L,U),X,K,S) :- 1164 KS is K*S, 1165 ( KS < 0 1166 -> intro_at(X,L,t_Lu(L,U)) 1167 ; KS > 0 1168 -> intro_at(X,U,t_lU(L,U)) 1169 ). 1170 1171% 1172% Careful when an indep turns into t_none !!! 1173% 1174 1175detach_bounds(V) :- 1176 get_attr(V,itf,Att), 1177 arg(2,Att,type(Type)), 1178 arg(4,Att,lin(Lin)), 1179 arg(5,Att,order(OrdV)), 1180 arg(6,Att,class(Class)), 1181 setarg(2,Att,type(t_none)), 1182 setarg(3,Att,strictness(0)), 1183 ( indep(Lin,OrdV) 1184 -> ( ub(Class,OrdV,Vub-Vb-_) 1185 -> % exchange against thightest 1186 class_basis_drop(Class,Vub), 1187 pivot(Vub,Class,OrdV,Vb,Type) 1188 ; lb(Class,OrdV,Vlb-Vb-_) 1189 -> class_basis_drop(Class,Vlb), 1190 pivot(Vlb,Class,OrdV,Vb,Type) 1191 ; true 1192 ) 1193 ; class_basis_drop(Class,V) 1194 ). 1195 1196detach_bounds_vlv(OrdV,Lin,Class,Var,NewLin) :- 1197 ( indep(Lin,OrdV) 1198 -> Lin = [_,R|_], 1199 ( ub(Class,OrdV,Vub-Vb-_) 1200 -> % in verify_lin, class might contain two occurrences of Var, 1201 % but it doesn't matter which one we delete 1202 class_basis_drop(Class,Var), 1203 pivot_vlv(Vub,Class,OrdV,Vb,R,NewLin) 1204 ; lb(Class,OrdV,Vlb-Vb-_) 1205 -> class_basis_drop(Class,Var), 1206 pivot_vlv(Vlb,Class,OrdV,Vb,R,NewLin) 1207 ; NewLin = Lin 1208 ) 1209 ; NewLin = Lin, 1210 class_basis_drop(Class,Var) 1211 ). 1212 1213% ----------------------------- manipulate the basis -------------------------- 1214 1215% basis_drop(X) 1216% 1217% Removes X from the basis of the class to which X belongs. 1218 1219basis_drop(X) :- 1220 get_attr(X,itf,Att), 1221 arg(6,Att,class(Cv)), 1222 class_basis_drop(Cv,X). 1223 1224% basis(X,Basis) 1225% 1226% Basis is the basis of the class to which X belongs. 1227 1228basis(X,Basis) :- 1229 get_attr(X,itf,Att), 1230 arg(6,Att,class(Cv)), 1231 class_basis(Cv,Basis). 1232 1233% basis_add(X,NewBasis) 1234% 1235% NewBasis is the result of adding X to the basis of the class to which X 1236% belongs. 1237 1238basis_add(X,NewBasis) :- 1239 get_attr(X,itf,Att), 1240 arg(6,Att,class(Cv)), 1241 class_basis_add(Cv,X,NewBasis). 1242 1243% basis_pivot(Leave,Enter) 1244% 1245% Removes Leave from the basis of the class to which it belongs, and adds 1246% Enter to that basis. 1247 1248basis_pivot(Leave,Enter) :- 1249 get_attr(Leave,itf,Att), 1250 arg(6,Att,class(Cv)), 1251 class_basis_pivot(Cv,Enter,Leave). 1252 1253% ----------------------------------- pivot ----------------------------------- 1254 1255% pivot(Dep,Indep) 1256% 1257% The linear equation of variable Dep, is transformed into one of variable 1258% Indep, containing Dep. Then, all occurrences of Indep in linear equations are 1259% substituted by this new definition 1260 1261% 1262% Pivot ignoring rhs and active states 1263% 1264 1265pivot(Dep,Indep) :- 1266 get_attr(Dep,itf,AttD), 1267 arg(4,AttD,lin(H)), 1268 arg(5,AttD,order(OrdDep)), 1269 get_attr(Indep,itf,AttI), 1270 arg(5,AttI,order(Ord)), 1271 arg(5,AttI,class(Class)), 1272 delete_factor(Ord,H,H0,Coeff), 1273 K is -1 rdiv Coeff, 1274 add_linear_ff(H0,K,[0,0,l(Dep* -1,OrdDep)],K,Lin), 1275 backsubst(Class,Ord,Lin). 1276 1277% pivot_a(Dep,Indep,IndepT,DepT) 1278% 1279% Removes Dep from the basis, puts Indep in, and pivots the equation of 1280% Dep to become one of Indep. The type of Dep becomes DepT (which means 1281% it gets deactivated), the type of Indep becomes IndepT (which means it 1282% gets activated) 1283 1284 1285pivot_a(Dep,Indep,Vb,Wd) :- 1286 basis_pivot(Dep,Indep), 1287 get_attr(Indep,itf,Att), 1288 arg(2,Att,type(Type)), 1289 arg(5,Att,order(Ord)), 1290 arg(6,Att,class(Class)), 1291 pivot(Dep,Class,Ord,Vb,Type), 1292 get_attr(Indep,itf,Att2), %changed? 1293 setarg(2,Att2,type(Wd)). 1294 1295pivot_b(Vub,V,Vb,Wd) :- 1296 ( Vub == V 1297 -> get_attr(V,itf,Att), 1298 arg(5,Att,order(Ord)), 1299 arg(6,Att,class(Class)), 1300 setarg(2,Att,type(Vb)), 1301 pivot_b_delta(Vb,Delta), % nonzero(Delta) 1302 backsubst_delta(Class,Ord,V,Delta) 1303 ; pivot_a(Vub,V,Vb,Wd) 1304 ). 1305 1306pivot_b_delta(t_Lu(L,U),Delta) :- Delta is L-U. 1307pivot_b_delta(t_lU(L,U),Delta) :- Delta is U-L. 1308 1309% select_active_bound(Type,Bound) 1310% 1311% Returns the bound that is active in Type (if such exists, 0 otherwise) 1312 1313select_active_bound(t_L(L),L). 1314select_active_bound(t_Lu(L,_),L). 1315select_active_bound(t_U(U),U). 1316select_active_bound(t_lU(_,U),U). 1317select_active_bound(t_none,0). 1318% 1319% for project.pl 1320% 1321select_active_bound(t_l(_),0). 1322select_active_bound(t_u(_),0). 1323select_active_bound(t_lu(_,_),0). 1324 1325 1326% pivot(Dep,Class,IndepOrd,DepAct,IndAct) 1327% 1328% See pivot/2. 1329% In addition, variable Indep with ordering IndepOrd has an active bound IndAct 1330 1331% 1332% 1333% Pivot taking care of rhs and active states 1334% 1335pivot(Dep,Class,IndepOrd,DepAct,IndAct) :- 1336 get_attr(Dep,itf,Att), 1337 arg(4,Att,lin(H)), 1338 arg(5,Att,order(DepOrd)), 1339 setarg(2,Att,type(DepAct)), 1340 select_active_bound(DepAct,AbvD), % New current value for Dep 1341 select_active_bound(IndAct,AbvI), % Old current value of Indep 1342 delete_factor(IndepOrd,H,H0,Coeff), % Dep = ... + Coeff*Indep + ... 1343 AbvDm is -AbvD, 1344 AbvIm is -AbvI, 1345 add_linear_f1([0,AbvIm],Coeff,H0,H1), 1346 K is -1 rdiv Coeff, 1347 add_linear_ff(H1,K,[0,AbvDm,l(Dep* -1,DepOrd)],K,H2), 1348 % Indep = -1/Coeff*... + 1/Coeff*Dep 1349 add_linear_11(H2,[0,AbvIm],Lin), 1350 backsubst(Class,IndepOrd,Lin). 1351 1352% Rewrite Dep = ... + Coeff*Indep + ... 1353% into Indep = ... + -1/Coeff*Dep + ... 1354% 1355% For backsubstitution, old current value of Indep must be removed from RHS 1356% New current value of Dep must be added to RHS 1357% For solving: old current value of Indep should be out of RHS 1358 1359pivot_vlv(Dep,Class,IndepOrd,DepAct,AbvI,Lin) :- 1360 get_attr(Dep,itf,Att), 1361 arg(4,Att,lin(H)), 1362 arg(5,Att,order(DepOrd)), 1363 setarg(2,Att,type(DepAct)), 1364 select_active_bound(DepAct,AbvD), % New current value for Dep 1365 delete_factor(IndepOrd,H,H0,Coeff), % Dep = ... + Coeff*Indep + ... 1366 AbvDm is -AbvD, 1367 AbvIm is -AbvI, 1368 add_linear_f1([0,AbvIm],Coeff,H0,H1), 1369 K is -1 rdiv Coeff, 1370 add_linear_ff(H1,K,[0,AbvDm,l(Dep* -1,DepOrd)],K,Lin), 1371 % Indep = -1/Coeff*... + 1/Coeff*Dep 1372 add_linear_11(Lin,[0,AbvIm],SubstLin), 1373 backsubst(Class,IndepOrd,SubstLin). 1374 1375% backsubst_delta(Class,OrdX,X,Delta) 1376% 1377% X with ordering attribute OrdX, is substituted in all linear equations of 1378% variables in the class Class, by linear equation [0,Delta,l(X*1,OrdX)]. This 1379% reflects the activation of a bound. 1380 1381backsubst_delta(Class,OrdX,X,Delta) :- 1382 backsubst(Class,OrdX,[0,Delta,l(X*1,OrdX)]). 1383 1384% backsubst(Class,OrdX,Lin) 1385% 1386% X with ordering OrdX is substituted in all linear equations of variables in 1387% the class Class, by linear equation Lin 1388 1389backsubst(Class,OrdX,Lin) :- 1390 class_allvars(Class,Allvars), 1391 bs(Allvars,OrdX,Lin). 1392 1393% bs(Vars,OrdV,Lin) 1394% 1395% In all linear equations of the variables Vars, variable V with ordering 1396% attribute OrdV is substituted by linear equation Lin. 1397% 1398% valid if nothing will go ground 1399% 1400 1401bs(Xs,_,_) :- 1402 var(Xs), 1403 !. 1404bs([X|Xs],OrdV,Lin) :- 1405 ( get_attr(X,itf,Att), 1406 arg(4,Att,lin(LinX)), 1407 nf_substitute(OrdV,Lin,LinX,LinX1) % does not change attributes 1408 -> setarg(4,Att,lin(LinX1)), 1409 bs(Xs,OrdV,Lin) 1410 ; bs(Xs,OrdV,Lin) 1411 ). 1412 1413% 1414% rank increasing backsubstitution 1415% 1416 1417% bs_collect_bindings(Deps,SelectedOrd,Lin,Bind,BindT) 1418% 1419% Collects bindings (of the form [X-I] where X = I is the binding) by 1420% substituting Selected in all linear equations of the variables Deps (which 1421% are of the same class), by Lin. Selected has ordering attribute SelectedOrd. 1422% 1423% E.g. when V = 2X + 3Y + 4, X = 3V + 2Z and Y = 4X + 3 1424% we can substitute V in the linear equation of X: X = 6X + 9Y + 2Z + 12 1425% we can't substitute V in the linear equation of Y of course. 1426 1427bs_collect_bindings(Xs,_,_,Bind0,BindT) :- 1428 var(Xs), 1429 !, 1430 Bind0 = BindT. 1431bs_collect_bindings([X|Xs],OrdV,Lin,Bind0,BindT) :- 1432 ( get_attr(X,itf,Att), 1433 arg(4,Att,lin(LinX)), 1434 nf_substitute(OrdV,Lin,LinX,LinX1) % does not change attributes 1435 -> setarg(4,Att,lin(LinX1)), 1436 LinX1 = [Inhom,_|Hom], 1437 bs_collect_binding(Hom,X,Inhom,Bind0,Bind1), 1438 bs_collect_bindings(Xs,OrdV,Lin,Bind1,BindT) 1439 ; bs_collect_bindings(Xs,OrdV,Lin,Bind0,BindT) 1440 ). 1441 1442% bs_collect_binding(Hom,Selected,Inhom,Bind,BindT) 1443% 1444% Collects binding following from Selected = Hom + Inhom. 1445% If Hom = [], returns the binding Selected-Inhom (=0) 1446% 1447bs_collect_binding([],X,Inhom) --> [X-Inhom]. 1448bs_collect_binding([_|_],_,_) --> []. 1449 1450% 1451% reconsider the basis 1452% 1453 1454% rcbl(Basis,Bind,BindT) 1455% 1456% 1457 1458rcbl([],Bind0,Bind0). 1459rcbl([X|Continuation],Bind0,BindT) :- 1460 ( rcb_cont(X,Status,Violated,Continuation,NewContinuation) % have a culprit 1461 -> rcbl_status(Status,X,NewContinuation,Bind0,BindT,Violated) 1462 ; rcbl(Continuation,Bind0,BindT) 1463 ). 1464 1465rcb_cont(X,Status,Violated,ContIn,ContOut) :- 1466 get_attr(X,itf,Att), 1467 arg(2,Att,type(Type)), 1468 arg(4,Att,lin([I,R|H])), 1469 ( Type = t_l(L) % case 1: lowerbound: R + I should always be larger 1470 % than the lowerbound 1471 -> R + I =< L, 1472 Violated = l(L), 1473 inc_step_cont(H,Status,ContIn,ContOut) 1474 ; Type = t_u(U) % case 2: upperbound: R + I should always be smaller 1475 % than the upperbound 1476 -> R + I >= U, 1477 Violated = u(U), 1478 dec_step_cont(H,Status,ContIn,ContOut) 1479 ; Type = t_lu(L,U) % case 3: check both 1480 -> At is R + I, 1481 ( At =< L 1482 -> Violated = l(L), 1483 inc_step_cont(H,Status,ContIn,ContOut) 1484 ; At >= U 1485 -> Violated = u(U), 1486 dec_step_cont(H,Status,ContIn,ContOut) 1487 ) 1488 ). % other types imply nonbasic variable or unbounded variable 1489 1490 1491 1492% 1493% reconsider one element of the basis 1494% later: lift the binds 1495% 1496reconsider(X) :- 1497 rcb(X,Status,Violated), 1498 !, 1499 rcbl_status(Status,X,[],Binds,[],Violated), 1500 export_binding(Binds). 1501reconsider(_). 1502 1503% 1504% Find a basis variable out of its bound or at its bound 1505% Try to move it into whithin its bound 1506% a) impossible -> fail 1507% b) optimum at the bound -> implied value 1508% c) else look at the remaining basis variables 1509% 1510% 1511% Idea: consider a variable V with linear equation Lin. 1512% When a bound on a variable X of Lin gets activated, its value, multiplied 1513% with the scalar of X, is added to the R component of Lin. 1514% When we consider the lowerbound of V, it must be smaller than R + I, since R 1515% contains at best the lowerbounds of the variables in Lin (but could contain 1516% upperbounds, which are of course larger). So checking this can show the 1517% violation of a bound of V. A similar case works for the upperbound. 1518 1519rcb(X,Status,Violated) :- 1520 get_attr(X,itf,Att), 1521 arg(2,Att,type(Type)), 1522 arg(4,Att,lin([I,R|H])), 1523 ( Type = t_l(L) % case 1: lowerbound: R + I should always be larger 1524 % than the lowerbound 1525 -> R + I =< L, 1526 Violated = l(L), 1527 inc_step(H,Status) 1528 ; Type = t_u(U) % case 2: upperbound: R + I should always be smaller 1529 % than the upperbound 1530 -> R + I >= U, 1531 Violated = u(U), 1532 dec_step(H,Status) 1533 ; Type = t_lu(L,U) % case 3: check both 1534 -> At is R + I, 1535 ( At =< L 1536 -> Violated = l(L), 1537 inc_step(H,Status) 1538 ; At >= U 1539 -> Violated = u(U), 1540 dec_step(H,Status) 1541 ) 1542 ). % other types imply nonbasic variable or unbounded variable 1543 1544% rcbl_status(Status,X,Continuation,[Bind|BindT],BindT,Violated) 1545% 1546% 1547 1548rcbl_status(optimum,X,Cont,B0,Bt,Violated) :- rcbl_opt(Violated,X,Cont,B0,Bt). 1549rcbl_status(applied,X,Cont,B0,Bt,Violated) :- rcbl_app(Violated,X,Cont,B0,Bt). 1550rcbl_status(unlimited(Indep,DepT),X,Cont,B0,Bt,Violated) :- 1551 rcbl_unl(Violated,X,Cont,B0,Bt,Indep,DepT). 1552 1553% 1554% Might reach optimum immediately without changing the basis, 1555% but in general we must assume that there were pivots. 1556% If the optimum meets the bound, we backsubstitute the implied 1557% value, solve will call us again to check for further implied 1558% values or unsatisfiability in the rank increased system. 1559% 1560rcbl_opt(l(L),X,Continuation,B0,B1) :- 1561 get_attr(X,itf,Att), 1562 arg(2,Att,type(Type)), 1563 arg(3,Att,strictness(Strict)), 1564 arg(4,Att,lin(Lin)), 1565 Lin = [I,R|_], 1566 Opt is R + I, 1567 ( L < Opt 1568 -> narrow_u(Type,X,Opt), % { X =< Opt } 1569 rcbl(Continuation,B0,B1) 1570 ; L =:= Opt, 1571 Strict /\ 2 =:= 0, % meets lower 1572 Mop is -Opt, 1573 normalize_scalar(Mop,MopN), 1574 add_linear_11(MopN,Lin,Lin1), 1575 Lin1 = [Inhom,_|Hom], 1576 ( Hom = [] 1577 -> rcbl(Continuation,B0,B1) % would not callback 1578 ; solve(Hom,Lin1,Inhom,B0,B1) 1579 ) 1580 ). 1581rcbl_opt(u(U),X,Continuation,B0,B1) :- 1582 get_attr(X,itf,Att), 1583 arg(2,Att,type(Type)), 1584 arg(3,Att,strictness(Strict)), 1585 arg(4,Att,lin(Lin)), 1586 Lin = [I,R|_], 1587 Opt is R + I, 1588 ( U > Opt 1589 -> narrow_l(Type,X,Opt), % { X >= Opt } 1590 rcbl(Continuation,B0,B1) 1591 ; U =:= Opt, 1592 Strict /\ 1 =:= 0, % meets upper 1593 Mop is -Opt, 1594 normalize_scalar(Mop,MopN), 1595 add_linear_11(MopN,Lin,Lin1), 1596 Lin1 = [Inhom,_|Hom], 1597 ( Hom = [] 1598 -> rcbl(Continuation,B0,B1) % would not callback 1599 ; solve(Hom,Lin1,Inhom,B0,B1) 1600 ) 1601 ). 1602 1603% 1604% Basis has already changed when this is called 1605% 1606rcbl_app(l(L),X,Continuation,B0,B1) :- 1607 get_attr(X,itf,Att), 1608 arg(4,Att,lin([I,R|H])), 1609 ( R + I > L % within bound now 1610 -> rcbl(Continuation,B0,B1) 1611 ; inc_step(H,Status), 1612 rcbl_status(Status,X,Continuation,B0,B1,l(L)) 1613 ). 1614rcbl_app(u(U),X,Continuation,B0,B1) :- 1615 get_attr(X,itf,Att), 1616 arg(4,Att,lin([I,R|H])), 1617 ( R + I < U % within bound now 1618 -> rcbl(Continuation,B0,B1) 1619 ; dec_step(H,Status), 1620 rcbl_status(Status,X,Continuation,B0,B1,u(U)) 1621 ). 1622% 1623% This is never called for a t_lu culprit 1624% 1625rcbl_unl(l(L),X,Continuation,B0,B1,Indep,DepT) :- 1626 pivot_a(X,Indep,t_L(L),DepT), % changes the basis 1627 rcbl(Continuation,B0,B1). 1628rcbl_unl(u(U),X,Continuation,B0,B1,Indep,DepT) :- 1629 pivot_a(X,Indep,t_U(U),DepT), % changes the basis 1630 rcbl(Continuation,B0,B1). 1631 1632% narrow_u(Type,X,U) 1633% 1634% Narrows down the upperbound of X (type Type) to U. 1635% Fails if Type is not t_u(_) or t_lu(_) 1636 1637narrow_u(t_u(_),X,U) :- 1638 get_attr(X,itf,Att), 1639 setarg(2,Att,type(t_u(U))). 1640narrow_u(t_lu(L,_),X,U) :- 1641 get_attr(X,itf,Att), 1642 setarg(2,Att,type(t_lu(L,U))). 1643 1644% narrow_l(Type,X,L) 1645% 1646% Narrows down the lowerbound of X (type Type) to L. 1647% Fails if Type is not t_l(_) or t_lu(_) 1648 1649narrow_l( t_l(_), X, L) :- 1650 get_attr(X,itf,Att), 1651 setarg(2,Att,type(t_l(L))). 1652 1653narrow_l( t_lu(_,U), X, L) :- 1654 get_attr(X,itf,Att), 1655 setarg(2,Att,type(t_lu(L,U))). 1656 1657% ----------------------------------- dump ------------------------------------ 1658 1659% dump_var(Type,Var,I,H,Dump,DumpTail) 1660% 1661% Returns in Dump a representation of the linear constraint on variable 1662% Var which has linear equation H + I and has type Type. 1663 1664dump_var(t_none,V,I,H) --> 1665 !, 1666 ( { 1667 H = [l(W*K,_)], 1668 V == W, 1669 I =:= 0, 1670 K =:= 1 1671 } 1672 -> % indep var 1673 [] 1674 ; {nf2sum(H,I,Sum)}, 1675 [V = Sum] 1676 ). 1677dump_var(t_L(L),V,I,H) --> 1678 !, 1679 dump_var(t_l(L),V,I,H). 1680% case lowerbound: V >= L or V > L 1681% say V >= L, and V = K*V1 + ... + I, then K*V1 + ... + I >= L 1682% and K*V1 + ... >= L-I and V1 + .../K = (L-I)/K 1683dump_var(t_l(L),V,I,H) --> 1684 !, 1685 { 1686 H = [l(_*K,_)|_], % avoid 1 >= 0 1687 get_attr(V,itf,Att), 1688 arg(3,Att,strictness(Strict)), 1689 Sm is Strict /\ 2, 1690 Kr is 1 rdiv K, 1691 Li is Kr*(L - I), 1692 mult_hom(H,Kr,H1), 1693 nf2sum(H1,0,Sum), 1694 ( K > 0 % K > 0 1695 -> dump_strict(Sm,Sum >= Li,Sum > Li,Result) 1696 ; dump_strict(Sm,Sum =< Li,Sum < Li,Result) 1697 ) 1698 }, 1699 [Result]. 1700dump_var(t_U(U),V,I,H) --> 1701 !, 1702 dump_var(t_u(U),V,I,H). 1703dump_var(t_u(U),V,I,H) --> 1704 !, 1705 { 1706 H = [l(_*K,_)|_], % avoid 0 =< 1 1707 get_attr(V,itf,Att), 1708 arg(3,Att,strictness(Strict)), 1709 Sm is Strict /\ 1, 1710 Kr is 1 rdiv K, 1711 Ui is Kr*(U-I), 1712 mult_hom(H,Kr,H1), 1713 nf2sum(H1,0.0,Sum), 1714 ( K > 0 1715 -> dump_strict(Sm,Sum =< Ui,Sum < Ui,Result) 1716 ; dump_strict(Sm,Sum >= Ui,Sum > Ui,Result) 1717 ) 1718 }, 1719 [Result]. 1720dump_var(t_Lu(L,U),V,I,H) --> 1721 !, 1722 dump_var(t_l(L),V,I,H), 1723 dump_var(t_u(U),V,I,H). 1724dump_var(t_lU(L,U),V,I,H) --> 1725 !, 1726 dump_var(t_l(L),V,I,H), 1727 dump_var(t_u(U),V,I,H). 1728dump_var(t_lu(L,U),V,I,H) --> 1729 !, 1730 dump_var(t_l(L),V,I,H), 1731 dump_var(t_U(U),V,I,H). 1732dump_var(T,V,I,H) --> % should not happen 1733 [V:T:I+H]. 1734 1735% dump_strict(FilteredStrictness,Nonstrict,Strict,Res) 1736% 1737% Unifies Res with either Nonstrict or Strict depending on FilteredStrictness. 1738% FilteredStrictness is the component of strictness related to the bound: 0 1739% means nonstrict, 1 means strict upperbound, 2 means strict lowerbound, 1740% 3 is filtered out to either 1 or 2. 1741 1742dump_strict(0,Result,_,Result). 1743dump_strict(1,_,Result,Result). 1744dump_strict(2,_,Result,Result). 1745 1746% dump_nz(V,H,I,Dump,DumpTail) 1747% 1748% Returns in Dump a representation of the nonzero constraint of variable V 1749% which has linear 1750% equation H + I. 1751 1752dump_nz(_,H,I) --> 1753 { 1754 H = [l(_*K,_)|_], 1755 Kr is 1 rdiv K, 1756 I1 is -Kr*I, 1757 mult_hom(H,Kr,H1), 1758 nf2sum(H1,0,Sum) 1759 }, 1760 [Sum =\= I1]. 1761