1% BEGIN LICENSE BLOCK 2% Version: CMPL 1.1 3% 4% The contents of this file are subject to the Cisco-style Mozilla Public 5% License Version 1.1 (the "License"); you may not use this file except 6% in compliance with the License. You may obtain a copy of the License 7% at www.eclipse-clp.org/license. 8% 9% Software distributed under the License is distributed on an "AS IS" 10% basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See 11% the License for the specific language governing rights and limitations 12% under the License. 13% 14% The Original Code is The ECLiPSe Constraint Logic Programming System. 15% The Initial Developer of the Original Code is Cisco Systems, Inc. 16% Portions created by the Initial Developer are 17% Copyright (C) 2006 Cisco Systems, Inc. All Rights Reserved. 18% 19% Contributor(s): 20% 21% END LICENSE BLOCK 22 23% TODO: 24% - many integer 0 are used, which should probably be 0.0 25% - replace +/2 expressions by lin/1 terms 26 27:- lib(constraint_pools). 28:- lib(linearize). 29:- lib(hash). 30:- lib(eplex). 31:- lib(dual_var). 32:- lib(bfs). 33 34:- import 35 lp_var_non_monotonic_set_bounds/4 36 from eplex_. 37 38:- inline(verify/1, t_verify/2). 39t_verify(verify(Goal), (Goal->true;writeln(error,verify_failed(Goal)),abort)). 40verify(Goal) :- verify(Goal). 41 42% ---------------------------------------------------------------------- 43% 44% Meta-attribute related stuff 45% This is the attribute of the generated MP variables 46% ---------------------------------------------------------------------- 47 48:- export cg_var_print/2. 49 50:- meta_attribute(colgen, [print:cg_var_print/2, unify:unify_colgen/2, 51 get_bounds:cg_get_bounds/3, 52 set_bounds:cg_set_bounds/3]). 53 54:- export struct( 55 colgen( 56 mp_val, % float (sometimes int 0?) 57 cost, % float: MP objective coefficient 58 coeffs, % [cstrid-coeff], keysorted 59 aux, % 'artificial', 'stabilisation', [], 60 % or defined by subproblem sp_sol{} 61 lo, % float, possibly int?? 62 hi, % float, possibly int?? 63 type, % 'real' or 'integer' 64 master_prob, % suspension list 65 solver, % cg_prob{} 66 next % colgen{} or [] 67 ) 68 ). 69 70% ---------------------------------------------------------------------- 71% 72% Pools 73% 74% ---------------------------------------------------------------------- 75 76:- export colgen_instance/1. % colgen_instance(+PoolName) 77 78% ---------------------------------------------------------------------- 79% 80% Predicates with pool argments (don't reexport these in colgen!) 81% 82% ---------------------------------------------------------------------- 83 84:- export var_dual1/7. % var_dual/6 85:- export get_dual1/3. % get_dual/2 86:- export get_coeff1/3. % get_coeff/2 87:- export get_idx1/3. % get_idx/2 88:- export get_rhs1/3. % get_rhs/2 89:- export always_set_dual1/3. % always_set_dual/2 90:- export set_dual1/3. % set_dual/2 91 92:- export bp_solve1/2. % solve/1 93:- export cg_solver_setup/3. % solver_setup/2 94:- export cg_solver_setup/4. % solver_setup/3 95:- export cg_integers1/2. % integers/1 96:- export add_cg_pool_constraint/3. % identified_constraint/2 97:- export cg_eq/3. % =:=/2 98:- export cg_ge/3. % >=/2 99:- export cg_le/3. % =</2 100:- export cg_sp_count1/2. % cg_subproblem_count/1 101:- export cg_sp_sol/2. % cg_subproblem_solution/1 102:- export cg_valid_columns1/2. % valid_columns/1 103:- export cg_sp_rc_sum/2. % cg_subproblem_rc_sum/1 104:- export cg_optimal_rc1/2. % optimal_rc/1 105:- export cg_minimize/5. % minimize/4 106:- export cg_minimize/4. % minimize/3 107:- export cg_var_get1/4. % var_get/3 108:- export cg_get1/3. % get/2 109:- export cg_branch1/2. % branch/1 110:- export cg_branch1/3. % branch/2 111 112% ---------------------------------------------------------------------- 113% 114% Predicates with handle argments 115% 116% ---------------------------------------------------------------------- 117 118:- reexport var_dual/7, % var_dual/6 119 get_dual/3, % get_dual/2 120 get_coeff/3, % get_coeff/2 121 get_idx/3, % get_idx/2 122 get_rhs/3, % get_rhs/2 123 always_set_dual/3, % always_set_dual/2 124 set_dual/3 % set_dual/2 125 from dual_var. 126 127:- export bp_solve/2. 128:- export cg_solver_setup/8. 129:- export cg_solver_setup/9. 130:- export cg_integers/2. % should this be a handle predicate? 131%:- export cg_add_constraints/3. 132:- export cg_sp_count/2. 133:- export cg_valid_columns/2. 134:- export cg_optimal_rc/2. 135:- export cg_minimize/9. 136:- export cg_minimize/8. 137:- export cg_var_get/4. 138:- export cg_get/3. 139:- export cg_branch/2. 140:- export cg_branch/3. 141 142% ---------------------------------------------------------------------- 143 144:- local struct(cg_constraint_type(mp_only, mp_sp, mp_branch)). 145 146:- export op(700, xfx, [$>=, $=, $=<, $::]). 147 148% ---------------------------------------------------------------------- 149% 150% Problem handle structure 151% 152% ---------------------------------------------------------------------- 153 154:- export struct( 155 cg_prob( 156 157 % Global State: 158 159 master_prob, % the master problem eplex handle 160 bfs_tree, % bfs_tree for branch-and-price 161 mp_susp, % suspension list containing the 162 % suspension for the MP solver 163 const_obj, % real: the constant part of 164 % the cost funciton 165 phase1_obj, % Expr: the artificial variable cost fn 166 % for phase 1 of two-phase 167 sp_solution_call, % the user supplied subproblem 168 % solution predicate 169 pool, % the associated constraint pool 170 tolerance, % float: tolerance for optimality 171 branch_tolerance, % float: tolerance for optimality 172 info_messages, % on, off: info message status 173 on_degeneracy, % stop, continue: should we halt when 174 % we find degeneracy, or continue 175 % (if so the sp solver is assumed 176 % to deal with it) 177 stabilisation, % off, 178 % on(BoundIter, BoundUpdate, 179 % CoeffIter, CoeffUpdate), 180 % stab_pred(UpdatePred, StoppingPred): 181 % the policy to perform basis 182 % stabilisation - if off then no 183 % stabilisation is performed, 184 % if on(...) then the default policy 185 % is used with var bounds/coefficients 186 % updated by BoundUpdate/CoeffUpdate 187 % after BoundIter/CoeffIter iterations 188 % resp, otherwise a user defined 189 % policy is employed and UpdatePred/ 190 % StoppingPred should be predicates 191 % that perform the updates and test 192 % for stopping conditions. 193 stab_terms, % [StabTerm1,...,StabTermM]: 194 % list of stabilisation terms 195 stab_iter_counts, % stab_counters struct: 196 % record of how many iterations 197 % since stabilisation var bound and 198 % coeff update 199 disallow, % lp, clp, off: 200 % policy for active prevention of 201 % duplicate columns 202 basis_perturbation,% on, off: 203 % should we try and perturb external 204 % solver basis when appear to be at 205 % optimal when external solver returns 206 % same basis after adding columns: 207 % off - no 208 % on - temporarily set the external 209 % solver to always perturb 210 upper_bound, % float: current bounds on solution 211 lower_bound, % objective value 212 integral_obj, % atom: yes or no, is the cost of all 213 % feasible solutions integral? (currently 214 % unused, but included for later) 215 duals, % array: array of current 216 % master problem dual values 217 idx_lookup, % hash table: lookup table to 218 % convert master problem constraint 219 % ids into external solver row ids 220 sp_obj_terms, % array: array of implicit sum terms 221 % forming the obj function of subproblems 222 % in same order as the duals array 223 % that are their coeffs in it 224 mp_cols_added, % int: total number of columns added to MP 225 mp_vars, % [Var1,...,Varj]: list of all mp vars 226 sep_call, % the user supplied separation predicate 227 module, % module: module in which to call the various 228 % user-defined predicates 229 230 % Column Management: 231 232 col_del, % atom: ?,none: 233 % column deletion strategy (currently 234 % unused, but included for later) 235 shelf, % cg_shelf{}: store for data that needs to 236 % be preserved across failure 237 phase, % indicator for current phase in two-phase 238 % method: 239 % 0 : phase 1 240 % -1 : phase 2 241 % this is also used as the "dual" multiplier 242 % for cost coefficient of a column in the 243 % subproblem objective 244 new_columns % record: columns waiting to be 245 % added to the master problem 246 ) 247 ). 248 249:- local struct(cg_shelf( 250 info, % data stored around a bp_node 251 optimal_rc, % best bound on the reduced cost 252 % of a new column (or 'none') 253 sp_sol_cnt % number SP solutions posted via 254 % cg_subproblem_solution/1 255 )). 256 257% ---------------------------------------------------------------------- 258% 259% Basis stabilisation term structure 260% 261% ---------------------------------------------------------------------- 262 263:- export struct( 264 stab_term( 265 idx, % the index of the associated 266 % constraint 267 plus_var, % the Y^{+} stabilisation var, 268 plus_coeff, % the current objective coeff 269 plus_bound, % and upper bound 270 minus_var, % the Y^{-} stabilisation var, 271 minus_coeff, % the current objective coeff 272 minus_bound % and upper bound 273 ) 274 ). 275 276% ---------------------------------------------------------------------- 277% 278% Basis stabilisation counters structure 279% 280% ---------------------------------------------------------------------- 281 282:- export struct( 283 stab_counters( 284 bound_counter, % # iterations 285 % since last bound 286 % updates 287 coeff_counter % # iterations 288 % since coeff bound 289 % updates 290 ) 291 ). 292% ---------------------------------------------------------------------- 293% 294% Subproblem handle structure 295% 296% Note: this info should be moved into the cg_prob structure and accessed 297% from there instead, so that we can quietly ignore the sp_prob struct 298% and eventually drop altogether. 299% 300% 301% ---------------------------------------------------------------------- 302 303:- export struct( 304 sp_prob( 305 master_pool, % atom: the MP pool to which to 306 % post cg_subproblem_solution/1 307 cutoff, % float: the bound for termination 308 % of column generation 309 cost, % var: dual_var for SP solution 310 % cost coefficient in MP 311 coeff_vars, % [Ai,...,An]: list of dual_var 312 % vars for SP solution 313 % coefficients in original 314 % constraints of MP 315 aux, % term: arbitrary additional 316 % data stored in colgen 317 % attribute of MP vars 318 disallow, % list of cstrs to post if 319 % actively preventing 320 % duplicate columns 321 % 2-elem list [Count,Templates] 322 status, % phase1, phase2, degenerate 323 module, % obsolete 324 lo, hi, % implicit variable default bounds 325 type % implicit variable default type 326 ) 327 ). 328 329% ---------------------------------------------------------------------- 330% 331% Subproblem solution structure 332% 333% Note: this should have a boolean "optimal" field defaulting to false 334% if not supplied to allow communication of sp_rc_sum without recourse 335% to separate predicate as now. 336% 337% ---------------------------------------------------------------------- 338 339:- export struct( 340 sp_sol( 341 cost, % number: cost coefficient in MP 342 coeff_vars, % [A1,...,An] list of reals 343 % or sparse list [Idi-Ai,...]: 344 % coefficients in original 345 % constraints of MP 346 aux, % term: arbitrary additional 347 % data 348 lo, % lower bound 349 hi, % upper bound 350 type % type integer or real 351 ) 352 ). 353 354% ---------------------------------------------------------------------- 355% 356% Temporary var info structure 357% 358% ---------------------------------------------------------------------- 359 360:- local struct( 361 cg_var_info( 362 lo, 363 hi, 364 val, 365 reduced_cost, 366 type, 367 attr 368 ) 369 ). 370 371% ---------------------------------------------------------------------- 372% 373% cg attribute handlers 374% 375% ---------------------------------------------------------------------- 376 377unify_colgen(_, Attr) :- 378 var(Attr). % Ignore if not a colgen var 379unify_colgen(Term, Attr) :- 380 compound(Attr), 381 unify_term_colgen(Term, Attr). 382 383:- mode unify_term_colgen(?, +). 384unify_term_colgen(X, Attr) :- 385 nonvar(X), % colgen var and NONVAR - instantiated 386 instantiation_deviates_for_handle(Attr, X). 387unify_term_colgen(Y{colgen:AttrY}, AttrX) :- 388 -?-> 389 unify_colgen_colgen(Y, AttrY, AttrX). 390 391unify_colgen_colgen(_, AttrY, AttrX) :- 392 var(AttrY), % No attribute for this extension 393 AttrX = AttrY. % Transfer the attribute 394unify_colgen_colgen(_, AttrY, AttrX) :- 395 nonvar(AttrY), % colgen var and colgen var 396 unify_handles(AttrX, AttrY). 397 398instantiation_deviates_for_handle(ThisAttr, X) :- 399 ( compound(ThisAttr) -> 400 ThisAttr = colgen{mp_val:Val, 401 cost:Cost, 402 solver:Handle, 403 next:NextAttr}, 404 ( float(X) =:= 0 -> 405 true 406 ; var(Cost) -> 407 printf(warning_output, 408 "Warning: instantiating a variable for" 409 " %w with unspecified cost coefficient to" 410 " %w - assuming zero cost%n", [Handle, X]) 411 ; 412 Handle = cg_prob{const_obj:Const*One}, 413 Const0 is Const + Cost * float(X), 414 setarg(const_obj of cg_prob, Handle, Const0*One) 415 ), 416 ( X =:= Val -> % instantiated to its mp_val 417 true 418 ; % otherwise wake the mp 419 % should probably post a constraint 420 % to the sp disallowing the 421 % corresponding sp solution 422 schedule_suspensions(master_prob of colgen, ThisAttr), 423 wake 424 ), 425 instantiation_deviates_for_handle(NextAttr, X) 426 ; 427 % chain terminated by atom 'end' 428 true 429 ). 430 431unify_handles(ThisAttrX, AttrY) :- 432 ThisAttrX = colgen{solver:Handle, 433 cost:Cost, 434 coeffs:Coeffs, 435 aux:Aux, 436 next:NextAttrX}, 437 remove_cg_attr_for_handle(Handle, AttrY, ThisAttrY, NextAttrY), 438 ( % if Y has an attribute for Handle they must match 439 ThisAttrY = colgen{cost:Cost, 440 coeffs:Coeffs, 441 aux:Aux}, 442 % two variables in the same solver are unified, 443 % send an equality constraint for the two columns 444 % to the external solver and wake it 445 schedule_suspensions(master_prob of colgen, ThisAttrX), 446 wake 447 ; 448 % Y has no attribute for 449 ThisAttrY = end 450 ), 451 % continue with the rest of X and Ys chains 452 unify_handles(NextAttrX, ThisAttrX, NextAttrY). 453 454unify_handles(ThisAttrX, Attr0, AttrY) :- 455 ( compound(ThisAttrX) -> 456 ( compound(AttrY) -> 457 ThisAttrX = colgen{solver:Handle, 458 cost:Cost, 459 coeffs:Coeffs, 460 aux:Aux, 461 next:NextAttrX}, 462 remove_cg_attr_for_handle(Handle, AttrY, ThisAttrY, NextAttrY), 463 ( % if Y has an attribute for Handle they must match 464 ThisAttrY = colgen{cost:Cost, 465 coeffs:Coeffs, 466 aux:Aux}, 467 % two variables in the same solver are unified, 468 % send an equality constraint for the two columns 469 % to the external solver and wake it 470 schedule_suspensions(master_prob of colgen, ThisAttrX), 471 wake 472 ; 473 % Y has no attribute for Handle 474 ThisAttrY = end 475 ), 476 % continue with the rest of X and Ys chains 477 unify_handles(NextAttrX, ThisAttrX, NextAttrY) 478 ; 479 % Ys chain terminated by atom'end' 480 true 481 ) 482 ; 483 % Xs chain terminated by atom 'end' 484 % put the rest of Ys chain here 485 setarg(next of colgen, Attr0, AttrY) 486 ). 487 488% From a cg_attr, removes the attribute corresponding to that for the 489% first argument form the chain and returns it. Fails if none found. 490remove_cg_attr_for_handle(Handle, ThisAttr, Attr, RestAttr) :- 491 % no need to test for var(ThisAttr) in chain 492 ThisAttr = colgen{solver:ThisHandle, next:NextAttr}, 493 (ThisHandle == Handle -> 494 RestAttr = NextAttr, 495 setarg(next of colgen, ThisAttr, end), 496 Attr = ThisAttr 497 ; 498 RestAttr = ThisAttr, 499 dechain_cg_attr_for_handle1(Handle, NextAttr, ThisAttr, Attr) 500 ). 501 502dechain_cg_attr_for_handle1(Handle, ThisAttr, Attr0, Attr) :- 503 % no need to test for var(ThisAttr) in chain 504 ( ThisAttr = colgen{solver:ThisHandle, next:NextAttr} -> 505 (ThisHandle == Handle -> 506 setarg(next of colgen, Attr0, NextAttr), 507 setarg(next of colgen, ThisAttr, end), 508 Attr = ThisAttr 509 ; 510 dechain_cg_attr_for_handle1(Handle, NextAttr, ThisAttr, Attr) 511 ) 512 ; % chain terminated by atom 'end' 513 ThisAttr = Attr 514 ). 515 516% get_bounds handler 517 518cg_get_bounds(_Var{colgen:Attr}, Lwb, Upb) ?- 519 ( var(Attr) -> true 520 ; cg_get_attr_bounds(Attr, -1.0Inf, 1.0Inf, Lwb, Upb) ). 521 522cg_get_attr_bounds(colgen{lo:Lo1, hi:Hi1, next:Next}, 523 Lo0, Hi0, Lo, Hi) ?- 524 Lo1 =< Hi0, 525 Hi1 >= Lo0, 526 Lo2 is max(Lo0, Lo1), 527 Hi2 is min(Hi0, Hi1), 528 cg_get_attr_bounds(Next, Lo2, Hi2, Lo, Hi). 529cg_get_attr_bounds(end, Lo0, Hi0, Lo, Hi) ?- 530 Lo0 = Lo, Hi0 = Hi. 531 532cg_set_bounds(Var{colgen:Attr}, Lwb, Upb) ?- 533 ( var(Attr) -> true 534 ; cg_set_attr_bounds(Var, Attr, Lwb, Upb) ). 535 536pos_inf(1e+20) :- !. 537pos_inf(1.0Inf) :- !. 538 539neg_inf(-1e+20) :- !. 540neg_inf(-1.0Inf) :- !. 541 542cg_set_attr_bounds(Var, Attr, Lwb, Upb) :- 543 ( compound(Attr) -> 544 Attr = colgen{coeffs:Coeffs, 545 solver:Handle, 546 lo:Lwb0, 547 hi:Upb0, 548 type:Type, 549 next:NextAttr}, 550 ( Type == real -> 551 Lwb1 = Lwb, 552 Upb1 = Upb 553 ; Type == integer -> 554 ( neg_inf(Lwb) -> Lwb1 = Lwb ; Lwb1 is fix(ceiling(Lwb)) ), 555 ( pos_inf(Upb) -> Upb1 = Upb ; Upb1 is fix(floor(Upb)) ) 556 ), 557 ( Lwb1 > Lwb0 -> 558 ( (neg_inf(Lwb0) ; neg_inf(Lwb1)) -> 559 Handle = cg_prob{master_prob:MPHandle}, 560 ( var(MPHandle) -> 561 % master problem lp not set up yet, will be 562 % taken care of in cg_solver_setup 563 true 564 ; lp_var_occurrence(Var, MPHandle, _I) -> 565 lp_var_set_bounds(MPHandle, Var, Lwb1, Upb1) 566 ; 567 % Var not yet added to master problem lp, will be 568 % taken care of when added 569 true 570 ), 571 LwbDiff = 0 572 ; 573 LwbDiff is Lwb1 - Lwb0 574 ), 575 setarg(lo of colgen, Attr, Lwb1) 576 ; 577 LwbDiff = 0 578 ), 579 ( Upb1 < Upb0 -> 580 ( (pos_inf(Upb0) ; pos_inf(Upb1)) -> 581 Handle = cg_prob{master_prob:MPHandle}, 582 ( var(MPHandle) -> 583 % master problem lp not set up yet, will be 584 % taken care of in cg_solver_setup 585 true 586 ; lp_var_occurrence(Var, MPHandle, _I) -> 587 lp_var_set_bounds(MPHandle, Var, Lwb1, Upb1) 588 ; 589 % Var not yet added to master problem lp, will be 590 % taken care of when added 591 true 592 ), 593 UpbDiff = 0 594 ; 595 UpbDiff is Upb1 - Upb0 596 ), 597 setarg(hi of colgen, Attr, Upb1) 598 ; 599 UpbDiff = 0 600 ), 601 ( LwbDiff =:= 0, UpbDiff =:= 0 -> 602 true 603 ; 604 Handle = cg_prob{master_prob:MPHandle, 605 sp_solution_call:SolveSubProblem}, 606 ( var(MPHandle) -> 607 % master problem lp not set up yet, will be 608 % taken care of in cg_solver_setup 609 true 610 ; lp_var_occurrence(Var, MPHandle, _I) -> 611 lp_var_set_bounds(MPHandle, Var, Lwb1, Upb1) 612 ; 613 % Var not yet added to master problem lp, will be 614 % taken care of when added 615 true 616 ), 617 ( var(SolveSubProblem) -> 618 true 619 ; var(Coeffs) -> 620 true 621 ; 622 arg(1, SolveSubProblem, sp_prob{coeff_vars:DualVars}), 623 ( 624 foreach(DualVar, DualVars), 625 param(Handle, Coeffs, LwbDiff, UpbDiff) 626 do 627 get_idx(DualVar, Ident, Handle), 628 ( once member(Ident-Val, Coeffs) -> 629 get_lhs_range(DualVar, Lo0, Hi0), 630 Lo1 is Lo0 + Val*LwbDiff, 631 Hi1 is Hi0 + Val*UpbDiff, 632 set_lhs_range(DualVar, Lo1, Hi1) 633 ; 634 true 635 ) 636 ) 637 ) 638 ), 639 cg_set_attr_bounds(Var, NextAttr, Lwb, Upb) 640 ; 641 % chain terminated by atom 'end' 642 true 643 ). 644 645cg_var_print(_{colgen:Attr}, Printed) ?- 646 nonvar(Attr), 647 printed_cg_attributes(Attr, Printed). 648 649printed_cg_attributes(Attr, Printed) :- 650 ( compound(Attr) -> 651 Attr = colgen{solver:Handle, 652 mp_val:Val, 653 cost:Cost, 654 coeffs:Coeffs, 655 aux:Aux, 656 lo:Lo, 657 hi:Hi, 658 next:NextAttr}, 659 Printed = [Handle:[mp_val:Val, cost:Cost, coeffs:Coeffs, 660 aux:Aux, lo:Lo, hi:Hi]|Rest], 661 printed_cg_attributes(NextAttr, Rest) 662 ; 663 % chain terminated by atom 'end' 664 Printed = [] 665 ). 666 667get_cg_attr(X{colgen:Attr0}, Handle, Attr) ?- 668 ( var(Attr0) -> 669 new_cg_attr(X, Handle, Attr) 670 ; 671 Attr0 = colgen{solver:Handle0,next:Next}, 672 % should not fail unless Attr0 incorrect 673 ( Handle0 == Handle -> 674 Attr = Attr0 675 ; 676 get_cg_attr1(Next, Attr0, Handle, Attr) 677 ) 678 ). 679get_cg_attr(X, Handle, Attr) :- % make a new colgen variable 680 free(X), 681 new_cg_attr(X, Handle, Attr). 682 683get_cg_attr1(ThisAttr, Attr0, Handle, Attr) :- 684 atom(ThisAttr), !, % chain terminated by atom 'end' 685 new_cg_attrstruct(Handle, Attr), 686 setarg(next of colgen, Attr0, Attr). 687get_cg_attr1(ThisAttr, _Attr0, Handle, Attr) :- 688 ThisAttr = colgen{solver:Handle0,next:Next}, 689 ( Handle0 == Handle -> 690 Attr = ThisAttr 691 ; 692 get_cg_attr1(Next, ThisAttr, Handle, Attr) 693 ). 694 695new_cg_attr(X, Handle, Attr) :- % make a new colgen variable: 696 new_cg_attrstruct(Handle, Attr), 697 add_attribute(X, Attr, colgen). % and add a colgen attribute 698 699:- mode new_cg_attrstruct(+,?). 700new_cg_attrstruct(Handle, Attr) :- 701 new_cg_attrstruct(Handle, _, _, _, _, _, _, Attr). 702 703 704%:- mode new_cg_attrstruct(+,?,?,?,?,?,?,-). 705new_cg_attrstruct(Handle, Obj, Coeffs, Lo, Hi, Type, Info, Attr) :- 706 ( Lo = -1.0Inf -> true ; true ), 707 ( Hi = 1.0Inf -> true ; true ), 708 ( Type = real -> true ; true ), 709 Attr = colgen{mp_val:0, 710 solver:Handle, 711 cost:Obj, 712 coeffs:Coeffs, 713 lo:Lo, 714 hi:Hi, 715 type:Type, 716 aux:Info, 717 next:end}, 718 init_suspension_list(master_prob of colgen, Attr). 719 720 721% From a cg_attr, searches for the attribute corresponding to that for the 722% first argument. Fails if none found. 723get_cg_attr_for_handle(Handle, Attr0, Attr) :- 724 compound(Attr0), 725 get_cg_attr_for_handle1(Handle, Attr0, Attr). 726 727get_cg_attr_for_handle1(Handle, Attr0, Attr) :- 728 % no need to test for var(Attr0) in chain 729 Attr0 = colgen{solver:Handle0,next:NextAttr}, 730 (Handle0 == Handle -> 731 Attr0 = Attr 732 ; 733 get_cg_attr_for_handle1(Handle, NextAttr, Attr) 734 ). 735 736% var_dual/6: for pools 737var_dual1(Var, Dual, Coeff, Idx, Type, Rhs, Pool) :- 738 get_pool_handle(Handle, Pool), 739 ( Handle == 0 -> 740 printf(error, "Colgen instance has no solver set up in %w%n", 741 [Pool:var_dual(Var, Dual, Coeff, Idx, Type, Rhs)]), 742 abort 743 ; 744 % make sure we consistently have floats in the attribute 745 FDual is float(Dual), 746 FRhs is float(Rhs), 747 var_dual(Var, FDual, Coeff, Idx, Type, FRhs, Handle) 748 ). 749 750% get_dual/2: for pools 751get_dual1(Var, Dual, Pool) :- 752 get_pool_handle(Handle, Pool), 753 ( Handle == 0 -> 754 printf(error, "Colgen instance has no solver set up in %w%n", 755 [Pool:get_dual(Var, Dual)]), 756 abort 757 ; 758 get_dual(Var, Dual, Handle) 759 ). 760 761% get_coeff/2: for pools 762get_coeff1(Var, Coeff, Pool) :- 763 get_pool_handle(Handle, Pool), 764 ( Handle == 0 -> 765 printf(error, "Colgen instance has no solver set up in %w%n", 766 [Pool:get_coeff(Var, Coeff)]), 767 abort 768 ; 769 get_coeff(Var, Coeff, Handle) 770 ). 771 772% get_idx/2: for pools 773get_idx1(Var, Idx, Pool) :- 774 get_pool_handle(Handle, Pool), 775 ( Handle == 0 -> 776 printf(error, "Colgen instance has no solver set up in %w%n", 777 [Pool:get_idx(Var, Idx)]), 778 abort 779 ; 780 get_idx(Var, Idx, Handle) 781 ). 782 783% get_rhs/2: for pools 784get_rhs1(Var, Rhs, Pool) :- 785 get_pool_handle(Handle, Pool), 786 ( Handle == 0 -> 787 printf(error, "Colgen instance has no solver set up in %w%n", 788 [Pool:get_rhs(Var, Rhs)]), 789 abort 790 ; 791 get_rhs(Var, Rhs, Handle) 792 ). 793 794% always_set_dual/2: for pools 795always_set_dual1(Var, Dual, Pool) :- 796 get_pool_handle(Handle, Pool), 797 ( Handle == 0 -> 798 printf(error, "Colgen instance has no solver set up in %w%n", 799 [Pool:always_set_dual(Var, Dual)]), 800 abort 801 ; 802 always_set_dual(Var, Dual, Handle) 803 ). 804 805% set_dual/2: for pools 806set_dual1(Var, Dual, Pool) :- 807 get_pool_handle(Handle, Pool), 808 ( Handle == 0 -> 809 printf(error, "Colgen instance has no solver set up in %w%n", 810 [Pool:set_dual(Var, Dual)]), 811 abort 812 ; 813 set_dual(Var, Dual, Handle) 814 ). 815 816% cg_var_get/4: for low-level handles 817cg_var_get(Handle, X, What, Val) :- 818 ( var(Handle) -> 819 error(4, cg_var_get(Handle, X, What, Val)) 820 ; Handle = cg_prob{} -> 821 cg_var_get_body(Handle, X, What, Val) 822 ; 823 error(5, cg_var_get(Handle, X, What, Val)) 824 ). 825 826% var_get/3: for pools 827cg_var_get1(X, What, Val, Pool) :- 828 get_pool_handle(Handle, Pool), 829 ( Handle == 0 -> 830 printf(error, "Colgen instance has no solver set up in %w%n", 831 [Pool:cg_var_get(X, What, Val)]), 832 abort 833 ; 834 cg_var_get_body(Handle, X, What, Val) 835 ). 836 837cg_var_get_body(Handle, X, node_val, Val) ?- !, 838 Handle = cg_prob{bfs_tree:BfsHandle}, 839 bfs_var_get(BfsHandle, X, node_val, Val). 840cg_var_get_body(Handle, X, reduced_cost, Val) ?- !, 841 Handle = cg_prob{bfs_tree:BfsHandle}, 842 bfs_var_get(BfsHandle, X, node_rc, Val). 843cg_var_get_body(Handle, X, mp_val, Val) ?- !, 844 cg_var_mp_val(Handle, X, Val). 845cg_var_get_body(Handle, X, cost, Val) ?- !, 846 cg_var_cost(Handle, X, Val). 847cg_var_get_body(Handle, X, coeffs, Val) ?- !, 848 cg_var_coeffs(Handle, X, Val). 849cg_var_get_body(Handle, X, aux, Val) ?- !, 850 cg_var_aux(Handle, X, Val). 851cg_var_get_body(Handle, X, What, Val) ?- 852 error(6, cg_var_get_body(Handle, X, What, Val)). 853 854cg_var_mp_val(Handle, _{colgen:Attr0}, Sol) ?- 855 get_cg_attr_for_handle(Handle, Attr0, colgen{mp_val:Sol}). 856cg_var_mp_val(_, X, Sol) :- 857 integer(X), 858 Sol is float(X). 859cg_var_mp_val(_, X, Sol) :- 860 real(X), 861 X = Sol. 862 863cg_var_cost(Handle, _{colgen:Attr0}, Cost) ?- 864 get_cg_attr_for_handle(Handle, Attr0, colgen{cost:Cost}). 865 866cg_var_coeffs(Handle, _{colgen:Attr0}, Coeffs) ?- 867 get_cg_attr_for_handle(Handle, Attr0, colgen{coeffs:Coeffs}). 868 869cg_var_aux(Handle, _{colgen:Attr0}, Aux) ?- 870 get_cg_attr_for_handle(Handle, Attr0, colgen{aux:Aux}). 871 872% cg_integers/2: for low-level handles 873cg_integers(Handle, Ints) :- 874 ( var(Handle) -> 875 error(4, cg_integers(Handle, Ints)) 876 ; Handle = cg_prob{} -> 877 cg_integers_body(Handle, Ints) 878 ; 879 error(5, cg_integers(Handle, Ints)) 880 ). 881 882% integers/1: for pools 883cg_integers1(Ints, Pool) :- 884 get_pool_handle(Handle, Pool), 885 ( Handle == 0 -> 886 printf(error, "Colgen instance has no solver set up in %w%n", 887 [Pool:integers(Ints)]), 888 abort 889 ; 890 cg_integers_body(Handle, Ints) 891 ). 892 893cg_integers_body(Handle, Ints) :- 894 ( Handle = cg_prob{bfs_tree:[],info_messages:OnOff} -> 895 bfs_solver_setup(BfsHandle, min, bp_node(Handle), 896 [info_messages(OnOff), 897 separation(bp_separate(Handle))]), 898 setarg(bfs_tree of cg_prob, Handle, BfsHandle) 899 ; Handle = cg_prob{bfs_tree:BfsHandle} ), 900 ( var(Ints) -> 901 get_cg_attr(Ints, Handle, Attr), 902 setarg(type of colgen, Attr, integer) 903 ; 904 ( 905 foreach(Int, Ints), 906 param(Handle) 907 do 908 get_cg_attr(Int, Handle, Attr), 909 setarg(type of colgen, Attr, integer) 910 ) 911 ), 912 bfs_integers(BfsHandle, Ints). 913 914% ---------------------------------------------------------------------- 915% The user-level constraints 916% ---------------------------------------------------------------------- 917 918% what should this do? look at cg_solver_setup_body 919%cg_add_constraints(Handle, Cstrs, Ids). 920 921cg_range(X, Lo..Hi, Pool) :- 922 get_pool_handle(Handle, Pool), 923 get_cg_attr(X, Handle, _Attr), % make sure it is a var for Pool 924 set_var_bounds(X, Lo, Hi). 925 926cg_eq(X, Y, Pool) :- add_cg_pool_constraint(X$=Y, _Id, Pool). 927cg_ge(X, Y, Pool) :- add_cg_pool_constraint(X$>=Y, _Id, Pool). 928cg_le(X, Y, Pool) :- add_cg_pool_constraint(X$=<Y, _Id, Pool). 929 930add_cg_pool_constraint(Cstr, Ident, Pool) :- 931 cg_normalise_cstr(Cstr, Norm0, CoeffVar, Coeff), 932 !, 933 get_pool_handle(Handle, Pool), 934 Handle = cg_prob{idx_lookup:Lookup,master_prob:MPHandle}, 935 ( nonvar(Ident) -> 936 true 937 ; 938 Id = Ident 939 ), 940 suspend(hash_set(Lookup, Ident, Id), 1, Id->inst), 941 ( nonvar(CoeffVar) -> 942 % only involves existing MP vars 943 try_propagate_bounds(Norm0, Norm), 944 ( var(Norm) -> 945 true % constraint simplified away 946 ; nonvar(MPHandle) -> 947 lp_add_indexed(MPHandle, [Norm], [], [Id]) 948 ; 949 post_typed_pool_constraint(Pool, 950 mp_only of cg_constraint_type, 951 Id:Norm) 952 ) 953 ; 954 % involves MP vars to be generated by SPs 955 % cannot propagate bounds yet 956 % give the CoeffVar of the Lambda vars to be generated 957 % a dual_var attribute 958 Norm0 = Sense:[Val*_One|_], 959 Rhs is float(-Val), 960 var_dual(CoeffVar, 0.0, Coeff, Ident, Sense, Rhs, Handle), 961 post_typed_pool_constraint(Pool, 962 mp_sp of cg_constraint_type, 963 [CoeffVar, Id]:Norm0) 964 ). 965add_cg_pool_constraint(Cstr, _Id, Pool) :- 966 error(5, Pool:Cstr). 967 968try_propagate_bounds(NormIn, NormOut) :- 969 NormIn = Sense:[Cst*_|Lhs], 970 ( Lhs == [] -> % ground: check immediately 971 satisfied(Sense, Cst) 972 ; Lhs = [C*X] -> % single var: update its bound 973 Bound is -Cst/C, 974 swap_sense(C, Sense, Sense1), 975 ( Sense1 == (=<) -> set_var_bounds(X, -1.0Inf, Bound) 976 ; Sense1 == (>=) -> set_var_bounds(X, Bound, 1.0Inf) 977 ; set_var_bounds(X, Bound, Bound) ), % may bind X! 978 wake 979 ; 980 NormIn = NormOut 981 ). 982 983 satisfied((=<), C) :- C =< 0. 984 satisfied((>=), C) :- C >= 0. 985 satisfied((=:=), C) :- C =:= 0. 986 987 swap_sense(C, (=<), (>=)) :- C < 0, !. 988 swap_sense(C, (>=), (=<)) :- C < 0, !. 989 swap_sense(_, S, S). 990 991% cg_sp_count/2: for low-level handles 992%:- deprecated(cg_sp_count/2, "No longer necessary"). 993cg_sp_count(Handle, P) :- 994 ( var(Handle) -> 995 error(4, cg_sp_count(Handle, P)) 996 ; var(P) -> 997 error(4, cg_sp_count(Handle, P)) 998 ; Handle = cg_prob{}, integer(P) -> 999 cg_sp_count_body(Handle, P) 1000 ; 1001 error(5, cg_sp_count(Handle, P)) 1002 ). 1003 1004% cg_subproblem_count/1: for pools 1005cg_sp_count1(P, Pool) :- 1006 get_pool_handle(Handle, Pool), 1007 ( Handle == 0 -> 1008 printf(error, "Colgen instance has no solver set up in %w%n", 1009 [Pool:cg_subproblem_count(P)]), 1010 abort 1011 ; 1012 cg_sp_count_body(Handle, P) 1013 ). 1014 1015cg_sp_count_body(_Handle, _P). 1016 1017:- deprecated(cg_sp_sol/2, "Use valid_columns/1"). 1018% cg_subproblem_solution/1: for pools 1019cg_sp_sol(ColSpecs, Pool) :- 1020 cg_valid_columns1(ColSpecs, Pool). 1021 1022% valid_columns/1: for pools 1023cg_valid_columns1(ColSpecs, Pool) :- 1024 get_pool_handle(Handle, Pool), 1025 ( Handle == 0 -> 1026 printf(error, "Colgen instance has no solver set up in %w%n", 1027 [Pool:valid_columns(ColSpecs)]), 1028 abort 1029 ; 1030 Handle = cg_prob{shelf:Shelf}, 1031 shelf_inc(Shelf, sp_sol_cnt of cg_shelf), % suppress auto-posting 1032 cg_valid_columns(Handle, ColSpecs) 1033 ). 1034 1035% cg_valid_columns/2: for low-level handles 1036cg_valid_columns(Handle, ColSpecs) :- 1037 ( var(Handle) -> 1038 error(4, cg_valid_columns(Handle, ColSpecs)) 1039 ; var(ColSpecs) -> 1040 error(4, cg_valid_columns(Handle, ColSpecs)) 1041 ; Handle = cg_prob{new_columns:NewColRec} -> 1042 % recorda for compatibility, should be recordz 1043 ( ColSpecs = sp_sol{} -> 1044 recorda(NewColRec, ColSpecs) 1045 ; 1046 ( foreach(ColSpec, ColSpecs), param(NewColRec) do 1047 ColSpec = sp_sol{}, 1048 recorda(NewColRec, ColSpecs) 1049 ) 1050 ) 1051 ; 1052 error(5, cg_valid_columns(Handle, ColSpecs)) 1053 ). 1054 1055 1056% Set optimal RC (during SP solution process): the SP solver should 1057% do that iff it computes an RC-optimal column, because we can derive 1058% an MP lower bound from it. 1059 1060% cg_optimal_rc/2: for low-level handles 1061cg_optimal_rc(Handle, RCOpt) :- 1062 ( var(Handle) -> 1063 error(4, cg_optimal_rc(Handle, RCOpt)) 1064 ; var(RCOpt) -> 1065 error(4, cg_optimal_rc(Handle, RCOpt)) 1066 ; Handle = cg_prob{shelf:Shelf}, number(RCOpt) -> 1067 shelf_set(Shelf, optimal_rc of cg_shelf, RCOpt) 1068 ; 1069 error(5, cg_optimal_rc(Handle, RCOpt)) 1070 ). 1071 1072% optimal_rc/1: for pools 1073cg_optimal_rc1(RCOpt, Pool) :- 1074 get_pool_handle(Handle, Pool), 1075 ( Handle == 0 -> 1076 printf(error, "Colgen instance has no solver set up in %w%n", 1077 [Pool:optimal_rc(RCOpt)]), 1078 abort 1079 ; 1080 cg_optimal_rc(Handle, RCOpt) 1081 ). 1082 1083:- deprecated(cg_sp_rc_sum/2, "Use optimal_rc/1"). 1084% subproblem_rc_sum/1: for pools 1085cg_sp_rc_sum(SPRCSum, Pool) :- 1086 cg_optimal_rc1(SPRCSum, Pool). 1087 1088 1089% cg_branch/2: for low-level handles 1090cg_branch(Handle, Branch) :- 1091 ( var(Handle) -> 1092 error(4, cg_branch(Handle, Branch)) 1093 ; Handle = cg_prob{} -> 1094 cg_branch_body(Handle, 0, Branch) 1095 ; 1096 error(5, cg_branch(Handle, Branch)) 1097 ). 1098 1099% cg_branch/3: for low-level handles 1100cg_branch(Handle, Score, Branch) :- 1101 ( var(Handle) -> 1102 error(4, cg_branch(Handle, Score, Branch)) 1103 ; Handle = cg_prob{} -> 1104 cg_branch_body(Handle, Score, Branch) 1105 ; 1106 error(5, cg_branch(Handle, Score, Branch)) 1107 ). 1108 1109% branch/1: for pools 1110cg_branch1(Branch, Pool) :- 1111 get_pool_handle(Handle, Pool), 1112 ( Handle == 0 -> 1113 printf(error, "Colgen instance has no solver set up in %w%n", 1114 [Pool:branch(Branch)]), 1115 abort 1116 ; 1117 cg_branch_body(Handle, 0, Branch) 1118 ). 1119 1120% branch/2: for pools 1121cg_branch1(Score, Branch, Pool) :- 1122 get_pool_handle(Handle, Pool), 1123 ( Handle == 0 -> 1124 printf(error, "Colgen instance has no solver set up in %w%n", 1125 [Pool:branch(Score, Branch)]), 1126 abort 1127 ; 1128 cg_branch_body(Handle, Score, Branch) 1129 ). 1130 1131cg_branch_body(cg_prob{shelf:Shelf}, Score, Branch) :- 1132 shelf_get(Shelf, info of cg_shelf, Branches), 1133 shelf_set(Shelf, info of cg_shelf, [Score:Branch|Branches]). 1134 1135cg_info_message(cg_prob{info_messages:OnOff}, String, Vars) :- 1136 ( OnOff == on -> printf(String, Vars), flush(output) ; true ). 1137 1138% cg_get/3: for low-level handles 1139cg_get(Handle, What, Val) :- 1140 ( var(Handle) -> 1141 error(4, cg_get(Handle, What, Val)) 1142 ; Handle = cg_prob{} -> 1143 cg_get_body(Handle, What, Val) 1144 ; 1145 error(5, cg_get(Handle, What, Val)) 1146 ). 1147 1148% get/2: for pools 1149cg_get1(What, Val, Pool) :- 1150 get_pool_handle(Handle, Pool), 1151 ( Handle == 0 -> 1152 printf(error, "Colgen instance has no solver set up in %w%n", 1153 [Pool:cg_get(What, Val)]), 1154 abort 1155 ; 1156 cg_get_body(Handle, What, Val) 1157 ). 1158 1159cg_get_body(Handle, obj_val, ObjVal) ?- !, 1160 Handle = cg_prob{master_prob:MPHandle, 1161 mp_vars:Vars, 1162 const_obj:Const*_One}, 1163 ( 1164 foreach(Var, Vars), 1165 fromto(Const, In, Out, ObjVal), 1166 param(Handle, MPHandle) 1167 do 1168 ( nonvar(Var) -> 1169 Out = In 1170 ; 1171 cg_var_cost(Handle, Var, Cost), 1172 lp_var_get(MPHandle, Var, solution, VarVal), 1173 Out is In + Cost*VarVal 1174 ) 1175 ). 1176cg_get_body(Handle, sp_obj, Val) :- !, 1177 cg_get_body(Handle, sp_obj(all), Val). 1178cg_get_body(Handle, sp_obj(Idents), Val) :- !, 1179 Handle = cg_prob{duals:DualArr, 1180 idx_lookup:Lookup, 1181 sp_obj_terms:TermArr}, 1182 ( Idents == all -> 1183 ( 1184 for(I,1,arity(TermArr)), 1185 fromto(Val, Out, In, []), 1186 param(TermArr,DualArr) 1187 do 1188 arg(I, TermArr, Term), 1189 arg(I, DualArr, Dual), 1190 ( nonvar(Term), Term =:= 0 -> Out = In 1191 ; Out = [Dual*Term|In] ) 1192 ) 1193 ; Idents = [_|_] -> 1194 ( 1195 foreach(Ident, Idents), 1196 foreach(Dual*Term, Val), 1197 param(Lookup, DualArr, TermArr) 1198 do 1199 hash_get(Lookup, Ident, Id), 1200 Id1 is Id + 1, 1201 arg(Id1, DualArr, Dual), 1202 arg(Id1, TermArr, Term) 1203 ) 1204 ; 1205 hash_get(Lookup, Idents, Id), 1206 Id1 is Id + 1, 1207 arg(Id1, DualArr, Dual), 1208 arg(Id1, TermArr, Term), 1209 Val = Dual*Term 1210 ). 1211cg_get_body(Handle, dual(Idents), Val) :- !, 1212 Handle = cg_prob{duals:DualArr, 1213 idx_lookup:Lookup}, 1214 ( Idents == all -> 1215 DualArr =.. [[]|Val] 1216 ; Idents = [_|_] -> 1217 ( 1218 foreach(Ident, Idents), 1219 foreach(V, Val), 1220 param(Lookup, DualArr) 1221 do 1222 hash_get(Lookup, Ident, Id), 1223 Id1 is Id + 1, 1224 arg(Id1, DualArr, V) 1225 ) 1226 ; 1227 hash_get(Lookup, Idents, Id), 1228 Id1 is Id + 1, 1229 arg(Id1, DualArr, Val) 1230 ). 1231cg_get_body(Handle, column_count, Val) ?- !, 1232 Handle = cg_prob{mp_cols_added:Val}. 1233cg_get_body(Handle, unsatisfiable_cstrs, Val) ?- !, 1234 Handle = cg_prob{sp_solution_call:SolveSubProblem}, 1235 arg(1, SolveSubProblem, sp_prob{coeff_vars:DualVars}), 1236 ( 1237 foreach(DualVar, DualVars), 1238 fromto(Val, Out, In, []), 1239 param(Handle) 1240 do 1241 ( satisfiable_primal_cstr(DualVar, Handle) -> Out = In ; Out = [DualVar|In] ) 1242 ). 1243cg_get_body(Handle, satisfiable_cstrs, Val) ?- !, 1244 Handle = cg_prob{sp_solution_call:SolveSubProblem}, 1245 arg(1, SolveSubProblem, sp_prob{coeff_vars:DualVars}), 1246 ( 1247 foreach(DualVar, DualVars), 1248 fromto(Val, Out, In, []), 1249 param(Handle) 1250 do 1251 ( satisfiable_primal_cstr(DualVar, Handle) -> Out = [DualVar|In] ; Out = In ) 1252 ). 1253cg_get_body(Handle, frac_vars, Val) ?- !, 1254 Handle = cg_prob{bfs_tree:BfsHandle}, 1255 bfs_get(BfsHandle, frac_vars, Val). 1256cg_get_body(Handle, generated_non_zero_vars, Val) ?- !, 1257 Handle = cg_prob{mp_vars:Vars}, 1258 ( 1259 foreach(Var, Vars), 1260 fromto(Val, Out, In, []), 1261 param(Handle) 1262 do 1263 ( nonvar(Var) -> 1264 ( Var > 1e-05 -> Out = [Var|In] ; Out = In ) 1265 ; 1266 get_cg_attr(Var, Handle, colgen{mp_val:Sol, coeffs:Coeffs}), 1267 ( abs(Sol) =< 1e-05 -> Out = In 1268 ; Coeffs == [] -> Out = In 1269 ; Out = [Var|In] ) 1270 ) 1271 ). 1272cg_get_body(Handle, non_zero_vars, Val) ?- !, 1273 Handle = cg_prob{mp_vars:Vars}, 1274 ( 1275 foreach(Var, Vars), 1276 fromto(Val, Out, In, []), 1277 param(Handle) 1278 do 1279 ( nonvar(Var) -> 1280 ( Var > 1e-05 -> Out = [Var|In] ; Out = In ) 1281 ; 1282 get_cg_attr(Var, Handle, colgen{mp_val:Sol}), 1283 ( abs(Sol) =< 1e-05 -> Out = In ; Out = [Var|In] ) 1284 ) 1285 ). 1286cg_get_body(Handle, generated_vars, Val) :- !, 1287 Handle = cg_prob{mp_vars:Vars}, 1288 ( 1289 foreach(Var, Vars), 1290 fromto(Val, Out, In, []), 1291 param(Handle) 1292 do 1293 ( nonvar(Var) -> 1294 Out = In 1295 ; 1296 get_cg_attr(Var, Handle, colgen{coeffs:Coeffs}), 1297 ( Coeffs == [] -> Out = In 1298 ; Out = [Var|In] ) 1299 ) 1300 ). 1301cg_get_body(Handle, vars, Val) ?- !, 1302 Handle = cg_prob{mp_vars:Val}. 1303cg_get_body(Handle, sep_goal, Val) ?- !, 1304 Handle = cg_prob{sep_call:(call(Val)@_Module)}. 1305cg_get_body(Handle, sp_solver, Val) ?- !, 1306 Handle = cg_prob{sp_solution_call:Val}. 1307cg_get_body(Handle, stab_coeff_minus(Ident), Val) :- !, 1308 Handle = cg_prob{stab_terms:StabTerms, 1309 idx_lookup:Lookup}, 1310 hash_get(Lookup, Ident, Id), 1311 get_stab_term(StabTerms, Id, stab_term{minus_coeff:Val}). 1312cg_get_body(Handle, stab_coeff_plus(Ident), Val) :- !, 1313 Handle = cg_prob{stab_terms:StabTerms, 1314 idx_lookup:Lookup}, 1315 hash_get(Lookup, Ident, Id), 1316 get_stab_term(StabTerms, Id, stab_term{plus_coeff:Val}). 1317cg_get_body(Handle, stab_bound_minus(Ident), Val) :- !, 1318 Handle = cg_prob{stab_terms:StabTerms, 1319 idx_lookup:Lookup}, 1320 hash_get(Lookup, Ident, Id), 1321 get_stab_term(StabTerms, Id, stab_term{minus_bound:Val}). 1322cg_get_body(Handle, stab_bound_plus(Ident), Val) :- !, 1323 Handle = cg_prob{stab_terms:StabTerms, 1324 idx_lookup:Lookup}, 1325 hash_get(Lookup, Ident, Id), 1326 get_stab_term(StabTerms, Id, stab_term{plus_bound:Val}). 1327cg_get_body(Handle, What, Val) :- 1328 error(6, cg_get_body(Handle, What, Val)). 1329 1330% cg_set/3: for low-level handles 1331cg_set(Handle, What, Val) :- 1332 ( var(Handle) -> 1333 error(4, cg_set(Handle, What, Val)) 1334 ; Handle = cg_prob{} -> 1335 cg_set_body(Handle, What, Val) 1336 ; 1337 error(5, cg_set(Handle, What, Val)) 1338 ). 1339 1340% set/2: for pools 1341cg_set1(What, Val, Pool) :- 1342 get_pool_handle(Handle, Pool), 1343 ( Handle == 0 -> 1344 printf(error, "Colgen instance has no solver set up in %w%n", 1345 [Pool:cg_set(What, Val)]), 1346 abort 1347 ; 1348 cg_set_body(Handle, What, Val) 1349 ). 1350 1351cg_set_body(Handle, disallow, Val) ?- !, 1352 setarg(disallow of cg_prob, Handle, Val). 1353cg_set_body(Handle, int_tolerance, Val) ?- !, 1354 setarg(tolerance of cg_prob, Handle, Val). 1355cg_set_body(Handle, basis_perturbation, Val) ?- !, 1356 setarg(basis_perturbation of cg_prob, Handle, Val). 1357cg_set_body(Handle, info_messages, Val) ?- !, 1358 setarg(info_messages of cg_prob, Handle, Val), 1359 Handle = cg_prob{bfs_tree:BfsHandle}, 1360 ( BfsHandle == [] -> true ; bfs_set(BfsHandle, info_messages, Val) ). 1361cg_set_body(Handle, on_degeneracy, Val) ?- !, 1362 setarg(on_degeneracy of cg_prob, Handle, Val). 1363cg_set_body(Handle, stabilisation, Val) ?- !, 1364 setarg(stabilisation of cg_prob, Handle, Val). 1365cg_set_body(Handle, stab_coeff_minus(Ident), Val) :- !, 1366 Handle = cg_prob{stab_terms:StabTerms, 1367 idx_lookup:Lookup}, 1368 hash_get(Lookup, Ident, Id), 1369 get_stab_term(StabTerms, Id, StabTerm), 1370 setarg(minus_coeff of stab_term, StabTerm, Val). 1371cg_set_body(Handle, stab_coeff_plus(Ident), Val) :- !, 1372 Handle = cg_prob{stab_terms:StabTerms, 1373 idx_lookup:Lookup}, 1374 hash_get(Lookup, Ident, Id), 1375 get_stab_term(StabTerms, Id, StabTerm), 1376 setarg(plus_coeff of stab_term, StabTerm, Val). 1377cg_set_body(Handle, stab_bound_minus(Ident), Val) :- !, 1378 Handle = cg_prob{master_prob:MPHandle, 1379 stab_terms:StabTerms, 1380 idx_lookup:Lookup}, 1381 hash_get(Lookup, Ident, Id), 1382 get_stab_term(StabTerms, Id, StabTerm), 1383 StabTerm = stab_term{minus_var:Yminus}, 1384 lp_var_non_monotonic_set_bounds(MPHandle, Yminus, 0, Val), 1385 setarg(minus_bound of stab_term, StabTerm, Val). 1386cg_set_body(Handle, stab_bound_plus(Ident), Val) :- !, 1387 Handle = cg_prob{master_prob:MPHandle, 1388 stab_terms:StabTerms, 1389 idx_lookup:Lookup}, 1390 hash_get(Lookup, Ident, Id), 1391 get_stab_term(StabTerms, Id, StabTerm), 1392 StabTerm = stab_term{plus_var:Yplus}, 1393 lp_var_non_monotonic_set_bounds(MPHandle, Yplus, 0, Val), 1394 setarg(plus_bound of stab_term, StabTerm, Val). 1395cg_set_body(Handle, What, Val) :- 1396 error(6, cg_set_body(Handle, What, Val)). 1397 1398% cg_statistics/1: for low-level handles 1399cg_statistics(Handle) :- 1400 ( var(Handle) -> 1401 error(4, cg_statistics(Handle)) 1402 ; Handle = cg_prob{} -> 1403 cg_statistics_body(Handle) 1404 ; 1405 error(5, cg_statistics(Handle)) 1406 ). 1407 1408% statistics/0: for pools 1409cg_statistics1(Pool) :- 1410 get_pool_handle(Handle, Pool), 1411 ( Handle == 0 -> 1412 printf(error, "Colgen instance has no solver set up in %w%n", 1413 [Pool:statistics]), 1414 abort 1415 ; 1416 cg_statistics_body(Handle) 1417 ). 1418 1419cg_statistics_body(Handle) :- 1420 Handle = cg_prob{bfs_tree:BfsHandle}, 1421 ( BfsHandle == [] -> true ; bfs_statistics(BfsHandle) ). 1422 1423get_stab_term([StabTerm|StabTerms], Id, Val) :- 1424 StabTerm = stab_term{idx:Idx}, 1425 ( Idx == Id -> 1426 Val = StabTerm 1427 ; 1428 get_stab_term(StabTerms, Id, Val) 1429 ). 1430 1431% ---------------------------------------------------------------------- 1432% The optimisation predicates 1433% ---------------------------------------------------------------------- 1434 1435% cg_minimize/8,9: for low-level handles 1436:- tool(cg_minimize/8, cg_minimize0/9). 1437 1438cg_minimize0(MPCstrs, MPIdxs, MPSPCstrs, MPSPIdxs, CoeffVars, SolveSubProblem, Obj, ObjVal, Module) :- 1439 cg_minimize(MPCstrs, MPIdxs, 1440 MPSPCstrs, MPSPIdxs, CoeffVars, 1441 SolveSubProblem, Obj, [], 1442 ObjVal, Module). 1443 1444:- tool(cg_minimize/9, cg_minimize/10). 1445 1446cg_minimize(MPCstrs, MPIdxs, MPSPCstrs, MPSPIdxs, CoeffVars, 1447 SolveSubProblem, Obj, OptionList, ObjVal, Module) :- 1448 cg_minimize_body(_Handle, MPCstrs, MPIdxs, 1449 MPSPCstrs, MPSPIdxs, CoeffVars, 1450 SolveSubProblem, Obj, OptionList, ObjVal, 1451 _, Module). 1452 1453% minimize/4,5: for pools 1454:- tool(cg_minimize/4, cg_minimize0/5). 1455 1456cg_minimize0(SolveSubProblem, Obj, ObjVal, Pool, Module) :- 1457 cg_minimize(SolveSubProblem, Obj, [], ObjVal, Pool, Module). 1458 1459:- tool(cg_minimize/5, cg_minimize/6). 1460 1461cg_minimize(SolveSubProblem, Obj, OptionList, ObjVal, Pool, Module) :- 1462 get_pool_item(Pool, Handle), 1463 ( Handle == 0 -> 1464 printf(error, "Colgen instance has no solver in %w%n", 1465 [Pool:minimize(SolveSubProblem, Obj, OptionList, ObjVal)]), 1466 abort 1467 ; 1468 cg_pool_collect_mp_sp_cstrs(Pool, MPCstrs, MPIdxs, 1469 MPSPCstrs, MPSPIdxs, CoeffVars), 1470 cg_minimize_body(Handle, MPCstrs, MPIdxs, 1471 MPSPCstrs, MPSPIdxs, CoeffVars, 1472 SolveSubProblem, Obj, OptionList, ObjVal, 1473 Pool, Module) 1474 ). 1475 1476cg_pool_collect_mp_sp_cstrs(Pool, MPCstrs, MPIdxs, 1477 MPSPCstrs, MPSPIdxs, CoeffVars) :- 1478 % collect any constraints only involving known MP vars 1479 collect_typed_pool_constraints(Pool, mp_only of cg_constraint_type, 1480 MPIdxCstrs), 1481 ( 1482 foreach(Id:MPCstr, MPIdxCstrs), 1483 foreach(Id, MPIdxs), 1484 foreach(MPCstr, MPCstrs) 1485 do 1486 true 1487 ), 1488 % collect any constraints constraints which involve generated vars 1489 collect_typed_pool_constraints(Pool, mp_sp of cg_constraint_type, 1490 MPSPIdxCstrs), 1491 ( 1492 foreach([CoeffVar, Id]:MPSPCstr, MPSPIdxCstrs), 1493 foreach(CoeffVar, CoeffVars), 1494 foreach(Id, MPSPIdxs), 1495 foreach(MPSPCstr, MPSPCstrs) 1496 do 1497 true 1498 ). 1499 1500cg_minimize_body(Handle, MPCstrs, MPIdxs, MPSPCstrs, MPSPIdxs, CoeffVars, 1501 SolveSubProblem, Obj, OptionList, ObjVal, Pool, Module) :- 1502 % setup original MP problem 1503 cg_solver_setup_body(Handle, MPCstrs, MPIdxs, 1504 MPSPCstrs, MPSPIdxs, CoeffVars, 1505 SolveSubProblem, Obj, OptionList, 1506 Pool, Module), 1507 % solve the initial MP 1508 cg_iteration(Handle), 1509 Handle = cg_prob{master_prob:MP, 1510 upper_bound:ObjVal}, 1511 lp_get(MP, vars, AllVarArr), 1512 ( 1513 foreacharg(Var, AllVarArr), 1514 param(MP, Handle) 1515 do 1516 ( nonvar(Var) -> 1517 true 1518 ; 1519 lp_var_get(MP, Var, solution, Val), 1520 get_cg_attr(Var, Handle, Attr), 1521 setarg(mp_val of colgen, Attr, Val) 1522 ) 1523 ). 1524 1525% bp_solve/2: for low-level handles 1526bp_solve(Handle, Obj) :- 1527 ( var(Handle) -> 1528 error(4, bp_solve(Handle, Obj)) 1529 ; Handle = cg_prob{} -> 1530 bp_solve_body(Handle, Obj) 1531 ; 1532 error(5, bp_solve(Handle, Obj)) 1533 ). 1534 1535% solve/1: for pools 1536bp_solve1(Obj, Pool) :- 1537 get_pool_item(Pool, Handle), 1538 ( Handle == 0 -> 1539 printf(error, "Colgen instance has no solver set up in %w%n", 1540 [Pool:solve(Obj)]), 1541 abort 1542 ; 1543 bp_solve_body(Handle, Obj) 1544 ). 1545 1546bp_solve_body(Handle, Obj) :- 1547 Handle = cg_prob{bfs_tree:BfsHandle0,info_messages:OnOff}, 1548 ( BfsHandle0 == [] -> 1549 bfs_solver_setup(BfsHandle, min, bp_node(Handle), 1550 [info_messages(OnOff), 1551 separation(bp_separate(Handle))]), 1552 setarg(bfs_tree of cg_prob, Handle, BfsHandle) 1553 ; 1554 BfsHandle = BfsHandle0 1555 ), 1556 bfs_solve(BfsHandle, Obj), 1557 Handle = cg_prob{mp_vars:Vars}, 1558 ( 1559 foreach(Var, Vars), 1560 param(Handle, BfsHandle) 1561 do 1562 ( var(Var) -> 1563 bfs_var_get(BfsHandle, Var, optimal_val, Val), 1564 ( var(Val) -> Val = 0 ; true ), 1565 get_cg_attr(Var, Handle, Attr), 1566 setarg(mp_val of colgen, Attr, Val) 1567 ; 1568 true 1569 ) 1570 ). 1571 1572bp_node(Handle) :- 1573 add_new_solver_rowcols(Handle), 1574 Handle = cg_prob{master_prob:MP, 1575 bfs_tree:BfsHandle, 1576 idx_lookup:Lookup, 1577 shelf:Shelf}, 1578 \+ \+ ( bfs_impose_node_state(other, BfsHandle), 1579 ( cg_get(Handle, unsatisfiable_cstrs, []) -> D = -1 ; D = 0 ), 1580 setarg(phase of cg_prob, Handle, D), 1581 cg_get(Handle, vars, OldVars), 1582 cg_iteration(Handle), 1583 cg_get(Handle, non_zero_vars, Vars), 1584 ( 1585 foreach(Var, OldVars), 1586 foreach(Info, OldVals), 1587 param(Handle, MP) 1588 do 1589 cg_var_get(Handle, Var, mp_val, Val), 1590 % note that if the optimisation has been stopped 1591 % because we have "close enough" bounds there may 1592 % be columns that have been added to the MP since 1593 % the last lp_solve and any attempt to get a 1594 % reduced cost will fail 1595 ( lp_var_get(MP, Var, reduced_cost, RC) -> 1596 true 1597 ; 1598 RC = 0 1599 ), 1600 get_var_bounds(Var, Lo, Hi), 1601 Info = cg_var_info{lo: Lo, 1602 hi: Hi, 1603 val: Val, 1604 reduced_cost: RC} 1605 ), 1606 ( 1607 foreach(Var, Vars), 1608 fromto(NewVals, Out, In, []), 1609 param(OldVars, Handle, MP) 1610 do 1611 % if it was added to the cg_prob before 1612 % this call to cg_iteration/1 its 1613 % solution info has been found above, 1614 % otherwise it must have been generated 1615 % in this call, so we have to save its 1616 % colgen attr as well as the solution info 1617 % to recreate the var after the \+ \+ 1618 ( var_member(Var, OldVars) -> 1619 Out = In 1620 ; 1621 cg_var_get(Handle, Var, mp_val, Val), 1622 % note that if the optimisation has been stopped 1623 % because we have "close enough" bounds there may 1624 % be columns that have been added to the MP since 1625 % the last lp_solve and any attempt to get a 1626 % reduced cost will fail 1627 ( lp_var_get(MP, Var, reduced_cost, RC) -> 1628 true 1629 ; 1630 RC = 0 1631 ), 1632 get_var_bounds(Var, Lo, Hi), 1633 get_cg_attr(Var, Handle, Attr), 1634 Attr = colgen{mp_val:Val, 1635 type:Type, 1636 cost:Cost, 1637 coeffs:Coeffs, 1638 aux:Aux}, 1639 Attr1 = colgen{cost:Cost, 1640 coeffs:Coeffs, 1641 aux:Aux}, 1642 Info = cg_var_info{lo: Lo, 1643 hi: Hi, 1644 val: Val, 1645 reduced_cost: RC, 1646 type: Type, 1647 attr: Attr1}, 1648 Out = [Info|In] 1649 ) 1650 ), 1651 Handle = cg_prob{upper_bound:NodeCost}, 1652 shelf_set(Shelf, info of cg_shelf, [NodeCost, OldVals, NewVals]) 1653 ), 1654 shelf_get(Shelf, info of cg_shelf, [NodeCost, OldVals, NewVals]), 1655 shelf_set(Shelf, info of cg_shelf, []), 1656 cg_get(Handle, vars, OldVars), 1657 ( 1658 foreach(Var, OldVars), 1659 foreach(Info, OldVals), 1660 param(BfsHandle) 1661 do 1662 ( var(Var) -> 1663 Info = cg_var_info{lo: Lo, 1664 hi: Hi, 1665 val: Val, 1666 reduced_cost: RC}, 1667 bfs_node_info(BfsHandle, Var, Lo, Hi, Val, RC) 1668 ; 1669 true 1670 ) 1671 ), 1672 ( 1673 foreach(Info, NewVals), 1674 foreach(Var:ObjCol, VarCols), 1675 fromto(NewVars, [Var|Vars], Vars, OldVars), 1676 param(BfsHandle, Lookup, Handle) 1677 do 1678 Info = cg_var_info{lo: Lo, 1679 hi: Hi, 1680 val: Val, 1681 reduced_cost: RC, 1682 type: Type, 1683 attr: Attr}, 1684 get_cg_attr(Var, Handle, Attr), 1685 Attr = colgen{cost:Obj, coeffs:Coeffs}, 1686 ( Obj =:= 0 -> ObjCol = BCol ; ObjCol = [obj:Obj|BCol] ), 1687 BCol = [lo:Lo, hi:Hi|Col], 1688 ( 1689 foreach(Ident-V, Coeffs), 1690 foreach(Id:V, Col), 1691 param(Lookup) 1692 do 1693 hash_get(Lookup, Ident, Id) 1694 ), 1695 ( Type == integer -> 1696 bfs_integers(BfsHandle, Var), cg_integers(Handle, Var) 1697 ; 1698 true 1699 ), 1700 bfs_node_info(BfsHandle, Var, Lo, Hi, Val, RC) 1701 ), 1702 setarg(mp_vars of cg_prob, Handle, NewVars), 1703 lp_add_columns(MP, VarCols), 1704 lp_get(MP, vars, VarArr), 1705 arity(VarArr, ColsAdded), 1706 setarg(mp_cols_added of cg_prob, Handle, ColsAdded), 1707 bfs_node_cost(BfsHandle, NodeCost). 1708 1709 1710var_member(Var, [H|T]) :- 1711 ( Var == H -> true ; var_member(Var, T) ). 1712 1713% cg_solver_setup/8,9: for low-level handles 1714:- tool(cg_solver_setup/8, cg_solver_setup0/9). 1715 1716cg_solver_setup0(Handle, MPCstrs, MPIdxs, MPSPCstrs, MPSPIdxs, CoeffVars, 1717 SubProblemSolver, Obj, Module) :- 1718 cg_solver_setup(Handle, MPCstrs, MPIdxs, MPSPCstrs, MPSPIdxs, CoeffVars, 1719 SubProblemSolver, Obj, [], Module). 1720 1721:- tool(cg_solver_setup/9, cg_solver_setup/10). 1722 1723cg_solver_setup(Handle, MPCstrs, MPIdxs, MPSPCstrs, MPSPIdxs, CoeffVars, 1724 SubProblemSolver, Obj, OptionList, Module) :- 1725 ( var(Handle) -> 1726 error(4, cg_solver_setup(Handle, MPCstrs, MPIdxs, MPSPCstrs, MPSPIdxs, CoeffVars, 1727 SubProblemSolver, Obj, OptionList)) 1728 ; Handle = cg_prob{} -> 1729 cg_solver_setup_body(Handle, MPCstrs, MPIdxs, 1730 MPSPCstrs, MPSPIdxs, CoeffVars, 1731 SubProblemSolver, Obj, OptionList, 1732 _, Module) 1733 ; 1734 error(5, cg_solver_setup(Handle, MPCstrs, MPIdxs, MPSPCstrs, MPSPIdxs, CoeffVars, 1735 SubProblemSolver, Obj, OptionList)) 1736 ). 1737 1738% solver_setup/2,3: for pools 1739:- tool(cg_solver_setup/3, cg_solver_setup0/4). 1740 1741cg_solver_setup0(SubProblemSolver, Obj, Pool, Module) :- 1742 cg_solver_setup(SubProblemSolver, Obj, [], Pool, Module). 1743 1744:- tool(cg_solver_setup/4, cg_solver_setup/5). 1745 1746cg_solver_setup(SubProblemSolver, Obj, OptionList, Pool, Module) :- 1747 get_pool_item(Pool, Handle), 1748 ( Handle == 0 -> 1749 printf(error, "Colgen instance has no solver in %w%n", 1750 [Pool:solver_setup(SubProblemSolver, Obj, OptionList)]), 1751 abort 1752 ; 1753 cg_pool_collect_mp_sp_cstrs(Pool, MPCstrs, MPIdxs, 1754 MPSPCstrs, MPSPIdxs, CoeffVars), 1755 cg_solver_setup_body(Handle, MPCstrs, MPIdxs, 1756 MPSPCstrs, MPSPIdxs, CoeffVars, 1757 SubProblemSolver, Obj, OptionList, 1758 Pool, Module) 1759 ). 1760 1761%:- mode cg_solver_setup_body(-, +, ?, +, ?, ?, +, ?, +, ?, ++). 1762cg_solver_setup_body(Handle, MPCstrs, MPIdxs, MPSPCstrs, MPSPIdxs, CoeffVars, 1763 SolveSubProblem, Obj, OptionList, Pool, Module) :- 1764 linearize(Obj, [ConstTerm|LinObj], NonLinObj), 1765 ( NonLinObj == [] -> 1766 % the variables to be generated do not appear in the 1767 % objective, set OptVar in the SPs to 0 1768 NormObj = LinObj, 1769 OptVar = 0 1770 ; 1771 % the variables to be generated do appear in the 1772 % objective, but give them a 0 dual val for now 1773 % in case we need to use two phase 1774 NonLinObj = [AuxVar = implicit_sum(OptVar)], 1775 filter_auxvar(AuxVar, LinObj, NormObj, ObjCoeff), 1776 var_dual(OptVar, 0.0, ObjCoeff, obj, _, 0.0, Handle) 1777 ), 1778 % process option list and fill in defaults 1779 process_options(OptionList, Handle, Module, SepGoal, EplexOptions), 1780 fill_in_defaults(Handle), 1781 arg(1, SolveSubProblem, SPHandle), 1782 fill_in_defaults_sp(Pool, SPHandle, Module), 1783 Handle = cg_prob{master_prob:MPHandle, 1784 module:Module, 1785 const_obj:ConstTerm, 1786 sp_solution_call:SolveSubProblem, 1787 phase1_obj:Phase1Obj, 1788 upper_bound:1.0Inf, 1789 lower_bound: -1.0Inf, 1790 mp_vars:Vars, 1791 stab_terms:StabTerms, 1792 stab_iter_counts:StabIterCounts, 1793 sep_call:(call(SepGoal)@Module), 1794 mp_cols_added:ColsAdded, 1795 pool:Pool, 1796 idx_lookup:Lookup, 1797 shelf:Shelf, 1798 phase:0}, 1799 StabIterCounts = stab_counters{bound_counter:1, 1800 coeff_counter:1}, 1801 hash_create(Lookup), 1802 verify(nonvar(Shelf)), 1803 % create the eplex handle and setup the fixed part of the obj 1804 ( NormObj == [] -> 1805 % need a dummy var in the obj fn to force 1806 % a CPLEX handle to be created properly 1807 % (or a dummy constraint) 1808 get_cg_attr(Dummy, Handle, DummyAttr), 1809 setarg(cost of colgen, DummyAttr, 0), 1810 setarg(coeffs of colgen, DummyAttr, []), 1811 setarg(aux of colgen, DummyAttr, artificial), 1812 MPObj=[ConstTerm, 1*Dummy] 1813 ; 1814 MPObj=[ConstTerm|NormObj] 1815 ), 1816 lp_setup([], min(sum(MPObj)), 1817 EplexOptions, 1818 MPHandle), 1819 lp_set(MPHandle, dual_solution, yes), 1820 lp_set(MPHandle, keep_basis, yes), 1821 lp_set(MPHandle, reduced_cost, yes), 1822 lp_set(MPHandle, slack, yes), 1823 % add any constraints only involving known MP vars: 1824 % if basis stabilisation is being used we will need to add 1825 % plus/minus stabilisation variables to each constraint, 1826 % we give them some default bounds of 0.0..1.0e+3 and 1827 % objective coefficient 0 initially 1828 % actually we always add them even if it is not used, to 1829 % make it easier to re-solve with it enabled if necessary 1830 % 1831 % note: the initial upper bound of 1.0e+3 is a bit arbitrary. 1832 % Adding these variables to the problem gives us a relaxation 1833 % of the problem. If the upper bounds are inifinite the 1834 % relaxation is always feasible, while fixing them to zero 1835 % gives the relaxation identical to the original problem. 1836 % Their intended purpose is to avoid the wildly oscillating 1837 % optimal dual values often encountered in column generation 1838 % we want is to solve 1839 ( 1840 foreach(MPCstr, MPCstrs), 1841 foreach(MPCstr1, MPNormCstrs), 1842 fromto(StabTerms, [StabTerm|Rest], 1843 Rest, StabTerms0), 1844 foreach(MPIdx, MPIdxs), 1845 param(Handle) 1846 do 1847 MPCstr = Type:[ConstTerm|LinTerms], 1848 get_cg_attr(Yminus, Handle, AttrM), 1849 setarg(cost of colgen, AttrM, 0), 1850 setarg(coeffs of colgen, AttrM, [Idx - -1]), 1851 setarg(aux of colgen, AttrM, stabilisation), 1852 get_cg_attr(Yplus, Handle, AttrP), 1853 setarg(cost of colgen, AttrP, 0), 1854 setarg(coeffs of colgen, AttrP, [Idx - 1]), 1855 setarg(aux of colgen, AttrP, stabilisation), 1856 MPCstr1 = Type:[ConstTerm, -1*Yminus, 1*Yplus|LinTerms], 1857 StabTerm = stab_term{idx: MPIdx, 1858 plus_var: Yplus, 1859 plus_coeff: 0, 1860 plus_bound: 1.0e+3, 1861 minus_var: Yminus, 1862 minus_coeff: 0, 1863 minus_bound: 1.0e+3} 1864 ), 1865 lp_add_indexed(MPHandle, MPNormCstrs, [], MPIdxs), 1866 % now add the constraints which involve generated vars and 1867 % setup the dual_var attributes of the subproblem vars with 1868 % index of MP constraint 1869 % note: we add in artificial variables in case the first 1870 % restricted MP at any node is infeasible 1871 % again we add plus/minus stabilisation variables to each 1872 % constraint and give them some default bounds of 0.0..1.0e+3 1873 % and objective coefficient 0 initially 1874 ( 1875 foreach(MPSPCstr, MPSPCstrs), 1876 foreach(StabTerm, StabTerms0), 1877 foreach(MPSPCstr1, MPSPNormCstrs), 1878 fromto(ArtVars, AVOut, AVIn, []), 1879 foreach(Idx, MPSPIdxs), 1880 param(Handle) 1881 do 1882 StabTerm = stab_term{idx:Idx, 1883 plus_var: Yplus, 1884 plus_coeff: 0, 1885 plus_bound: 1.0e+3, 1886 minus_var: Yminus, 1887 minus_coeff: 0, 1888 minus_bound: 1.0e+3}, 1889 MPSPCstr = Type:[ConstTerm|LinTerms], 1890 ( Type == (=<) -> 1891 get_cg_attr(Art, Handle, Attr), 1892 setarg(cost of colgen, Attr, 0), 1893 setarg(coeffs of colgen, Attr, [Idx - -1]), 1894 setarg(aux of colgen, Attr, artificial), 1895 get_cg_attr(Yminus, Handle, AttrM), 1896 setarg(cost of colgen, AttrM, 0), 1897 setarg(coeffs of colgen, AttrM, [Idx - -1]), 1898 setarg(aux of colgen, AttrM, stabilisation), 1899 get_cg_attr(Yplus, Handle, AttrP), 1900 setarg(cost of colgen, AttrP, 0), 1901 setarg(coeffs of colgen, AttrP, [Idx - 1]), 1902 setarg(aux of colgen, AttrP, stabilisation), 1903 MPSPCstr1 = Type:[ConstTerm, -1*Art, 1904 -1*Yminus, 1*Yplus|LinTerms], 1905 AVOut = [Art|AVIn] 1906 ; Type == (>=) -> 1907 get_cg_attr(Art, Handle, Attr), 1908 setarg(cost of colgen, Attr, 0), 1909 setarg(coeffs of colgen, Attr, [Idx - 1]), 1910 setarg(aux of colgen, Attr, artificial), 1911 get_cg_attr(Yminus, Handle, AttrM), 1912 setarg(cost of colgen, AttrM, 0), 1913 setarg(coeffs of colgen, AttrM, [Idx - -1]), 1914 setarg(aux of colgen, AttrM, stabilisation), 1915 get_cg_attr(Yplus, Handle, AttrP), 1916 setarg(cost of colgen, AttrP, 0), 1917 setarg(coeffs of colgen, AttrP, [Idx - 1]), 1918 setarg(aux of colgen, AttrP, stabilisation), 1919 MPSPCstr1 = Type:[ConstTerm, 1*Art, 1920 -1*Yminus, 1*Yplus|LinTerms], 1921 AVOut = [Art|AVIn] 1922 ; Type == (=:=) -> 1923 get_cg_attr(Art1, Handle, Attr1), 1924 setarg(cost of colgen, Attr1, 0), 1925 setarg(coeffs of colgen, Attr1, [Idx - -1]), 1926 setarg(aux of colgen, Attr1, artificial), 1927 get_cg_attr(Art2, Handle, Attr2), 1928 setarg(cost of colgen, Attr2, 0), 1929 setarg(coeffs of colgen, Attr2, [Idx - 1]), 1930 setarg(aux of colgen, Attr2, artificial), 1931 get_cg_attr(Yminus, Handle, AttrM), 1932 setarg(cost of colgen, AttrM, 0), 1933 setarg(coeffs of colgen, AttrM, [Idx - -1]), 1934 setarg(aux of colgen, AttrM, stabilisation), 1935 get_cg_attr(Yplus, Handle, AttrP), 1936 setarg(cost of colgen, AttrP, 0), 1937 setarg(coeffs of colgen, AttrP, [Idx - 1]), 1938 setarg(aux of colgen, AttrP, stabilisation), 1939 MPSPCstr1 = Type:[ConstTerm, -1*Art1, 1*Art2, 1940 -1*Yminus, 1*Yplus|LinTerms], 1941 AVOut = [Art1, Art2|AVIn] 1942 ) 1943 ), 1944 Phase1Obj = min(sum(ArtVars)), 1945 lp_add_indexed(MPHandle, MPSPNormCstrs, [], MPSPIdxs), 1946 expand_sp_obj_terms(MPIdxs, MPSPIdxs, CoeffVars, Handle), 1947 % give the known MP vars a colgen attribute 1948 % and set their objective cost 1949 % this has to be done after both types of constraints 1950 % above are added in case there are known MP vars 1951 % appearing in the type 2 constraints that do not 1952 % appear in the type 1 constraints 1953 % however, now the artificial variables have been added 1954 % to the MP so we have to check if each var has an attribute 1955 % already and only add attribute/include in vars of cg_prob 1956 % if it does not 1957 lp_get(MPHandle, vars, VarArr), 1958 ( 1959 foreacharg(Var, VarArr), 1960 fromto(Vars, Out, In, []), 1961 param(MPHandle, Handle) 1962 do 1963 get_cg_attr(Var, Handle, Attr), 1964 Attr = colgen{aux:Aux, lo:Lo, hi:Hi}, 1965 % set initial mp eplex bounds 1966 lp_var_set_bounds(MPHandle, Var, Lo, Hi), 1967 ( Aux == artificial -> 1968 % artificial variable 1969 Out = In 1970 ; Aux == stabilisation -> 1971 % stabilisation variable 1972 Out = In 1973 ; 1974 % known MP variable 1975 Attr = colgen{cost:0, coeffs:[], aux:[]}, 1976 Out = [Var|In] 1977 ) 1978 ), 1979 ( 1980 foreach(Const*Var, NormObj), 1981 param(Handle) 1982 do 1983 get_cg_attr(Var, Handle, Attr), 1984 setarg(cost of colgen, Attr, Const) 1985 ), 1986 % finally add initial SP solution column set to MP 1987 cg_new_MP_columns(Handle, VarCols), 1988 ( VarCols == [] -> 1989 % no initial solution columns 1990 true 1991 ; 1992 % add initial solution columns to MP 1993 lp_add_columns(MPHandle, VarCols), 1994 % record the total number of columns now in the mp 1995 lp_get(MPHandle, vars, MPVarArr), 1996 arity(MPVarArr, ColsAdded) 1997 ), 1998 ( 1999 foreach(Var, ArtVars) 2000 do 2001 set_var_bounds(Var, 0, 1.0Inf) 2002 ), 2003 ( 2004 foreach(StabTerm, StabTerms), 2005 param(MPHandle) 2006 do 2007 StabTerm = stab_term{plus_var: Yplus, 2008 plus_bound: Boundplus, 2009 minus_var: Yminus, 2010 minus_bound: Boundminus}, 2011 lp_var_non_monotonic_set_bounds(MPHandle, Yplus, 0, Boundplus), 2012 lp_var_non_monotonic_set_bounds(MPHandle, Yminus, 0, Boundminus) 2013 ), 2014 ( NormObj == [] -> Dummy = 0 ; true ), 2015 % make a suspension for the MP iterator 2016 % and insert it in the suspension list of Handle 2017 make_suspension(cg_iteration(Handle), 0, MPSusp), 2018 enter_suspension_list(mp_susp of cg_prob, Handle, MPSusp), 2019 % make a suspension for the SP iterator 2020 % and insert it in the suspension lists of the dual val vars 2021 make_suspension(solveSPs(Handle), 0, SPSusp), 2022 SPHandle = sp_prob{cost:OptVar, coeff_vars:DualVars}, 2023 insert_suspension([OptVar|DualVars], SPSusp, susps of dual_var, 2024 dual_var). 2025 2026add_new_solver_rowcols(Handle) :- 2027 Handle = cg_prob{master_prob:MPHandle, 2028 phase1_obj:Phase1Obj, 2029 mp_vars:OldVars, 2030 stab_terms:OldStabTerms, 2031 pool:Pool}, 2032 ( var(Pool) -> 2033 % there is no associated colgen pool 2034 MPIdxCstrs = [], 2035 MPSPNormCstrs = [] 2036 ; 2037 % collect any new constraints from the pool 2038 collect_typed_pool_constraints(Pool, mp_only of cg_constraint_type, 2039 MPIdxCstrs), 2040 collect_typed_pool_constraints(Pool, mp_sp of cg_constraint_type, 2041 MPSPNormCstrs) 2042 ), 2043 ( ( MPIdxCstrs == [], MPSPNormCstrs == [] ) -> 2044 % there were no new constraints, 2045 % done 2046 true 2047 ; 2048 % there were some, add them to the MP solver 2049 get_pool_handle(Handle, Pool), 2050 % process any constraints only involving known MP vars 2051 ( 2052 foreach(MPIdx:MPIdxCstr, MPIdxCstrs), 2053 fromto(NormExprs, [Expr|RestExprs], RestExprs, NormExprs0), 2054 foreach(MPIdxCstr1, MPIdxNormCstrs), 2055 fromto(StabTerms, [StabTerm|Rest], 2056 Rest, StabTerms0), 2057 fromto(NewStabTerms, [StabTerm|Rest0], 2058 Rest0, StabTerms00), 2059 foreach(MPIdx, MPIdxs), 2060 param(Handle) 2061 do 2062 % add a stabilisation term to the constraint before 2063 % sending to MPHandle 2064 MPIdxCstr = Type:Expr, 2065 MPIdxCstr1 = Type:Expr1, 2066 add_stabilisation_term(MPIdx, Handle, Expr, Expr1, StabTerm) 2067 ), 2068 % now process the constraints which involve generated vars and 2069 % setup the dual_var attributes of the subproblem vars with 2070 % index of MP constraint 2071 % note: we add in artificial variables in case the first 2072 % restricted MP at any node is infeasible 2073 Phase1Obj = min(sum(OldArtVars)), 2074 ( 2075 foreach([CoeffVar, Idx]:NormCstr, MPSPNormCstrs), 2076 foreach(Expr, NormExprs0), 2077 fromto(StabTerms0, [StabTerm|Rest], 2078 Rest, OldStabTerms), 2079 foreach(StabTerm, StabTerms00), 2080 foreach(CoeffVar, CoeffVars), 2081 foreach(Cstr, Cstrs), 2082 fromto(NewArtVars, NAVOut, NAVIn, []), 2083 fromto(ArtVars, AVOut, AVIn, OldArtVars), 2084 foreach(Idx, Idxs), 2085 param(Handle) 2086 do 2087 % add a stabilisation term and artificial vars to the 2088 % constraint before sending to MPHandle 2089 NormCstr = Type:Expr, 2090 add_stabilisation_term(Idx, Handle, Expr, [ConstTerm|LinTerms], StabTerm), 2091 ( Type == (=<) -> 2092 get_cg_attr(Art, Handle, Attr), 2093 setarg(cost of colgen, Attr, 0), 2094 setarg(coeffs of colgen, Attr, [Idx - -1]), 2095 setarg(aux of colgen, Attr, artificial), 2096 Cstr = Type:[ConstTerm, -1*Art|LinTerms], 2097 NAVOut = [Art|NAVIn], 2098 AVOut = [Art|AVIn] 2099 ; Type == (>=) -> 2100 get_cg_attr(Art, Handle, Attr), 2101 setarg(cost of colgen, Attr, 0), 2102 setarg(coeffs of colgen, Attr, [Idx - 1]), 2103 setarg(aux of colgen, Attr, artificial), 2104 Cstr = Type:[ConstTerm, 1*Art|LinTerms], 2105 NAVOut = [Art|NAVIn], 2106 AVOut = [Art|AVIn] 2107 ; Type == (=:=) -> 2108 get_cg_attr(Art1, Handle, Attr1), 2109 setarg(cost of colgen, Attr1, 0), 2110 setarg(coeffs of colgen, Attr1, [Idx - -1]), 2111 setarg(aux of colgen, Attr1, artificial), 2112 get_cg_attr(Art2, Handle, Attr2), 2113 setarg(cost of colgen, Attr2, 0), 2114 setarg(coeffs of colgen, Attr2, [Idx - 1]), 2115 setarg(aux of colgen, Attr2, artificial), 2116 Cstr = Type:[ConstTerm, -1*Art1, 1*Art2|LinTerms], 2117 NAVOut = [Art1, Art2|NAVIn], 2118 AVOut = [Art1, Art2|AVIn] 2119 ) 2120 ), 2121 NewPhase1Obj = min(sum(ArtVars)), 2122 setarg(phase1_obj of cg_prob, Handle, NewPhase1Obj), 2123 setarg(stab_terms of cg_prob, Handle, StabTerms), 2124 % give the new MP variables a colgen attribute 2125 % this has to be done for variables in both types of constraints 2126 % however, some variables may already have been added 2127 % to the MP so we have to check if each variable occurs 2128 % already and only add attribute/include in vars of cg_prob 2129 % if it does not 2130 term_variables(NormExprs, ExprVars), 2131 ( 2132 foreach(ExprVar, ExprVars), 2133 fromto(Vars, Out, In, OldVars), 2134 param(MPHandle, Handle) 2135 do 2136 ( lp_var_occurrence(ExprVar, MPHandle, _) -> 2137 % already added to MPHandle, 2138 % it must be a pre-existing variable for Handle 2139 Out = In 2140 ; 2141 % not yet in MPHandle, 2142 % it is a new variable for Handle 2143 get_cg_attr(ExprVar, Handle, colgen{lo:Lo, hi:Hi}), 2144 % set initial mp eplex bounds 2145 lp_var_set_bounds(MPHandle, ExprVar, Lo, Hi), 2146 Out = [ExprVar|In] 2147 ) 2148 ), 2149 lp_add_indexed(MPHandle, MPIdxNormCstrs, [], MPIdxs), 2150 lp_add_indexed(MPHandle, Cstrs, [], Idxs), 2151 expand_sp_obj_terms(MPIdxs, Idxs, CoeffVars, Handle), 2152 % finally add any new SP solution columns to MP 2153 cg_new_MP_columns(Handle, VarCols), 2154 ( VarCols == [] -> 2155 % no initial solution columns 2156 true 2157 ; 2158 % add initial solution columns to MP 2159 lp_add_columns(MPHandle, VarCols), 2160 ( 2161 foreach(Var:_, VarCols), 2162 fromto(NewVars, [Var|Rest], Rest, Vars) 2163 do 2164 true 2165 ), 2166 % record the total number of columns now in the mp 2167 lp_get(MPHandle, vars, MPVarArr), 2168 arity(MPVarArr, NewColsAdded), 2169 setarg(mp_vars of cg_prob, Handle, NewVars), 2170 setarg(mp_cols_added of cg_prob, Handle, NewColsAdded) 2171 ), 2172 ( 2173 foreach(Var, NewArtVars) 2174 do 2175 set_var_bounds(Var, 0, 1.0Inf) 2176 ), 2177 ( 2178 foreach(StabTerm, NewStabTerms), 2179 param(MPHandle) 2180 do 2181 StabTerm = stab_term{plus_var: Yplus, 2182 plus_bound: Boundplus, 2183 minus_var: Yminus, 2184 minus_bound: Boundminus}, 2185 lp_var_non_monotonic_set_bounds(MPHandle, Yplus, 0, Boundplus), 2186 lp_var_non_monotonic_set_bounds(MPHandle, Yminus, 0, Boundminus) 2187 ), 2188 % insert the SP iterator suspension in the suspension lists of 2189 % the dual val vars 2190 % TODO: this looks wrong - why a new suspension??? 2191 % also: sp_prob's coeff_vars field hasn't been updated? 2192 make_suspension(solveSPs(Handle), 0, SPSusp), 2193 insert_suspension(CoeffVars, SPSusp, susps of dual_var, 2194 dual_var) 2195 ). 2196 2197add_stabilisation_term(Idx, Handle, Expr, StabExpr, StabTerm) :- 2198 Expr = [ConstTerm|LinTerms], 2199 % add stabilisation var terms to the linear expression 2200 StabExpr = [ConstTerm, -1*Yminus, 1*Yplus|LinTerms], 2201 % create a new stab_term structure for them 2202 StabTerm = stab_term{idx: Idx, 2203 plus_var: Yplus, 2204 plus_coeff: 0, 2205 plus_bound: 1.0e+3, 2206 minus_var: Yminus, 2207 minus_coeff: 0, 2208 minus_bound: 1.0e+3}, 2209 % add colgen attributes to the new stabilisation vars 2210 get_cg_attr(Yminus, Handle, AttrM), 2211 setarg(cost of colgen, AttrM, 0), 2212 setarg(coeffs of colgen, AttrM, [Idx - -1]), 2213 setarg(aux of colgen, AttrM, stabilisation), 2214 get_cg_attr(Yplus, Handle, AttrP), 2215 setarg(cost of colgen, AttrP, 0), 2216 setarg(coeffs of colgen, AttrP, [Idx - 1]), 2217 setarg(aux of colgen, AttrP, stabilisation). 2218 2219expand_sp_obj_terms([], [], [], _Handle) :- !. 2220expand_sp_obj_terms(MPIdxs, Idxs, CoeffVars, Handle) :- 2221 ( MPIdxs == [] -> 2222 NCstrs is max(Idxs) 2223 ; Idxs == [] -> 2224 NCstrs is max(MPIdxs) 2225 ; 2226 NCstrs is max(max(MPIdxs), max(Idxs)) 2227 ), 2228 NCstrs1 is NCstrs + 1, 2229 dim(ExprArr, [NCstrs1]), 2230 Handle = cg_prob{sp_obj_terms:OldExprArr}, 2231 ( var(OldExprArr) -> 2232 OldExprArr = ExprArr 2233 ; 2234 dim(OldExprArr, [OldNCstrs]), 2235 NCstrs1 > OldNCstrs, 2236 ( 2237 for(Idx, 1, OldNCstrs), 2238 param(OldExprArr, ExprArr) 2239 do 2240 arg(Idx, OldExprArr, Arg), 2241 arg(Idx, ExprArr, Arg) 2242 ), 2243 setarg(sp_obj_terms of cg_prob, Handle, ExprArr) 2244 ), 2245 ( 2246 foreach(CoeffVar, CoeffVars), 2247 foreach(Idx, Idxs), 2248 param(ExprArr) 2249 do 2250 Idx1 is Idx + 1, 2251 arg(Idx1, ExprArr, CoeffVar) 2252 ), 2253 ( 2254 foreach(MPIdx, MPIdxs), 2255 param(ExprArr) 2256 do 2257 MPIdx1 is MPIdx + 1, 2258 arg(MPIdx1, ExprArr, 0) 2259 ). 2260 2261process_options([], _, _, Separation, []) :- !, 2262 ( var(Separation) -> Separation = true ; true ). 2263process_options([O|Os], Handle, Module, Separation, EplexOptions) :- !, 2264 process_option(O, Handle, Module, Separation, EplexOptions, EplexOptions0), 2265 process_options(Os, Handle, Module, Separation, EplexOptions0). 2266process_options(_NonList, Handle, _, _, _) :- 2267 Handle = cg_prob{pool:Pool}, 2268 cg_info_message(Handle, "%w : options not proper list." 2269 " Ignored.%n,", [Pool]). 2270 2271process_option(separate(Separation0), Handle, Module, Separation, EplexOptions, EplexOptions0) ?- !, 2272 ( var(Separation0) -> 2273 error(4, process_option(separate(Separation0), Handle, Module, Separation, 2274 EplexOptions, EplexOptions0)) 2275 ; 2276 Separation = Separation0, 2277 EplexOptions = EplexOptions0 2278 ). 2279process_option(node_select(Val), Handle, Module, Separation, EplexOptions, EplexOptions0) ?- !, 2280 ( var(Val) -> 2281 error(4, process_option(node_select(Val), Handle, Module, Separation, 2282 EplexOptions, EplexOptions0)) 2283 ; 2284 ( Handle = cg_prob{bfs_tree:[], info_messages:OnOff} -> 2285 bfs_solver_setup(BfsHandle, min, bp_node(Handle), 2286 [info_messages(OnOff), 2287 node_select(Val), 2288 separation(bp_separate(Handle))]), 2289 setarg(bfs_tree of cg_prob, Handle, BfsHandle) 2290 ; 2291 Handle = cg_prob{bfs_tree:BfsHandle}, 2292 bfs_set(BfsHandle, node_select, Val) 2293 ), 2294 EplexOptions = EplexOptions0 2295 ). 2296process_option(int_tolerance(Val), Handle, Module, Separation, EplexOptions, EplexOptions0) ?- !, 2297 ( var(Val) -> 2298 error(4, process_option(int_tolerance(Val), Handle, Module, Separation, 2299 EplexOptions, EplexOptions0)) 2300 ; float(Val) -> 2301 ( 0 =< Val, Val =< 0.5 -> 2302 cg_set(Handle, int_tolerance, Val), 2303 EplexOptions = EplexOptions0 2304 ; error(6, process_option(int_tolerance(Val), Handle, Module, Separation, 2305 EplexOptions, EplexOptions0)) 2306 ) 2307 ; error(5, process_option(int_tolerance(Val), Handle, Module, Separation, 2308 EplexOptions, EplexOptions0)) 2309 ). 2310process_option(disallow(Val), Handle, Module, Separation, EplexOptions, EplexOptions0) ?- !, 2311 ( var(Val) -> 2312 error(4, process_option(disallow(Val), Handle, Module, Separation, 2313 EplexOptions, EplexOptions0)) 2314 ; atom(Val) -> 2315 ( valid_disallow_setting(Val) -> 2316 cg_set(Handle, disallow, Val), 2317 EplexOptions = EplexOptions0 2318 ; error(6, process_option(disallow(Val), Handle, Module, Separation, 2319 EplexOptions, EplexOptions0)) 2320 ) 2321 ; error(5, process_option(disallow(Val), Handle, Module, Separation, 2322 EplexOptions, EplexOptions0)) 2323 ). 2324process_option(info_messages(Val), Handle, Module, Separation, EplexOptions, EplexOptions0) ?- !, 2325 ( var(Val) -> 2326 error(4, process_option(info_messages(Val), Handle, Module, Separation, 2327 EplexOptions, EplexOptions0)) 2328 ; atom(Val) -> 2329 ( valid_setting(Val) -> 2330 cg_set(Handle, info_messages, Val), 2331 EplexOptions = EplexOptions0 2332 ; error(6, process_option(info_messages(Val), Handle, Module, Separation, 2333 EplexOptions, EplexOptions0)) 2334 ) 2335 ; error(5, process_option(info_messages(Val), Handle, Module, Separation, 2336 EplexOptions, EplexOptions0)) 2337 ). 2338process_option(on_degeneracy(Val), Handle, Module, Separation, EplexOptions, EplexOptions0) ?- !, 2339 ( var(Val) -> 2340 error(4, process_option(on_degeneracy(Val), Handle, Module, Separation, 2341 EplexOptions, EplexOptions0)) 2342 ; atom(Val) -> 2343 ( valid_degeneracy_setting(Val) -> 2344 cg_set(Handle, on_degeneracy, Val), 2345 EplexOptions = EplexOptions0 2346 ; error(6, process_option(on_degeneracy(Val), Handle, Module, Separation, 2347 EplexOptions, EplexOptions0)) 2348 ) 2349 ; error(5, process_option(on_degeneracy(Val), Handle, Module, Separation, 2350 EplexOptions, EplexOptions0)) 2351 ). 2352process_option(stabilisation(Val), Handle, Module, Separation, EplexOptions, EplexOptions0) ?- !, 2353 ( var(Val) -> 2354 error(4, process_option(stabilisation(Val), Handle, Module, Separation, 2355 EplexOptions, EplexOptions0)) 2356 ; valid_stabilisation_setting(Val) -> 2357 cg_set(Handle, stabilisation, Val), 2358 EplexOptions = EplexOptions0 2359 ; error(6, process_option(stabilisation(Val), Handle, Module, Separation, 2360 EplexOptions, EplexOptions0)) 2361 ). 2362process_option(basis_perturbation(Val), Handle, Module, Separation, EplexOptions, EplexOptions0) ?- !, 2363 ( var(Val) -> 2364 error(4, process_option(basis_perturbation(Val), Handle, Module, Separation, 2365 EplexOptions, EplexOptions0)) 2366 ; atom(Val) -> 2367 ( valid_setting(Val) -> 2368 cg_set(Handle, basis_perturbation, Val), 2369 EplexOptions = EplexOptions0 2370 ; error(6, process_option(basis_perturbation(Val), Handle, Module, Separation, 2371 EplexOptions, EplexOptions0)) 2372 ) 2373 ; error(5, process_option(basis_perturbation(Val), Handle, Module, Separation, 2374 EplexOptions, EplexOptions0)) 2375 ). 2376process_option(eplex_option(Val), Handle, Module, Separation, EplexOptions, EplexOptions0) ?- !, 2377 ( var(Val) -> 2378 error(4, process_option(eplex_option(Val), Handle, Module, Separation, 2379 EplexOptions, EplexOptions0)) 2380 ; 2381 EplexOptions = [Val|EplexOptions0] 2382 ). 2383process_option(Option, Handle, Module, Separation, EplexOptions, EplexOptions0) :- 2384 ( var(Option) -> 2385 error(4, process_option(Option, Handle, Module, Separation, 2386 EplexOptions, EplexOptions0)) 2387 ; 2388 error(6, process_option(Option, Handle, Module, Separation, 2389 EplexOptions, EplexOptions0)) 2390 ). 2391 2392valid_degeneracy_setting(continue). 2393valid_degeneracy_setting(stop). 2394 2395valid_stabilisation_setting(on(BoundIter, BoundUpdate, 2396 CoeffIter, CoeffUpdate)) ?- 2397 integer(BoundIter), 2398 BoundIter >= 1, 2399 number(BoundUpdate), 2400 BoundUpdate > 0, 2401 integer(CoeffIter), 2402 CoeffIter >= 1, 2403 number(CoeffUpdate), 2404 CoeffUpdate > 0. 2405valid_stabilisation_setting(off). 2406valid_stabilisation_setting(stab_pred(_, _)). 2407 2408valid_disallow_setting(clp). 2409valid_disallow_setting(lp). 2410valid_disallow_setting(off). 2411 2412valid_setting(on). 2413valid_setting(off). 2414 2415fill_in_defaults(Handle) :- 2416 Handle = cg_prob{disallow:DisOnOff, 2417 tolerance:Tol, 2418 branch_tolerance:BranchTol, 2419 info_messages:MessageOnOff, 2420 on_degeneracy:OnDegeneracy, 2421 stabilisation:Stabilisation, 2422 basis_perturbation:PerturbOnOff}, 2423 ( var(DisOnOff) -> DisOnOff = off ; true ), 2424 ( var(Tol) -> Tol = 1e-05 ; true ), 2425 ( var(BranchTol) -> BranchTol = 0 ; true ), 2426 ( var(MessageOnOff) -> MessageOnOff = off ; true ), 2427 ( var(OnDegeneracy) -> OnDegeneracy = stop ; true ), 2428 ( var(Stabilisation) -> Stabilisation = off ; true ), 2429 ( var(PerturbOnOff) -> PerturbOnOff = off ; true ). 2430 2431fill_in_defaults_sp(Pool, SPHandle, Module) :- 2432 SPHandle = sp_prob{master_pool:Pool, 2433 cutoff:CutOff, 2434 disallow:Disallow, 2435 status: phase1, 2436 module:Module, 2437 lo:Lo, hi:Hi, type:Type 2438 }, 2439 ( var(CutOff) -> CutOff = 1e-05 ; true ), 2440 ( var(Disallow) -> Disallow = [0, []] ; true ), 2441 ( var(Lo) -> Lo = 0.0 ; true ), % or -1.0Inf, but 0 more likely 2442 ( var(Hi) -> Hi = 1.0Inf ; true ), 2443 ( var(Type) -> Type = real ; true ). 2444 2445 2446% ---------------------------------------------------------------------- 2447% most general instantiation templates to disallow duplicate columns 2448% ---------------------------------------------------------------------- 2449 2450% repeatedly amalgamate a new coefficient list with 2451% the coefficient list of a member of a list of 2452% linearisation-coefficient list pairs, extract 2453% linearisation of new most general amalgamation 2454% and return new linearisation-coefficient list pair list 2455lin_amalgamate(N, List, Coeffs, Vars, NewN, [[N1, Cstrs]-NewCoeffs|NewList]) :- 2456 % try to amalgamate Coeffs with others in list 2457 % repeatedly until we have most general disallow template 2458 ( 2459 fromto(no, _, Done, yes), 2460 fromto(Coeffs, CIn, COut, NewCoeffs), 2461 fromto(N, NIn, NOut, N0), 2462 fromto(List, LIn, LOut, NewList) 2463 do 2464 lin_amalgamate(LIn, CIn, NIn, LOut, COut, NOut, Done) 2465 ), 2466 % linearise new template 2467 linearize_template(Vars, NewCoeffs, 1, N1, LCstrs, Flags, Done), 2468 ( Done == yes -> NewN = N0, N1 = 0, Cstrs = [] 2469 ; NewN is N0+N1, Cstrs = [(>=):[-1*1|Flags]|LCstrs] ). 2470 2471linearize_template([], [], N, N, [], [], no). 2472linearize_template([Var|Vars], [Coeff|Coeffs], N0, N, Cstrs, Flags, Done) :- 2473 get_var_bounds(Var, Lo, Hi), 2474 ( Coeff = [_|_] -> CList = Coeff ; CList = [Coeff] ), 2475 lin_not_in(CList, Var, Lo, Hi, N0, N1, Cstrs, Cstrs1, Flags, Flags1, Done), 2476 ( Done == yes -> % Var cannot match, no Cstrs needed 2477 true 2478 ; 2479 linearize_template(Vars, Coeffs, N1, N, Cstrs1, Flags1, Done) 2480 ). 2481 2482lin_not_in([], _Var, _VLo, _VHi, _N, _N1, _Cstrs, _Cstrs1, _Flags, _Flags1, yes). 2483lin_not_in([H|T], Var, VLo, VHi, N, N1, Cstrs, Cstrs1, Flags, Flags1, Done) :- 2484 ( integer(H) -> Lo = H, Hi = H ; H = Lo..Hi ), 2485 ( VHi < Lo -> % cannot be in this or later ranges 2486 Done = yes 2487 ; VLo < Lo, VHi =< Hi -> % may be in this range, cannot be in later 2488 % Var =< VHi-(VHi+1-Lo)*Flag 2489 Alpha is fix(VHi+1-Lo), 2490 Beta is fix(-1*VHi), 2491 N1 is N+1, 2492 Cstrs = [(=<):[Beta*1, 1*Var, Alpha*Flag]|Cstrs1], 2493 Flags = [1*Flag|Flags1], 2494 Flag::0..1 2495 ; VHi =< Hi -> % is in this range, no cstr for this var 2496 N = N1, 2497 Cstrs = Cstrs1, 2498 Flags = Flags1 2499 ; VLo =< VHi -> % may be in this or later but cannot be earlier 2500 lin_not_in1(T, Hi, Var, VLo, VHi, N, N1, Cstrs, Cstrs1, Flags, Flags1) 2501 ; % may be in later ranges but cannot be in this 2502 lin_not_in(T, Var, VLo, VHi, N, N1, Cstrs, Cstrs1, Flags, Flags1, Done) 2503 ). 2504 2505lin_not_in1([], Hi, Var, VLo, _VHi, N, N1, 2506 [(>=):[Delta*1,1*Var,Gamma*Flag]|Cstrs1], Cstrs1, 2507 [1*Flag|Flags1], Flags1) :- 2508 % Var >= VLo+(Hi+1-VLo)*Flag 2509 N1 is N+1, 2510 Flag::0..1, 2511 Gamma is fix(VLo-1-Hi), 2512 Delta is fix(-1*VLo). 2513lin_not_in1([H|T], Hi, Var, VLo, VHi, N, N1, 2514 [(>=):[Delta*1,1*Var,Gamma*Flag], 2515 (=<):[Beta*1,1*Var,Alpha*Flag]|Cstrs], Cstrs1, 2516 [1*Flag|Flags], Flags1) :- 2517 % Var >= VLo+(Hi+1-VLo)*Flag 2518 % Var =< VHi-(VHi+1-Lo1)*Flag 2519 ( integer(H) -> Lo1 = H, Hi1 = H ; H = Lo1..Hi1 ), 2520 N0 is N+2, 2521 Flag::0..1, 2522 Alpha is fix(VHi+1-Lo1), 2523 Beta is fix(-1*VHi), 2524 Gamma is fix(VLo-1-Hi), 2525 Delta is fix(-1*VLo), 2526 lin_not_in1(T, Hi1, Var, VLo, VHi, N0, N1, Cstrs, Cstrs1, Flags, Flags1). 2527 2528% amalgamate a coefficient list with the coefficient 2529% list of a member of a list of linearisation-coefficient 2530% list pairs and return remaining list and amalgamated coefficients 2531lin_amalgamate([], Coeffs, N, [], Coeffs, N, yes). 2532lin_amalgamate([[L, Cstrs]-Coeffs|List], NewCoeffs, N, 2533 NewList, NewCoeffs1, NewN, Done) :- 2534 ( amalg(NewCoeffs, Coeffs, 0, 0, NewCoeffs1) -> 2535 % new coeff list can be amalgamated with this one 2536 % to create a more general disallow template 2537 % return rest of list and new template 2538 NewList = List, 2539 NewN is N-L, 2540 Done = no 2541 ; 2542 % could not be amalgamated 2543 NewList = [[L, Cstrs]-Coeffs|NewList1], 2544 lin_amalgamate(List, NewCoeffs, N, NewList1, NewCoeffs1, NewN, Done) 2545 ). 2546 2547% repeatedly amlagamate a new coefficient list with 2548% the coefficient list of a member of a list of 2549% suspension-coefficient list pairs killing any 2550% associated suspensions, suspend disallow demons for 2551% the most general amalgamation and return new 2552% suspension-coefficient list pair list 2553% The coefficients are numbers or ../2 ranges. 2554amalgamate(N, List, Coeffs, Vars, NewN, [[N1, Susps]-NewCoeffs|NewList]) :- 2555 % try to amalgamate Coeffs with others in list 2556 % repeatedly until we have most general disallow template 2557 ( 2558 fromto(no, _, Done, yes), 2559 fromto(Coeffs, CIn, COut, NewCoeffs), 2560 fromto(N, NIn, NOut, N0), 2561 fromto(List, LIn, LOut, NewList) 2562 do 2563 amalgamate(LIn, CIn, NIn, LOut, COut, NOut, Done) 2564 ), 2565 % suspend disallow demons for new template 2566 ( 2567 foreach(Var, Vars), 2568 foreach(Coeff, NewCoeffs), 2569 foreach(Flag, AndFlags), 2570 foreach(Susp, AndSusps), 2571 count(_,2,NSusps), 2572 param(OrFlag) 2573 do 2574 ( Coeff = [_|_] -> ListTerm = l(Coeff) ; ListTerm = l([Coeff]) ), 2575 suspend(not_among(Var, ListTerm, Flag, OrFlag, Susp), 0, [Var->constrained, [Flag, OrFlag]->inst], Susp), 2576 not_among(Var, ListTerm, Flag, OrFlag, Susp) 2577 ), 2578 ( AndFlags = [AndFlag] -> 2579 AndFlag = 1, 2580 NewN is N0+1, 2581 N1 = 1, 2582 Susps = AndSusps 2583 ; 2584 NewN is N0+NSusps, 2585 N1 = NSusps, 2586 Susps = [OrSusp|AndSusps], 2587 AndFlagTerm = l(AndFlags), 2588 suspend(or(OrFlag, AndFlagTerm, OrSusp), 0, [OrFlag|AndFlags]->inst, OrSusp), 2589 or(OrFlag, AndFlagTerm, OrSusp) 2590 ). 2591 2592% amalgamate a coefficient list with the coefficient list of 2593% a member of a list of suspension-coefficient list pairs 2594% killing any associated suspensions and return remaining 2595% list and amalgamated coefficients 2596amalgamate([], Coeffs, N, [], Coeffs, N, yes). 2597amalgamate([[L, Susps]-Coeffs|List], NewCoeffs, N, 2598 NewList, NewCoeffs1, NewN, Done) :- 2599 ( amalg(NewCoeffs, Coeffs, 0, 0, NewCoeffs1) -> 2600 % new coeff list can be amalgamated with this one 2601 % to create a more general disallow template 2602 % kill the suspension for this disallow goal 2603 % return rest of list and new template 2604 ( 2605 foreach(Susp, Susps) 2606 do 2607 kill_suspension(Susp) 2608 ), 2609 NewList = List, 2610 NewN is N-L, 2611 Done = no 2612 ; 2613 % could not be amalgamated 2614 NewList = [[L, Susps]-Coeffs|NewList1], 2615 amalgamate(List, NewCoeffs, N, NewList1, NewCoeffs1, NewN, Done) 2616 ). 2617 2618% amalgamate two coefficient lists to create a _most general amalgamation_ 2619% succeed if: they are pairwise disjoint at exactly one coeff and 2620% pairwise indentical at all others 2621% abort if: they are equal/intersecting at all coeffs (since we have a 2622% duplicate coeff list which should have been disallowed) 2623% fail: otherwise 2624amalg([], [], D, I, []) :- 2625 ( D = 0 -> printf("lists not disjoint in amalg/5%n", []), 2626 flush(output), abort 2627 ; D = 1, I = 0 ). 2628amalg([A|As], [B|Bs], D, I, Amalg) :- 2629 ( A == B -> 2630 % pairwise identical here 2631 D1 = D, I1 = I, Amalg = [A|Amalg1] 2632 ; 2633 % otherwise make both coeff templates into 2634 % disallowed domain lists 2635 % and check for disjointness/intersection 2636 ( A = [_|_] -> List1 = A ; List1 = [A] ), 2637 ( B = [_|_] -> List2 = B ; List2 = [B] ), 2638 amalg(List1, List2, D, D1, I, I1, List3), Amalg = [List3|Amalg1] 2639 ), 2640 amalg(As, Bs, D1, I1, Amalg1). 2641 2642amalg(List1, List2, D, D1, I, I1, Amalg) :- 2643 amalg1(List1, List2, D, D1, I, I1, List3), 2644 ( List3 = [Val] -> Amalg = Val 2645 ; Amalg = List3 ). 2646 2647amalg1([], List, 0, 1, 0, 0, List) :- !. 2648amalg1(List, [], 0, 1, 0, 0, List) :- !. 2649amalg1([H1|T1], [H2|T2], D, D1, I, I1, Amalg) :- 2650 ( integer(H1) -> Lo1 = H1, Hi1 = H1 ; H1 = Lo1..Hi1 ), 2651 ( integer(H2) -> Lo2 = H2, Hi2 = H2 ; H2 = Lo2..Hi2 ), 2652 ( overlapping_ranges(Lo1, Hi1, Lo2, Hi2) -> 2653 % lists intersect here, 2654 % fail immediately if already disjoint 2655 % otherwise succeed with 2656 % no need to continue amalgamating 2657 D = 0, D1 = 0, I1 = 1 2658 ; 2659 % no intersection between H1 and H2, 2660 % continue with T1 and T2, 2661 % amalgamating H1 and H2 if necessary 2662 ( Lo1 is Hi2 + 1 -> T3 = [Lo2..Hi1|T1], T4 = T2, Amalg = Amalg1 2663 ; Lo1 > Hi2 + 1 -> T3 = [H1|T1], T4 = T2, Amalg = [H2|Amalg1] 2664 ; Lo2 is Hi1 + 1 -> T3 = T1, T4 = [Lo1..Hi2|T2], Amalg = Amalg1 2665 ; T3 = T1, T4 = [H2|T2], Amalg = [H1|Amalg1] ), 2666 amalg1(T3, T4, D, D1, I, I1, Amalg1) 2667 ). 2668 2669overlapping_ranges(Lo1, Hi1, Lo2, Hi2) :- 2670 (Lo1 >= Lo2, Lo1 =< Hi2) ; (Lo2 >= Lo1, Lo2 =< Hi1). 2671 2672% suspend this on [OrFlag|AndFlags]->inst 2673:- demon or/3. 2674:- set_flag(or/3, priority, 3). 2675or(OrFlag, AndFlagTerm, Susp) :- 2676 ( OrFlag == 1 -> 2677 kill_suspension(Susp) 2678 ; 2679 % an AndFlag was set to 0 2680 AndFlagTerm =.. [_, AndFlags], 2681 ( 2682 foreach(AndFlag, AndFlags), 2683 fromto(NewFlags, Out, In, []) 2684 do 2685 ( AndFlag == 0 -> Out = In ; Out = [AndFlag|In] ) 2686 ), 2687 ( NewFlags = [AndFlag] -> AndFlag = 1, kill_suspension(Susp) 2688 ; setarg(1, AndFlagTerm, NewFlags) ) 2689 ). 2690 2691% suspend this on Var->any, [Flag, OrFlag]->inst 2692:- demon not_among/5. 2693:- set_flag(not_among/5, priority, 4). 2694not_among(Var, ListTerm, Flag, OrFlag, Susp) :- 2695 ( OrFlag == 1 -> 2696 % some other variable has been 2697 % instantiated to a non-matching 2698 % value, done 2699 kill_suspension(Susp) 2700 ; 2701 % change in domain of Var 2702 % or all other variables have been 2703 % instantiated to matching values, 2704 % Var cannot match ListTerm 2705 get_var_bounds(Var, Lo, Hi), 2706 ListTerm =.. [_, List], 2707 not_among_body(List, Lo, Hi, NewList, Flag), 2708 ( Flag == 0 -> 2709 % matches the list, done 2710 kill_suspension(Susp) 2711 ; NewList == [] -> 2712 % does not match, done 2713 OrFlag = 1, Flag = 1, kill_suspension(Susp) 2714 ; Flag == 1 -> 2715 % all other variables have been 2716 % instantiated to matching values, 2717 % Var cannot match ListTerm 2718 2719 % is there a way to remove intervals from ic var domains? 2720 2721 % for now 2722 2723 % remove the top and bottom overlaps 2724 % if it can be done without creating holes 2725 NewList = [H|Rest], 2726 ( integer(H) -> Lo1 = H, Hi1 = H ; H = Lo1..Hi1 ), 2727 ( Lo1 =< Lo -> NewLo is Hi1 + 1 ; NewLo = Lo ), 2728 ( Rest == [] -> NewHi = Hi 2729 ; 2730 ( fromto(Rest, [_|Out], Out, [T]) do true ), 2731 ( integer(T) -> Lo2 = T, Hi2 = T ; T = Lo2..Hi2 ), 2732 ( Hi2 >= Hi -> NewHi is Lo2 - 1 ; NewHi = Hi ) 2733 ), 2734 set_var_bounds(Var, NewLo, NewHi), 2735 setarg(1, ListTerm, NewList) 2736 ; setarg(1, ListTerm, NewList) ) 2737 ). 2738 2739% not_among_body(+RangeList, +VLo, +VHi, ?List, ?Flag) 2740% succeeds if VLo..VHi is a subset of a member of RangeList and Flag = 0 2741% or if List is the intersection of VLo..VHi and RangeList 2742not_among_body([], _VLo, _VHi, [], _Flag). 2743not_among_body([H|T], VLo, VHi, List, Flag) :- 2744 ( integer(H) -> Lo = H, Hi = H ; H = Lo..Hi ), 2745 ( VHi < Lo -> % cannot be in this or later ranges 2746 List = [] 2747 ; VLo < Lo, VHi =< Hi -> % may be in this range, cannot be in later 2748 List = [H] 2749 ; VHi =< Hi -> % is in this range, set flag false 2750 Flag = 0 2751 ; % may be in this and/or later ranges 2752 ( VLo =< Hi -> List = [H|List1] ; List = List1 ), 2753 not_among_body(T, VLo, VHi, List1, Flag) 2754 ). 2755 2756:- mode cg_new_MP_columns(+,-). 2757cg_new_MP_columns(Handle, VarCols) :- 2758 Handle = cg_prob{sp_solution_call:SolveSubProblem, 2759 disallow:DisOnOff, 2760 mp_vars:OldVars, 2761 idx_lookup:Lookup, 2762 new_columns:NewColRec}, 2763 recorded_list(NewColRec, Solns), 2764 erase_all(NewColRec), 2765 arg(1, SolveSubProblem, SPHandle), 2766 SPHandle = sp_prob{coeff_vars:DualVars, 2767 disallow:[Count, Templates]}, 2768 ( 2769 foreach(Soln, Solns), 2770 fromto(VarCols, [Var:ObjCol|Rest], Rest, []), 2771 fromto(NewVars, [Var|Vars], Vars, OldVars), 2772 2773 % most general inst templates for disallowing cols 2774 % for thesis results tables 2775 fromto(Count, CountIn, CountOut, Count0), 2776 fromto(Templates, TemplatesIn, TemplatesOut, NewTemplates), 2777 param(DisOnOff), 2778 2779 param(Handle, DualVars, Lookup) 2780 do 2781 Soln = sp_sol{cost:Obj, 2782 coeff_vars:Coeffs, 2783 aux:Info, 2784 lo:Lo, 2785 hi:Hi, 2786 type:Type}, 2787 ( Obj =:= 0 -> ObjCol = BCol ; ObjCol = [obj:Obj|BCol] ), 2788 BCol = [lo:Lo, hi:Hi|Col], 2789 new_cg_attrstruct(Handle, Obj, HashCol, Lo, Hi, Type, Info, Attr), 2790 add_attribute(Var, Attr, colgen), 2791 ( Coeffs = [_Id-_V|_] -> 2792 % Coeffs is a sparse coefficient list (unordered) 2793 2794 % most general inst templates for disallowing cols 2795 % for thesis results tables 2796 ( DisOnOff == off -> 2797 CountOut = CountIn, TemplatesOut = TemplatesIn 2798 ; 2799 % reconstruct full coefficient list from sparse one 2800 ( 2801 foreach(Var, DualVars), %+ 2802 foreach(TVal, TCoeffs), %- 2803 param(Handle, Coeffs) 2804 do 2805 get_idx(Var, Id, Handle), 2806 ( member(Id-TVal, Coeffs) -> true 2807 ; TVal = 0 ) 2808 ), 2809 ( DisOnOff == lp -> 2810 lin_amalgamate(CountIn, TemplatesIn, 2811 TCoeffs, DualVars, 2812 CountOut, TemplatesOut) 2813 ; DisOnOff == clp -> 2814 amalgamate(CountIn, TemplatesIn, 2815 TCoeffs, DualVars, 2816 CountOut, TemplatesOut) 2817 ) 2818 ), 2819 2820 % make a sparse column representation Col1 2821 ( 2822 foreach(Id-V, Coeffs), %+ 2823 foreach(I:V, Col1), %- 2824 param(Lookup) 2825 do 2826 hash_get(Lookup, Id, I) 2827 ), 2828 keysort(Coeffs, HashCol) 2829 2830 ; 2831 % Coeffs is a full list of coefficients (not sparse) 2832 2833 % most general inst templates for disallowing cols 2834 % for thesis results tables 2835 ( DisOnOff == off -> 2836 CountOut = CountIn, TemplatesOut = TemplatesIn 2837 ; 2838 ( DisOnOff == lp -> 2839 lin_amalgamate(CountIn, TemplatesIn, 2840 Coeffs, DualVars, 2841 CountOut, TemplatesOut) 2842 ; DisOnOff == clp -> 2843 amalgamate(CountIn, TemplatesIn, 2844 Coeffs, DualVars, 2845 CountOut, TemplatesOut) 2846 ) 2847 ), 2848 2849 % make a sparse column representation Col1 2850 ( 2851 foreach(Var, DualVars), 2852 foreach(V, Coeffs), 2853 fromto(HC, HCOut, HCIn, []), 2854 fromto(Col1, Out, In, []), 2855 param(Handle, Lookup) 2856 do 2857 (V =:= 0 -> 2858 HCOut = HCIn, 2859 Out = In 2860 ; 2861 get_idx(Var, Id, Handle), 2862 HCOut = [Id-V|HCIn], 2863 hash_get(Lookup, Id, I), 2864 Out = [I:V|In] 2865 ) 2866 ), 2867 keysort(HC, HashCol) 2868 ), 2869 keysort(Col1, Col) 2870 ), 2871 2872 % most general inst templates for disallowing cols 2873 % for thesis results tables 2874 NewCount is max(Count, Count0), 2875 setarg(disallow of sp_prob,SPHandle, [NewCount, NewTemplates]), 2876 2877 setarg(mp_vars of cg_prob, Handle, NewVars). 2878 2879 2880:- demon cg_iteration/1. 2881:- set_flag(cg_iteration/1, priority, 7). 2882:- set_flag(cg_iteration/1, run_priority, 7). 2883%cg_iteration(+Handle) 2884% 2885% perform one column generation iteration: 2886% a) add any new constraints/variables and solve master problem, 2887% b) retrieve dual values and solve subproblems, 2888% c) retrieve beneficial columns and add to master problem 2889cg_iteration(Handle) :- 2890 % collect any new constraints or variables added since last 2891 % iteration 2892 %cg_info_message(Handle, "%tExtending restricted master problem ... ", []), 2893 add_new_solver_rowcols(Handle), 2894 % solve the new master problem 2895 cg_info_message(Handle, "%tSolving restricted master problem ... ", []), 2896 ( in_phase_1(Handle) -> 2897 solve_masterproblem_phase1(Handle, Status, ObjVal) 2898 ; 2899 solve_masterproblem_phase2(Handle, Status, ObjVal) 2900 ), 2901 % solve the subproblems 2902 solve_subproblems(Status, Handle), 2903 % retrieve beneficial columns 2904 cg_new_MP_columns(Handle, VarCols), 2905 ( stabilisation_stopping_criteria(Handle, VarCols) -> 2906 % cannot improve solution by adding new columns, 2907 % if we are in phase 1 (phase == 0) and we have hit 2908 % the stopping conditions we do not have a feasible 2909 % solution yet so we fail, while if we are in phase 2 2910 % (phase == -1) and we have hit the stopping 2911 % conditions then we have a feasible solution that we 2912 % cannot improve upon 2913 \+ in_phase_1(Handle) 2914 ; 2915 % update the master problem stabilisation terms 2916 update_stabilisation(Handle, VarCols), 2917 % update the degeneracy status 2918 update_degeneracy_status(Handle), 2919 % update the master problem lp 2920 update_masterproblem_lp(Handle, VarCols), 2921 % update the master problem lower bound 2922 update_lower_bound(Handle, ObjVal) 2923 ). 2924 2925in_phase_1(cg_prob{phase:0}). 2926 2927%solve_masterproblem_phase1(+Handle, -Status, -ObjVal) 2928% 2929% solve phase 1 of the two-phase method, i.e. with artifical cost 2930% function: we want to change to phase 2 iff the current 2931% problem has objective value 0. In order to avoid problems with 2932% rounding we do not solve and test solution value, but rather try to 2933% do the phase change and solve the phase 2 problem first. If that 2934% succeeds, then the phase 1 problem must have had objective value 0, 2935% so we return a phase2 Status with the value. If it fails, then the 2936% phase 1 problem must have had strictly positive objective value, so 2937% we try to solve the phase 1 problem. If that succeeds, we are still 2938% in phase 1 and return a phase 1 flag with the value. If both have 2939% failed then the phase 1 problem has become infeasible: since all 2940% constraints involving generated vars have an artificial variable 2941% they should always be satisfiable, thus any infeasibility is due to 2942% constraints involving no generated vars. If such a constraint is 2943% infeasible now then it will always be infeasible in future, so we 2944% fail. 2945% NB: for safety we really should check which constraints were 2946% unsatisfiable and abort if they contain a constraint involving 2947% generated vars 2948solve_masterproblem_phase1(Handle, Status, ObjVal) ?- 2949 % try phase change and phase 2 solution now 2950 Handle = cg_prob{master_prob:MP, 2951 phase1_obj:Phase1Fn, 2952 stabilisation:Stabilisation, 2953 stab_terms:StabTerms, 2954 mp_vars:Vars}, 2955 phase_change(Phase1Fn), 2956 lp_get(MP, objective, ObjExpr), 2957% ( lp_get(MP, cbasisarr, CB) -> true ; CB=[] ), 2958 solve_perturbed_objective_function(Stabilisation, StabTerms, 2959 MP, ObjExpr, ObjVal), !, 2960% ( lp_get(MP, cbasisarr, CB1) -> true ; CB1=[] ), 2961% write_term(CB, [compact(true)]), nl, 2962% write_term(CB1, [compact(true)]), nl, 2963 % phase change and phase 2 solution successful, 2964 % we have a feasible solution change the phase setting 2965 setarg(phase of cg_prob, Handle, -1), 2966 % any optimal solution of a phase 2 master problem is an upper 2967 % bound on the true optimal value, so update the upper bound ... 2968 setarg(upper_bound of cg_prob, Handle, ObjVal), 2969 % ... and set optimal vals 2970 set_optimal_mp_vals(Vars, Handle, MP), 2971 % no need to perturb the dual values since at the very least 2972 % the objective coefficient "dual" has changed from 0 to -1 2973 cg_info_message(Handle, "done, z_mp = %w%n", [ObjVal]), 2974 Status = phase2. 2975solve_masterproblem_phase1(Handle, Status, ObjVal) ?- 2976 % phase change and phase 2 solution was unsuccessful, 2977 % try phase 1 solution 2978 Handle = cg_prob{master_prob:MP, 2979 phase1_obj:Phase1Fn, 2980 on_degeneracy:OnDegeneracy, 2981 stabilisation:Stabilisation, 2982 stab_terms:StabTerms}, 2983 ( lp_get(MP, cbasis, StartCBasis) -> true ; StartCBasis = [] ), 2984 ( lp_get(MP, rbasis, StartRBasis) -> true ; StartRBasis = [] ), 2985% ( lp_get(MP, cbasisarr, CB) -> true ; CB=[] ), 2986 solve_perturbed_objective_function(Stabilisation, StabTerms, 2987 MP, Phase1Fn, ObjVal), !, 2988 % phase 1 solution successful, still restricted infeasible 2989 % if we have retained the same optimal basis try to perturb 2990 % the dual values so that the subproblem can generate 2991 % different columns 2992 perturb_if_necessary(Handle, Phase1Fn, StartCBasis, StartRBasis), 2993 cg_info_message(Handle, "infeasible%n", []), 2994% ( lp_get(MP, cbasisarr, CB1) -> true ; CB1=[] ), 2995% write_term(CB, [compact(true)]), nl, 2996% write_term(CB1, [compact(true)]), nl, 2997 ( lp_get(MP, cbasis, StartCBasis), 2998 lp_get(MP, rbasis, StartRBasis) -> 2999 % degenerate so Status is infeasible if not handled, 3000 % and degenerate if it is handled by the subproblem 3001 ( OnDegeneracy == stop -> Status = infeasible 3002 ; OnDegeneracy == continue -> Status = degenerate ) 3003 ; 3004 % not degenerate so Status is phase1 3005 Status = phase1 3006 ). 3007solve_masterproblem_phase1(Handle, _Status, _ObjVal) ?- 3008 % phase 1 solution unsuccessful, fail 3009 cg_info_message(Handle, "failed: inconsistent constraints posted%n", []), 3010 fail. 3011 3012%solve_masterproblem_phase2(+Handle, -Status, -ObjVal) 3013% 3014% solve phase 2 of the two-phase method, i.e. with true objective 3015% function and all artificial variables fixed to 0 3016solve_masterproblem_phase2(Handle, Status, ObjVal) ?- 3017 % ensure the artificial variables are ground just in case 3018 Handle = cg_prob{master_prob:MP, 3019 phase1_obj:Phase1Fn, 3020 on_degeneracy:OnDegeneracy, 3021 stabilisation:Stabilisation, 3022 stab_terms:StabTerms, 3023 mp_vars:Vars, 3024 tolerance:Tolerance, 3025 lower_bound:LowerBound}, 3026 ( lp_get(MP, cbasis, StartCBasis) -> true ; StartCBasis = [] ), 3027 ( lp_get(MP, rbasis, StartRBasis) -> true ; StartRBasis = [] ), 3028% ( lp_get(MP, cbasisarr, CB) -> true ; CB=[] ), 3029 % TODO: don't zero out when already done 3030 phase_change(Phase1Fn), 3031 % TODO: Don't retrieve and reconstruct potentially big objective 3032 % TODO: at least use a sum(List) instead of nested +/2 3033 lp_get(MP, objective, ObjExpr), 3034 solve_perturbed_objective_function(Stabilisation, StabTerms, 3035 MP, ObjExpr, ObjVal), 3036 % any optimal solution of a phase 2 master problem is an upper 3037 % bound on the true optimal value, so update the upper bound ... 3038 setarg(upper_bound of cg_prob, Handle, ObjVal), 3039 % ... and set optimal vals 3040 set_optimal_mp_vals(Vars, Handle, MP), 3041 % if we have retained the same optimal basis try to perturb 3042 % the dual values so that the subproblem can generate 3043 % different columns 3044 perturb_if_necessary(Handle, ObjExpr, StartCBasis, StartRBasis), 3045 cg_info_message(Handle, "done, z_mp = %w%n", [ObjVal]), 3046% ( lp_get(MP, cbasisarr, CB1) -> true ; CB1=[] ), 3047% write_term(CB, [compact(true)]), nl, 3048% write_term(CB1, [compact(true)]), nl, 3049 ( ObjVal - (Tolerance*ObjVal) =< LowerBound -> 3050 % within Tolerance of lower bound so status is Optimal 3051 Status = optimal 3052 ; lp_get(MP, cbasis, StartCBasis), 3053 lp_get(MP, rbasis, StartRBasis) -> 3054 % degenerate so Status is suboptimal if not handled, 3055 % and degenerate if it is handled by the subproblem 3056 ( OnDegeneracy == stop -> Status = suboptimal 3057 ; OnDegeneracy == continue -> Status = degenerate ) 3058 ; 3059 % not degenerate so Status is phase2 3060 Status = phase2 3061 ). 3062 3063phase_change(min(sum(ArtVars))) :- 3064 % zero out the artificial variables 3065 ( 3066 foreach(0, ArtVars) 3067 do 3068 true 3069 ). 3070 3071solve_perturbed_objective_function(off, StabTerms, 3072 MP, ObjExpr, ObjVal) ?- 3073 % no stabilisation policy, zero out the stabilisation vars ... 3074 %TODO: don't do that repeatedly 3075 ( 3076 foreach(StabTerm, StabTerms) 3077 do 3078 StabTerm = stab_term{plus_var:0.0,minus_var:0.0} 3079 ), 3080 % ... and solve for the existing objective function 3081 lp_probe(MP, ObjExpr, ObjVal). 3082% lp_get(MP, constraints, _), 3083% lp_probe(MP, ObjExpr, ObjVal), 3084% lp_get(MP, typed_solution, Sol), 3085% writeln(mp_sol:Sol). 3086solve_perturbed_objective_function(on(_,_,_,_), StabTerms, 3087 MP, ObjExpr, ObjVal) ?- 3088 solve_perturbed_objective_function_(StabTerms, MP, ObjExpr, ObjVal). 3089solve_perturbed_objective_function(stab_pred(_,_), StabTerms, 3090 MP, ObjExpr, ObjVal) ?- 3091 solve_perturbed_objective_function_(StabTerms, MP, ObjExpr, ObjVal). 3092 3093solve_perturbed_objective_function_(StabTerms, MP, ObjExpr, ObjVal) :- 3094 % some form of stabilisation policy (either default or 3095 % user-defined) is in use, perturb the objective function ... 3096 ( 3097 foreach(StabTerm, StabTerms), 3098 fromto(SExpr, Out, In, Expr) 3099 do 3100 StabTerm = stab_term{plus_var: Yplus, 3101 plus_coeff: CoeffPlus, 3102 minus_var: Yminus, 3103 minus_coeff: CoeffMinus}, 3104 Out = CoeffPlus*Yplus - CoeffMinus*Yminus + In 3105 ), 3106 ObjExpr =.. [Sense, Expr], 3107 StabExpr =.. [Sense, SExpr], 3108 % ... and solve for the perturbed objective function 3109 lp_probe(MP, StabExpr, ObjVal). 3110 3111set_optimal_mp_vals(Vars, Handle, MP) :- 3112 ( 3113 foreach(Var, Vars), 3114 param(Handle, MP) 3115 do 3116 ( nonvar(Var) -> true 3117 ; lp_var_get(MP, Var, solution, Sol), 3118 get_cg_attr(Var, Handle, Attr), 3119 setarg(mp_val of colgen, Attr, Sol) 3120 ) 3121 ). 3122 3123%solve_subproblems(++Status, +Handle) 3124% 3125% Status in {phase1,phase2,optimal,suboptimal,infeasible,degenerate} 3126% 3127% solve subproblems dealing with degeneracy appropriately for the 3128% current phase and parameters 3129solve_subproblems(infeasible, Handle) ?- !, 3130 % apparently degenerate and infeasible (in phase 1), 3131 % no degeneracy handling - 3132 % fail 3133 cg_info_message(Handle, "%t... detected identical" 3134 " external solver basis after mp" 3135 " optimisation%n%t terminating with" 3136 " no solution where one may exist%n", []), 3137 fail. 3138solve_subproblems(suboptimal, Handle) ?- !, 3139 % apparently degenerate and feasible (in phase 2), 3140 % no degeneracy handling - 3141 % terminate column generation 3142 cg_info_message(Handle, "%t... detected identical" 3143 " external solver basis after mp" 3144 " optimisation%n%t terminating with" 3145 " potentially suboptimal solution%n", []). 3146solve_subproblems(optimal, _Handle) ?- !. 3147 % within tolerance of optimal, no need to solve subproblems 3148solve_subproblems(Status, Handle) :- 3149 % either degenerate where the user-defined subproblem solution 3150 % predicate will try to deal with degeneracy, or not 3151 % degenerate - 3152 % flag status, get dual values and update dual_var attributes, 3153 % triggers SP solution 3154 cg_info_message(Handle, "%tSolving subproblems ... ", []), 3155 Handle = cg_prob{master_prob:MP, 3156 sp_solution_call:SolveSubProblem, 3157 idx_lookup:Lookup, 3158 shelf:Shelf, 3159 new_columns:NewColRec, 3160 phase:OptDual}, 3161 lp_get(MP, dual_solution, Duals), 3162 DualArr =.. [[]|Duals], 3163 setarg(duals of cg_prob, Handle, DualArr), 3164 arg(1, SolveSubProblem, SPHandle), 3165 setarg(status of sp_prob, SPHandle, Status), 3166 shelf_set(Shelf, optimal_rc of cg_shelf, none), 3167 verify(recorded_list(NewColRec,[])), 3168 call_priority(update_duals(Handle, SPHandle, OptDual, Lookup, DualArr), 3169 2), 3170 % Subproblem solver woken here!!! 3171 cg_info_message(Handle, "done%n", []). 3172 3173 3174stabilisation_stopping_criteria(Handle, VarCols) :- 3175 Handle = cg_prob{duals:Duals, 3176 module:Module, 3177 master_prob: MPHandle, 3178 sp_solution_call: SolveSubProblem, 3179 stabilisation: Stabilisation, 3180 stab_terms: StabTerms, 3181 tolerance: Tolerance, 3182 phase:OptDual}, 3183 ( Stabilisation == off -> 3184 VarCols = [] 3185 ; Stabilisation = stab_pred(_, Goal) -> 3186 call(Goal)@Module 3187 ; Stabilisation = on(_BoundIter, _BoundUpdate, 3188 _CoeffIter, _CoeffUpdate) -> 3189 ( 3190 foreach(_Var:ObjCol, VarCols), 3191 param(Duals, OptDual, Tolerance) 3192 do 3193 % check all cols have wrong reduced cost 3194 ( ObjCol = [obj:Obj|Col] -> true 3195 ; Obj = 0, ObjCol = Col ), 3196 Cost is OptDual*Obj, 3197 ( 3198 foreach(I:V, Col), 3199 fromto(Cost, In, Out, RC), 3200 param(Duals) 3201 do 3202 I1 is I + 1, 3203 arg(I1, Duals, Dual), 3204 Out is In + Dual*V 3205 ), 3206 RC =< -Tolerance 3207 ), 3208 ( 3209 foreach(StabTerm, StabTerms), 3210 param(MPHandle, Tolerance) 3211 do 3212 StabTerm = stab_term{plus_var: Yplus, 3213 plus_bound: Boundplus, 3214 minus_var: Yminus, 3215 minus_bound: Boundminus}, 3216 lp_var_get(MPHandle, Yplus, solution, YplusVal), 3217 abs(YplusVal) =< Tolerance, 3218 Boundplus =< Tolerance, 3219 lp_var_get(MPHandle, Yminus, solution, YminusVal), 3220 abs(YminusVal) =< Tolerance, 3221 Boundminus =< Tolerance 3222 ) 3223 ). 3224 3225update_stabilisation(Handle, VarCols) :- 3226 Handle = cg_prob{duals:Duals, 3227 module:Module, 3228 sp_solution_call: SolveSubProblem, 3229 master_prob: MPHandle, 3230 stabilisation: Stabilisation, 3231 stab_terms: StabTerms, 3232 stab_iter_counts: StabCounts, 3233 tolerance: Tolerance, 3234 phase:OptDual}, 3235 ( Stabilisation == off -> 3236 true 3237 ; Stabilisation = stab_pred(Goal, _) -> 3238 call(Goal)@Module 3239 ; Stabilisation = on(BoundIter, BoundUpdate, 3240 CoeffIter, CoeffUpdate) -> 3241 StabCounts = stab_counters{bound_counter:BoundCount, 3242 coeff_counter:CoeffCount}, 3243 ( BoundCount < BoundIter -> 3244 BoundCount1 is BoundCount + 1, 3245 setarg(bound_counter of stab_counters, StabCounts, 3246 BoundCount1) 3247 ; 3248 update_stab_bounds(StabTerms, MPHandle, BoundUpdate), 3249 setarg(bound_counter of stab_counters, StabCounts, 3250 1) 3251 ), 3252 ( CoeffCount < CoeffIter -> 3253 CoeffCount1 is CoeffCount + 1, 3254 setarg(coeff_counter of stab_counters, StabCounts, 3255 CoeffCount1) 3256 ; update_stab_coeffs(StabTerms, VarCols, OptDual, 3257 Duals, Tolerance, CoeffUpdate) -> 3258 setarg(coeff_counter of stab_counters, StabCounts, 3259 1) 3260 ; 3261 true 3262 ) 3263 ). 3264 3265update_stab_bounds(StabTerms, MPHandle, BoundUpdate) :- 3266 ( 3267 foreach(StabTerm, StabTerms), 3268 param(MPHandle, BoundUpdate) 3269 do 3270 StabTerm = stab_term{plus_var: Yplus, 3271 plus_bound: Boundplus, 3272 minus_var: Yminus, 3273 minus_bound: Boundminus}, 3274 Boundplus1 is max(0, Boundplus - BoundUpdate), 3275 lp_var_non_monotonic_set_bounds(MPHandle, Yplus, 0, Boundplus1), 3276 setarg(plus_bound of stab_term, StabTerm, Boundplus1), 3277 Boundminus1 is max(0, Boundminus - BoundUpdate), 3278 lp_var_non_monotonic_set_bounds(MPHandle, Yminus, 0, Boundminus1), 3279 setarg(minus_bound of stab_term, StabTerm, Boundminus1) 3280 ). 3281 3282update_stab_coeffs(StabTerms, VarCols, OptDual, 3283 Duals, Tolerance, 3284 CoeffUpdate) :- 3285 ( 3286 foreach(_Var:ObjCol, VarCols), 3287 param(OptDual, Duals, Tolerance) 3288 do 3289 % check all cols have wrong reduced cost 3290 ( ObjCol = [obj:Obj|Col] -> true 3291 ; Obj = 0, ObjCol = Col ), 3292 Cost is OptDual*Obj, 3293 ( 3294 foreach(I:V, Col), 3295 fromto(Cost, In, Out, RC), 3296 param(Duals) 3297 do 3298 I1 is I + 1, 3299 arg(I1, Duals, Dual), 3300 Out is In + Dual*V 3301 ), 3302 RC < Tolerance 3303 ), 3304 ( 3305 foreach(StabTerm, StabTerms), 3306 param(Duals, CoeffUpdate) 3307 do 3308 StabTerm = stab_term{idx:Idx}, 3309 Idx1 is Idx + 1, 3310 arg(Idx1, Duals, Dual), 3311 Coeffplus is Dual + CoeffUpdate, 3312 setarg(plus_coeff of stab_term, StabTerm, Coeffplus), 3313 Coeffminus is Dual - CoeffUpdate, 3314 setarg(minus_coeff of stab_term, StabTerm, Coeffminus) 3315 ). 3316 3317update_degeneracy_status(cg_prob{sp_solution_call:SolveSubProblem, phase:CD}) :- 3318 ( CD == 0 -> Status = phase1 ; Status = phase2 ), 3319 arg(1, SolveSubProblem, SPHandle), 3320 setarg(status of sp_prob, SPHandle, Status). 3321 3322update_masterproblem_lp(Handle, VarCols) :- 3323 Handle = cg_prob{master_prob:MP}, 3324 % add new columns to MP 3325 lp_add_columns(MP, VarCols), 3326 lp_get(MP, vars, MPVarArr), 3327 arity(MPVarArr, NewColsAdded), 3328 %%% WHY not reuse eplex's array here 3329 MPVarArr =.. [_|NewMPVars], 3330 setarg(mp_vars of cg_prob, Handle, NewMPVars), 3331 setarg(mp_cols_added of cg_prob, Handle, NewColsAdded). 3332 3333update_lower_bound(Handle, ObjVal) :- 3334 % calculate the Lasdon bound 3335 ( lasdon_bound(Handle, ObjVal, LasdonBound) -> 3336 Handle = cg_prob{tolerance:Tolerance, upper_bound:UpperBound}, 3337 ( LasdonBound >= UpperBound - (Tolerance * UpperBound) -> 3338 % cannot improve solution, done 3339 true 3340 ; 3341 % schedule the next MP iteration 3342 schedule_suspensions(mp_susp of cg_prob, Handle), 3343 wake 3344 ) 3345 ; 3346 % no bound available 3347 % 3348 % (note: if the subproblem solver is guaranteed to return 3349 % an optimal column we could just find the best reduced 3350 % cost as in update_stab_coeffs/6 and use that to 3351 % calculate the bound, but we allow solvers that return 3352 % potentially suboptimal columns; an improvement would be 3353 % to allow solvers to post an "optimality" flag along 3354 % with the column so that on any given iteration we could 3355 % update the bound iff the flag has been encountered for a 3356 % column; for multiple subproblems we would need to take 3357 % the maxium reduced cost of the "optimal" flagged 3358 % columns) 3359 % 3360 % schedule the next MP iteration 3361 schedule_suspensions(mp_susp of cg_prob, Handle), 3362 wake 3363 ). 3364 3365lasdon_bound(Handle, ObjVal, LasdonBound) :- 3366 Handle = cg_prob{lower_bound:LowerBound,shelf:Shelf}, 3367 \+ in_phase_1(Handle), 3368 shelf_get(Shelf, optimal_rc of cg_shelf, RCSum), 3369 number(RCSum), 3370 RCSum < 1.0Inf, 3371 LasdonBound is max(ObjVal - RCSum, LowerBound), 3372 setarg(lower_bound of cg_prob, Handle, LasdonBound), 3373 cg_info_message(Handle, "New lower bound: %w%n", [LasdonBound]). 3374 3375%perturb_if_necessary(+Handle, +ObjExpr, ++CBasis, ++RBasis) 3376% 3377% iff the new basis retrieved from the LP is the same as the [CR]Basis 3378% args, and basis_perturbation is on, attempt to perturb the lp 3379% solution, so that we avoid identical dual values and repeated 3380% columns returned from the subproblems, which would result in a loop. 3381perturb_if_necessary(cg_prob{basis_perturbation:off}, _ObjExpr, _CBasis, _RBasis) :- !. 3382perturb_if_necessary(Handle, ObjExpr, CBasis, RBasis) :- !, 3383 Handle = cg_prob{master_prob:MPHandle}, 3384 ( lp_get(MPHandle, cbasis, CBasis), 3385 lp_get(MPHandle, rbasis, RBasis) -> 3386 lp_get(optimizer, Optimizer), 3387 % the actual parameters to change are optimizer specific 3388 optimizer_specific_perturbation(Optimizer, ObjExpr, MPHandle), 3389 lp_get(MPHandle, cbasis, NewCBasis), 3390 lp_get(MPHandle, rbasis, NewRBasis), 3391 ( CBasis = NewCBasis, 3392 RBasis = NewRBasis -> writeln(perturbation-failed) ; true ) 3393 ; 3394 true 3395 ). 3396 3397optimizer_specific_perturbation(cplex, ObjExpr, MPHandle) ?- !, 3398 lp_get(optimizer_param(perind), Ind), 3399 lp_get(optimizer_param(perlim), Lim), 3400 lp_set(optimizer_param(perind), 1), 3401 lp_set(optimizer_param(perlim), 1), 3402 lp_probe(MPHandle, ObjExpr, _), 3403 lp_set(optimizer_param(perind), Ind), 3404 lp_set(optimizer_param(perlim), Lim). 3405optimizer_specific_perturbation(xpress, ObjExpr, MPHandle) ?- !, 3406 lp_get(optimizer_param(perturb), Val), 3407 lp_set(optimizer_param(perturb), 1.0), 3408 lp_probe(MPHandle, ObjExpr, _), 3409 lp_set(optimizer_param(perturb), Val). 3410%TODO: parameters for COIN 3411optimizer_specific_perturbation(Solver, _ObjExpr, _MPHandle) :- 3412 printf(error, "Colgen: %w not supported%n", [Solver]), 3413 abort. 3414 3415 3416update_duals(Handle, SPHandle, OptDual, Lookup, Duals) :- 3417 SPHandle = sp_prob{cost:OptVar, 3418 coeff_vars:DualVars, 3419 cutoff:Cutoff}, 3420 ( nonvar(OptVar) -> 3421 Cutoff1 is Cutoff + OptDual*OptVar, 3422 setarg(cutoff of sp_prob, SPHandle, Cutoff1) 3423 ; 3424 always_set_dual(OptVar, OptDual, Handle) 3425 ), 3426 ( 3427 foreach(DualVar, DualVars), 3428 param(Handle, Lookup, Duals) 3429 do 3430 get_idx(DualVar, Ident, Handle), 3431 hash_get(Lookup, Ident, Idx), 3432 Idx1 is Idx + 1, 3433 arg(Idx1, Duals, Dual), 3434 always_set_dual(DualVar, Dual, Handle) 3435 ). 3436/* 3437update_duals(SPHandle, Mu, Nu, OptDual, Duals, ObjVal, ZInc) :- 3438 SPHandle = sp_prob{cost:OptVar, 3439 coeff_vars:DualVars, 3440 gub_coeff_vars:Branches}, 3441 set_dual(OptVar, OptDual), 3442 3443 ( get_idx(Mu, MuIdx) -> 3444 % calculate Vanderbeck-Wolsey bound 3445 % Delta1 = Z - sum(PiB) - sum(MuK) - sum(NuL), 3446 % Delta2 = ZInc - sum(PiB) - sum(MuK) - sum(NuL), 3447 % Rho1 = max(Delta1/K0, Delta1/L0), 3448 % Rho2 = max(Delta2/K0, Delta2/L0), 3449 % Rho = min(Rho1, Rho2), 3450 % SPOpt > Mu0 + Nu0 - Rho 3451 % approximated by 3452 % SPOpt >= Mu0 + Nu0 - Rho + 1e-05 3453 ( 3454 foreach(DualVar, DualVars), 3455 fromto(0, In, Out, PiB), 3456 param(Duals) 3457 do 3458 get_idx(DualVar, Idx), 3459 get_rhs(DualVar, Rhs), 3460 Idx1 is Idx + 1, 3461 arg(Idx1, Duals, Dual), 3462 Out is In + Dual*Rhs, 3463 set_dual(DualVar, Dual) 3464 ), 3465 ( 3466 foreach(Branch, Branches), 3467 3468 % foreach(cg_gub{dual_var:GUBDualVar}, Branches), 3469 fromto(PiB, In, Out, ConstTerm), 3470 param(Duals) 3471 do 3472 3473 Branch = cg_gub{type:Type, 3474 dual_var:GUBDualVar}, 3475 3476 ( atom(Type) -> 3477 3478 get_rhs(GUBDualVar, Rhs), 3479 3480 ( (Rhs = 0, Type \= (>=)) -> 3481 Out = In 3482 ; 3483 3484 get_idx(GUBDualVar, Idx), 3485 Idx1 is Idx + 1, 3486 arg(Idx1, Duals, Dual), 3487 Out is In + Dual*Rhs, 3488 set_dual(GUBDualVar, Dual) 3489 3490 ) 3491 3492 ; 3493 ( 3494 foreach(T, Type), 3495 foreach(GUBDV, GUBDualVar), 3496 fromto(In, I, O, Out), 3497 param(Duals) 3498 do 3499 get_rhs(GUBDV, Rhs), 3500 ( (Rhs = 0, T \= (>=)) -> 3501 O = I, 3502 % the constraint was not explicitly 3503 % added to the lp because it just 3504 % fixes vars to 0, 3505 % no constraint, no dual val, 3506 % no contribution to ConstTerm 3507 % but we MUST still give it a big 3508 % negative dual val for the sps 3509 set_dual(GUBDV, -1000000000) 3510 ; 3511 get_idx(GUBDV, Idx), 3512 ( var(Idx) -> 3513 % the constraint was not 3514 % added to the lp because 3515 % it contained no vars yet 3516 % is it safe to set Dual = 0? 3517 Dual = 0 3518 ; 3519 3520 Idx1 is Idx + 1, 3521 arg(Idx1, Duals, Dual), 3522 3523 true 3524 ), 3525 3526 O is I + Dual*Rhs, 3527 set_dual(GUBDV, Dual) 3528 ) 3529 ) 3530 ) 3531 3532 ), 3533 % get_idx(Mu, MuIdx), 3534 get_rhs(Mu, K), 3535 MuIdx1 is MuIdx + 1, 3536 arg(MuIdx1, Duals, Mu0), 3537 set_dual(Mu, Mu0), 3538 get_idx(Nu, NuIdx), 3539 get_rhs(Nu, L), 3540 NuIdx1 is NuIdx + 1, 3541 arg(NuIdx1, Duals, Nu0), 3542 set_dual(Nu, Nu0), 3543 Delta1 is ObjVal - ConstTerm, 3544 Rho1 is max(Delta1/K, Delta1/L), 3545 Delta2 is ZInc - ConstTerm, 3546 Rho2 is max(Delta2/K, Delta2/L), 3547 ( Rho2 =< Rho1 -> 3548 % should prune node if no SP solutions 3549 Rho = Rho2 3550 ; 3551 % should terminate CG in node if no SP solutiuons 3552 Rho = Rho1 3553 ), 3554 % but since 3555 % SPOpt = max{(pi-e)X - fW + sum(muZ) + sum(nuZ) + mu0 + nu0} 3556 % we can just leave Mu0, Nu0 out? 3557 % CGCutoff is Mu0 + Nu0 - Rho + 1e-05, 3558 CGCutoff is 1e-05 - Rho, 3559 setarg(cutoff of sp_prob, SPHandle, CGCutoff) 3560 ; 3561 3562 ( 3563 foreach(DualVar, DualVars), 3564 param(Duals) 3565 do 3566 get_idx(DualVar, Idx), 3567 Idx1 is Idx + 1, 3568 arg(Idx1, Duals, Dual), 3569 set_dual(DualVar, Dual) 3570 ) 3571 3572 ). 3573*/ 3574 3575bp_separate(Handle) :- 3576 Handle = cg_prob{bfs_tree:BfsHandle, 3577 sep_call:Goal, 3578 shelf:Shelf}, 3579 Goal = call(SepGoal)@Module, 3580 ( SepGoal == true -> 3581 true 3582 ; 3583 \+ \+ ( bfs_impose_node_state(other, BfsHandle), 3584 call(Goal) 3585 ), 3586 shelf_get(Shelf, info of cg_shelf, Branches), 3587 shelf_set(Shelf, info of cg_shelf, []), 3588 ( 3589 foreach(Score:Branch, Branches), 3590 param(BfsHandle, Module) 3591 do 3592 bfs_branch(BfsHandle, [Score, call(Branch)@Module]) 3593 ) 3594 ). 3595 3596 3597fractional_vars(Vars, FracVars, Vals, Diffs, Fracs, L, U, Pool) :- 3598 Vars = [X|_], 3599 get_pool_handle(Handle, Pool), 3600 cg_var_get(Handle, X, mp_val, Sol), 3601 ( var(Sol) -> 3602 Handle = cg_prob{master_prob:MP}, 3603 ( 3604 foreach(Var, Vars), 3605 fromto(FracVars, VarsIn, VarsOut, []), 3606 fromto(Vals, ValsIn, ValsOut, []), 3607 fromto(Diffs, DiffsIn, DiffsOut, []), 3608 fromto(Fracs, FracsIn, FracsOut, []), 3609 fromto(0.0, LIn, LOut, L), 3610 fromto(1.0, UIn, UOut, U), 3611 param(MP) 3612 do 3613 ( nonvar(Var) -> 3614 VarsIn = VarsOut, 3615 ValsIn = ValsOut, 3616 DiffsIn = DiffsOut, 3617 FracsIn = FracsOut, 3618 LIn = LOut, 3619 UIn = UOut 3620 ; 3621 lp_var_get(MP, Var, solution, Val), 3622 Diff is abs(round(Val) - Val), 3623 ( Diff =< 1e-5 -> 3624 VarsIn = VarsOut, 3625 ValsIn = ValsOut, 3626 DiffsIn = DiffsOut, 3627 FracsIn = FracsOut, 3628 LIn = LOut, 3629 UIn = UOut 3630 ; 3631 Frac is Val - floor(Val), 3632 ( Frac < 0.5 -> 3633 LOut is max(Frac, LIn), 3634 UOut = UIn 3635 ; 3636 Frac > 0.5 -> 3637 LOut = LIn, 3638 UOut is min(Frac, UIn) 3639 ; 3640 LOut = 0.5, 3641 UOut = 0.5 3642 ), 3643 VarsIn = [Var|VarsOut], 3644 ValsIn = [Val|ValsOut], 3645 DiffsIn = [Diff|DiffsOut], 3646 FracsIn = [Frac|FracsOut] 3647 ) 3648 ) 3649 3650 ) 3651 ; 3652 ( 3653 foreach(Var, Vars), 3654 fromto(FracVars, VarsIn, VarsOut, []), 3655 fromto(Vals, ValsIn, ValsOut, []), 3656 fromto(Diffs, DiffsIn, DiffsOut, []), 3657 fromto(Fracs, FracsIn, FracsOut, []), 3658 fromto(0.0, LIn, LOut, L), 3659 fromto(1.0, UIn, UOut, U), 3660 param(Handle) 3661 do 3662 ( nonvar(Var) -> 3663 VarsIn = VarsOut, 3664 ValsIn = ValsOut, 3665 DiffsIn = DiffsOut, 3666 FracsIn = FracsOut, 3667 LIn = LOut, 3668 UIn = UOut 3669 ; 3670 cg_var_get(Handle, Var, mp_val, Val), 3671 Diff is abs(round(Val) - Val), 3672 ( Diff =< 1e-5 -> 3673 VarsIn = VarsOut, 3674 ValsIn = ValsOut, 3675 DiffsIn = DiffsOut, 3676 FracsIn = FracsOut, 3677 LIn = LOut, 3678 UIn = UOut 3679 ; 3680 Frac is Val - floor(Val), 3681 ( Frac < 0.5 -> 3682 LOut is max(Frac, LIn), 3683 UOut = UIn 3684 ; 3685 Frac > 0.5 -> 3686 LOut = LIn, 3687 UOut is min(Frac, UIn) 3688 ; 3689 LOut = 0.5, 3690 UOut = 0.5 3691 ), 3692 VarsIn = [Var|VarsOut], 3693 ValsIn = [Val|ValsOut], 3694 DiffsIn = [Diff|DiffsOut], 3695 FracsIn = [Frac|FracsOut] 3696 ) 3697 ) 3698 ) 3699 ). 3700 3701:- demon solveSPs/1. 3702:- set_flag(solveSPs/1, priority, 6). 3703:- set_flag(solveSPs/1, run_priority, 6). 3704solveSPs(Handle) :- 3705 ( 3706 Handle = cg_prob{sp_solution_call:SolveSubProblem, 3707 shelf:Shelf,module:Module}, 3708 arg(1, SolveSubProblem, SPHandle), 3709 % note SolveSubProblem MUST post at least one positive 3710 % reduced cost solution to the MP pool with 3711 % cg_subproblem_solution if any exist for any subproblem 3712 % otherwise we will terminate early with a suboptimal solution 3713 shelf_set(Shelf, sp_sol_cnt of cg_shelf, 0), 3714 3715 call(SolveSubProblem)@Module, 3716 3717 % If the SP has not posted solution(s) itself, post the current one 3718 ( shelf_get(Shelf, sp_sol_cnt of cg_shelf, 0) -> 3719 SPHandle = sp_prob{cost:Cost, coeff_vars:Coeffs, 3720 lo:Lo, hi:Hi, type:Type}, 3721 ( ground(Cost-Coeffs) -> 3722 cg_valid_columns(Handle, sp_sol{ 3723 cost:Cost, coeff_vars:Coeffs, lo:Lo, hi:Hi, type:Type}) 3724 ; 3725 writeln(warning_output, "Warning: subproblem solver " 3726 "succeeded with nonground variables (ignored)") 3727 ) 3728 ; 3729 shelf_set(Shelf, sp_sol_cnt of cg_shelf, 0) 3730 ), 3731 fail 3732 ; 3733 true 3734 ). 3735 3736 3737%----------------------------------------------------------------------- 3738% Pools 3739%----------------------------------------------------------------------- 3740 3741:- local record(colgen_pools). % list of colgen pool names 3742 3743create_colgen_pool(Pool) :- 3744 create_constraint_pool(Pool, property(arity) of cg_constraint_type, 3745 [ 3746 var_dual/6 -> var_dual1/7, 3747 get_dual/2 -> get_dual1/3, 3748 get_coeff/2 -> get_coeff1/3, 3749 get_idx/2 -> get_idx1/3, 3750 get_rhs/2 -> get_rhs1/3, 3751 always_set_dual/2 -> always_set_dual1/3, 3752 set_dual/2 -> set_dual1/3, 3753 solve/1 -> bp_solve1/2, 3754 solver_setup/3 -> cg_solver_setup/4, 3755 solver_setup/2 -> cg_solver_setup/3, 3756 var_get/3 -> cg_var_get1/4, 3757 get/2 -> cg_get1/3, 3758 set/2 -> cg_set/3, 3759 statistics/0 -> cg_statistics/1, 3760 integers/1 -> cg_integers1/2, 3761 identified_constraint/2 -> add_cg_pool_constraint/3, 3762 (::)/2 -> cg_range/3, 3763 (=:=)/2 -> cg_eq/3, 3764 (>=)/2 -> cg_ge/3, 3765 (=<)/2 -> cg_le/3, 3766 ($::)/2 -> cg_range/3, 3767 ($=)/2 -> cg_eq/3, 3768 ($>=)/2 -> cg_ge/3, 3769 ($=<)/2 -> cg_le/3, 3770 branch/1 -> cg_branch1/2, 3771 branch/2 -> cg_branch1/3, 3772 cg_subproblem_count/1 -> cg_sp_count1/2, 3773 cg_subproblem_solution/1 -> cg_valid_columns1/2, 3774 valid_columns/1 -> cg_valid_columns1/2, 3775 cg_subproblem_rc_sum/1 -> cg_sp_rc_sum/2, 3776 optimal_rc/1 -> cg_optimal_rc1/2, 3777 minimize/4 -> cg_minimize/5, 3778 minimize/3 -> cg_minimize/4 3779 ]). 3780 3781 3782colgen_instance(PoolName) :- 3783 ( is_constraint_pool(PoolName), 3784 recorded(colgen_pools, PoolName) % is a colgen pool 3785 -> 3786 % if pool exists, check if it is currently empty 3787 ( pool_is_empty(PoolName), 3788 get_pool_item(PoolName, 0) % has no associated solver 3789 -> 3790 true 3791 ; 3792 printf(error, "Colgen instance still in use in %w%n", [colgen_instance(PoolName)]), 3793 abort 3794 ) 3795 ; 3796% ( current_module(PoolName) -> 3797% error(6, colgen_instance(PoolName)) 3798% ; 3799 record(colgen_pools, PoolName), 3800 create_colgen_pool(PoolName) 3801% ) 3802 ). 3803 3804get_pool_handle(Handle, Pool) :- 3805 ( get_pool_item(Pool, Handle), Handle = cg_prob{} -> 3806 true 3807 ; 3808 get_pool_item(Pool, 0), 3809 init_cg_prob(Pool, Handle), 3810 set_pool_item(Pool, Handle) 3811 ). 3812 3813init_cg_prob(Pool, Handle) :- 3814 Handle = cg_prob{ 3815 pool:Pool, 3816 bfs_tree:[], 3817 new_columns:NewColRec, 3818 shelf:Shelf }, 3819 record_create(NewColRec), 3820 shelf_create(cg_shelf{ 3821 info:[], 3822 optimal_rc:none, 3823 sp_sol_cnt:0 3824 }, Shelf), 3825 init_suspension_list(mp_susp of cg_prob, Handle). 3826 3827 3828% ---------------------------------------------------------------------- 3829% Expression simplifier (using lib(linearize)) 3830% 3831% A linear expression is normalised into a list (sum) of the form 3832% [C0*1, C1*X1, C2*X2, ...] 3833% where Ci are numbers and Xi are distinct variables. 3834% The first (constant) term is always present, Ci (i>=1) are nonzero. 3835% The expression must be built from variables, numbers, +/2, */2, -/2, -/1. 3836% The simplifier fails if the expression is not linear. 3837% 3838% renormalise/2 renormalises a normal form expression after variable 3839% bindings, unifications. 3840% 3841% A normalised constraint has the form 3842% Sense:NormExpr 3843% where Sense is one of the atoms =:=, >= or =< and NormExpr is a 3844% normalised expression as above. E.g. (>=):[-5*1,3*X] encodes 3845% the constraint -5 + 3*X >= 0. 3846% ---------------------------------------------------------------------- 3847 3848cg_normalise_cstr(E1 =:= E2, Norm, CoeffVar, Coeff) :- !, 3849 cg_normalise_cstr(E1 $= E2, Norm, CoeffVar, Coeff). 3850cg_normalise_cstr(E1 >= E2, Norm, CoeffVar, Coeff) :- !, 3851 cg_normalise_cstr(E1 $>= E2, Norm, CoeffVar, Coeff). 3852cg_normalise_cstr(E1 =< E2, Norm, CoeffVar, Coeff) :- !, 3853 cg_normalise_cstr(E1 $=< E2, Norm, CoeffVar, Coeff). 3854cg_normalise_cstr(E1 $= E2, (=:=):Norm, CoeffVar, Coeff) :- !, 3855 linearize(E1-E2, Lin, NonLin), 3856 ( NonLin == [] -> 3857 Norm = Lin, 3858 CoeffVar = 0, 3859 Coeff = 0 3860 ; 3861 NonLin = [AuxVar = implicit_sum(CoeffVar)], 3862 filter_auxvar(AuxVar, Lin, Norm, Coeff) 3863 ). 3864cg_normalise_cstr(E1 $>= E2, (>=):Norm, CoeffVar, Coeff) :- !, 3865 linearize(E1-E2, Lin, NonLin), 3866 ( NonLin == [] -> 3867 Norm = Lin, 3868 CoeffVar = 0, 3869 Coeff = 0 3870 ; 3871 NonLin = [AuxVar = implicit_sum(CoeffVar)], 3872 filter_auxvar(AuxVar, Lin, Norm, Coeff) 3873 ). 3874cg_normalise_cstr(E1 $=< E2, (=<):Norm, CoeffVar, Coeff) :- !, 3875 linearize(E1-E2, Lin, NonLin), 3876 ( NonLin == [] -> 3877 Norm = Lin, 3878 CoeffVar = 0, 3879 Coeff = 0 3880 ; 3881 NonLin = [AuxVar = implicit_sum(CoeffVar)], 3882 filter_auxvar(AuxVar, Lin, Norm, Coeff) 3883 ). 3884cg_normalise_cstr(Cstr, _, _, _) :- 3885 writeln(error, "Error: Unknown constraint":Cstr), 3886 abort. 3887 3888 filter_auxvar(Var, [C*X|Terms], Norm, Coeff) :- 3889 ( Var == X -> 3890 Norm = Terms, 3891 Coeff = C 3892 ; 3893 Norm = [C*X|Norm0], 3894 filter_auxvar(Var, Terms, Norm0, Coeff) 3895 ). 3896