1/* Part of SWI-Prolog 2 3 Author: Markus Triska 4 E-mail: triska@metalevel.at 5 WWW: http://www.swi-prolog.org 6 Copyright (C): 2005-2016, Markus Triska 7 All rights reserved. 8 9 Redistribution and use in source and binary forms, with or without 10 modification, are permitted provided that the following conditions 11 are met: 12 13 1. Redistributions of source code must retain the above copyright 14 notice, this list of conditions and the following disclaimer. 15 16 2. Redistributions in binary form must reproduce the above copyright 17 notice, this list of conditions and the following disclaimer in 18 the documentation and/or other materials provided with the 19 distribution. 20 21 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 22 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 23 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 24 FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 25 COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 26 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 27 BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 28 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 29 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 30 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 31 ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 32 POSSIBILITY OF SUCH DAMAGE. 33*/ 34 35 36:- module(simplex, 37 [ 38 assignment/2, 39 constraint/3, 40 constraint/4, 41 constraint_add/4, 42 gen_state/1, 43 gen_state_clpr/1, 44 gen_state_clpr/2, 45 maximize/3, 46 minimize/3, 47 objective/2, 48 shadow_price/3, 49 transportation/4, 50 variable_value/3 51 ]). 52 53:- if(exists_source(library(clpr))). 54:- use_module(library(clpr)). 55:- endif. 56:- use_module(library(assoc)). 57:- use_module(library(pio)). 58 59:- autoload(library(apply), [foldl/4, maplist/2, maplist/3]). 60:- autoload(library(lists), [flatten/2, member/2, nth0/3, reverse/2, select/3]).

Solve linear programming problems

Introduction

A linear programming problem or simply linear program (LP) consists of:

a set of linear constraints
a set of variables
a linear objective function.

The goal is to assign values to the variables so as to maximize (or minimize) the value of the objective function while satisfying all constraints.

Many optimization problems can be modeled in this way. As one basic example, consider a knapsack with fixed capacity C, and a number of items with sizes s(i) and values v(i). The goal is to put as many items as possible in the knapsack (not exceeding its capacity) while maximizing the sum of their values.

As another example, suppose you are given a set of coins with certain values, and you are to find the minimum number of coins such that their values sum up to a fixed amount. Instances of these problems are solved below.

All numeric quantities are converted to rationals via rationalize/1, and rational arithmetic is used throughout solving linear programs. In the current implementation, all variables are implicitly constrained to be non-negative. This may change in future versions, and non-negativity constraints should therefore be stated explicitly.

Example 1

This is the "radiation therapy" example, taken from Introduction to Operations Research by Hillier and Lieberman.

Prolog DCG notation is used to implicitly thread the state through posting the constraints:

:- use_module(library(simplex)).

radiation(S) :-
        gen_state(S0),
        post_constraints(S0, S1),
        minimize([0.4*x1, 0.5*x2], S1, S).

post_constraints -->
        constraint([0.3*x1, 0.1*x2] =< 2.7),
        constraint([0.5*x1, 0.5*x2] = 6),
        constraint([0.6*x1, 0.4*x2] >= 6),
        constraint([x1] >= 0),
        constraint([x2] >= 0).

An example query:

?- radiation(S), variable_value(S, x1, Val1),
                 variable_value(S, x2, Val2).
Val1 = 15 rdiv 2,
Val2 = 9 rdiv 2.

Example 2

Here is an instance of the knapsack problem described above, where C = 8, and we have two types of items: One item with value 7 and size 6, and 2 items each having size 4 and value 4. We introduce two variables, x(1) and x(2) that denote how many items to take of each type.

:- use_module(library(simplex)).

knapsack(S) :-
        knapsack_constraints(S0),
        maximize([7*x(1), 4*x(2)], S0, S).

knapsack_constraints(S) :-
        gen_state(S0),
        constraint([6*x(1), 4*x(2)] =< 8, S0, S1),
        constraint([x(1)] =< 1, S1, S2),
        constraint([x(2)] =< 2, S2, S).

An example query yields:

?- knapsack(S), variable_value(S, x(1), X1),
                variable_value(S, x(2), X2).
X1 = 1
X2 = 1 rdiv 2.

That is, we are to take the one item of the first type, and half of one of the items of the other type to maximize the total value of items in the knapsack.

If items can not be split, integrality constraints have to be imposed:

knapsack_integral(S) :-
        knapsack_constraints(S0),
        constraint(integral(x(1)), S0, S1),
        constraint(integral(x(2)), S1, S2),
        maximize([7*x(1), 4*x(2)], S2, S).

Now the result is different:

?- knapsack_integral(S), variable_value(S, x(1), X1),
                         variable_value(S, x(2), X2).

X1 = 0
X2 = 2

That is, we are to take only the two items of the second type. Notice in particular that always choosing the remaining item with best performance (ratio of value to size) that still fits in the knapsack does not necessarily yield an optimal solution in the presence of integrality constraints.

Example 3

We are given:

3 coins each worth 1 unit
20 coins each worth 5 units and
10 coins each worth 20 units.

The task is to find a minimal number of these coins that amount to 111 units in total. We introduce variables c(1), c(5) and c(20) denoting how many coins to take of the respective type:

:- use_module(library(simplex)).

coins(S) :-
        gen_state(S0),
        coins(S0, S).

coins -->
        constraint([c(1), 5*c(5), 20*c(20)] = 111),
        constraint([c(1)] =< 3),
        constraint([c(5)] =< 20),
        constraint([c(20)] =< 10),
        constraint([c(1)] >= 0),
        constraint([c(5)] >= 0),
        constraint([c(20)] >= 0),
        constraint(integral(c(1))),
        constraint(integral(c(5))),
        constraint(integral(c(20))),
        minimize([c(1), c(5), c(20)]).

An example query:

?- coins(S), variable_value(S, c(1), C1),
             variable_value(S, c(5), C5),
             variable_value(S, c(20), C20).

C1 = 1,
C5 = 2,
C20 = 5.

author: - Markus Triska */

238/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 239 CLP(R) bindings 240 the (unsolved) state is stored as a structure of the form 241 clpr_state(Options, Cs, Is) 242 Options: list of Option=Value pairs, currently only eps=Eps 243 Cs: list of constraints, i.e., structures of the form 244 c(Name, Left, Op, Right) 245 anonymous constraints have Name == 0 246 Is: list of integral variables 247- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 248 249gen_state_clpr(State) :- gen_state_clpr([], State). 250 251gen_state_clpr(Options, State) :- 252 ( memberchk(eps=_, Options) -> Options1 = Options 253 ; Options1 = [eps=1.0e-6|Options] 254 ), 255 State = clpr_state(Options1, [], []). 256 257clpr_state_options(clpr_state(Os, _, _), Os). 258clpr_state_constraints(clpr_state(_, Cs, _), Cs). 259clpr_state_integrals(clpr_state(_, _, Is), Is). 260clpr_state_add_integral(I, clpr_state(Os, Cs, Is), clpr_state(Os, Cs, [I|Is])). 261clpr_state_add_constraint(C, clpr_state(Os, Cs, Is), clpr_state(Os, [C|Cs], Is)). 262clpr_state_set_constraints(Cs, clpr_state(Os,_,Is), clpr_state(Os, Cs, Is)). 263 264clpr_constraint(Name, Constraint, S0, S) :- 265 ( Constraint = integral(Var) -> clpr_state_add_integral(Var, S0, S) 266 ; Constraint =.. [Op,Left,Right], 267 coeff_one(Left, Left1), 268 clpr_state_add_constraint(c(Name, Left1, Op, Right), S0, S) 269 ). 270 271clpr_constraint(Constraint, S0, S) :- 272 clpr_constraint(0, Constraint, S0, S). 273 274clpr_shadow_price(clpr_solved(_,_,Duals,_), Name, Value) :- 275 memberchk(Name-Value0, Duals), 276 Value is Value0. 277 %( var(Value0) -> 278 % Value = 0 279 %; 280 % Value is Value0 281 %). 282 283 284 285clpr_make_variables(Cs, Aliases) :- 286 clpr_constraints_variables(Cs, Variables0, []), 287 sort(Variables0, Variables1), 288 pairs_keys(Aliases, Variables1). 289 290clpr_constraints_variables([]) --> []. 291clpr_constraints_variables([c(_, Left, _, _)|Cs]) --> 292 variables(Left), 293 clpr_constraints_variables(Cs). 294 295clpr_set_up(Aliases, c(_Name, Left, Op, Right)) :- 296 clpr_translate_linsum(Left, Aliases, LinSum), 297 CLPRConstraint =.. [Op, LinSum, Right], 298 clpr:{ CLPRConstraint }. 299 300clpr_set_up_noneg(Aliases, Var) :- 301 memberchk(Var-CLPVar, Aliases), 302 { CLPVar >= 0 }. 303 304clpr_translate_linsum([], _, 0). 305clpr_translate_linsum([Coeff*Var|Ls], Aliases, LinSum) :- 306 memberchk(Var-CLPVar, Aliases), 307 LinSum = Coeff*CLPVar + LinRest, 308 clpr_translate_linsum(Ls, Aliases, LinRest). 309 310clpr_dual(Objective0, S0, DualValues) :- 311 clpr_state_constraints(S0, Cs0), 312 clpr_constraints_variables(Cs0, Variables0, []), 313 sort(Variables0, Variables1), 314 maplist(clpr_standard_form, Cs0, Cs1), 315 clpr_include_all_vars(Cs1, Variables1, Cs2), 316 clpr_merge_into(Variables1, Objective0, Objective, []), 317 clpr_unique_names(Cs2, 0, Names), 318 maplist(clpr_constraint_coefficient, Cs2, Coefficients), 319 lists_transpose(Coefficients, TCs), 320 maplist(clpr_dual_constraints(Names), TCs, Objective, DualConstraints), 321 phrase(clpr_nonneg_constraints(Cs2, Names), DualNonNeg), 322 append(DualConstraints, DualNonNeg, DualConstraints1), 323 maplist(clpr_dual_objective, Cs2, Names, DualObjective), 324 clpr_make_variables(DualConstraints1, Aliases), 325 maplist(clpr_set_up(Aliases), DualConstraints1), 326 clpr_translate_linsum(DualObjective, Aliases, LinExpr), 327 minimize(LinExpr), 328 Aliases = DualValues. 329 330 331 332clpr_dual_objective(c(_, _, _, Right), Name, Right*Name). 333 334clpr_nonneg_constraints([], _) --> []. 335clpr_nonneg_constraints([C|Cs], [Name|Names]) --> 336 { C = c(_, _, Op, _) }, 337 ( { Op == (=<) } -> [c(0, [1*Name], (>=), 0)] 338 ; [] 339 ), 340 clpr_nonneg_constraints(Cs, Names). 341 342 343clpr_dual_constraints(Names, Coeffs, O*_, Constraint) :- 344 maplist(clpr_dual_linsum, Coeffs, Names, Linsum), 345 Constraint = c(0, Linsum, (>=), O). 346 347clpr_dual_linsum(Coeff, Name, Coeff*Name). 348 349clpr_constraint_coefficient(c(_, Left, _, _), Coeff) :- 350 maplist(all_coeffs, Left, Coeff). 351 352all_coeffs(Coeff*_, Coeff). 353 354clpr_unique_names([], _, []). 355clpr_unique_names([C0|Cs0], Num, [N|Ns]) :- 356 C0 = c(Name, _, _, _), 357 ( Name == 0 -> N = Num, Num1 is Num + 1 358 ; N = Name, Num1 = Num 359 ), 360 clpr_unique_names(Cs0, Num1, Ns). 361 362clpr_include_all_vars([], _, []). 363clpr_include_all_vars([C0|Cs0], Variables, [C|Cs]) :- 364 C0 = c(Name, Left0, Op, Right), 365 clpr_merge_into(Variables, Left0, Left, []), 366 C = c(Name, Left, Op, Right), 367 clpr_include_all_vars(Cs0, Variables, Cs). 368 369clpr_merge_into([], _, Ls, Ls). 370clpr_merge_into([V|Vs], Left, Ls0, Ls) :- 371 ( member(Coeff*V, Left) -> 372 Ls0 = [Coeff*V|Rest] 373 ; 374 Ls0 = [0*V|Rest] 375 ), 376 clpr_merge_into(Vs, Left, Rest, Ls). 377 378 379 380 381clpr_maximize(Expr0, S0, S) :- 382 coeff_one(Expr0, Expr), 383 clpr_state_constraints(S0, Cs), 384 clpr_make_variables(Cs, Aliases), 385 maplist(clpr_set_up(Aliases), Cs), 386 clpr_constraints_variables(Cs, Variables0, []), 387 sort(Variables0, Variables1), 388 maplist(clpr_set_up_noneg(Aliases), Variables1), 389 clpr_translate_linsum(Expr, Aliases, LinExpr), 390 clpr_state_integrals(S0, Is), 391 ( Is == [] -> 392 maximize(LinExpr), 393 Sup is LinExpr, 394 clpr_dual(Expr, S0, DualValues), 395 S = clpr_solved(Sup, Aliases, DualValues, S0) 396 ; 397 clpr_state_options(S0, Options), 398 memberchk(eps=Eps, Options), 399 maplist(clpr_fetch_var(Aliases), Is, Vars), 400 bb_inf(Vars, -LinExpr, Sup, Vertex, Eps), 401 pairs_keys_values(Values, Is, Vertex), 402 % what about the dual in MIPs? 403 Sup1 is -Sup, 404 S = clpr_solved(Sup1, Values, [], S0) 405 ). 406 407clpr_minimize(Expr0, S0, S) :- 408 coeff_one(Expr0, Expr1), 409 maplist(linsum_negate, Expr1, Expr2), 410 clpr_maximize(Expr2, S0, S1), 411 S1 = clpr_solved(Sup, Values, Duals, S0), 412 Inf is -Sup, 413 S = clpr_solved(Inf, Values, Duals, S0). 414 415clpr_fetch_var(Aliases, Var, X) :- memberchk(Var-X, Aliases). 416 417clpr_variable_value(clpr_solved(_, Aliases, _, _), Variable, Value) :- 418 memberchk(Variable-Value0, Aliases), 419 Value is Value0. 420 %( var(Value0) -> 421 % Value = 0 422 %; 423 % Value is Value0 424 %). 425 426clpr_objective(clpr_solved(Obj, _, _, _), Obj). 427 428clpr_standard_form(c(Name, Left, Op, Right), S) :- 429 clpr_standard_form_(Op, Name, Left, Right, S). 430 431clpr_standard_form_((=), Name, Left, Right, c(Name, Left, (=), Right)). 432clpr_standard_form_((>=), Name, Left, Right, c(Name, Left1, (=<), Right1)) :- 433 Right1 is -Right, 434 maplist(linsum_negate, Left, Left1). 435clpr_standard_form_((=<), Name, Left, Right, c(Name, Left, (=<), Right)). 436 437 438/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 439 General Simplex Algorithm 440 Structures used: 441 442 tableau(Objective, Variables, Indicators, Constraints) 443 *) objective function, represented as row 444 *) list of variables corresponding to columns 445 *) indicators denoting which variables are still active 446 *) constraints as rows 447 448 row(Var, Left, Right) 449 *) the basic variable corresponding to this row 450 *) coefficients of the left-hand side of the constraint 451 *) right-hand side of the constraint 452- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 453 454 455find_row(Variable, [Row|Rows], R) :- 456 Row = row(V, _, _), 457 ( V == Variable -> R = Row 458 ; find_row(Variable, Rows, R) 459 ).

466variable_value(State, Variable, Value) :- 467 functor(State, F, _), 468 ( F == solved -> 469 solved_tableau(State, Tableau), 470 tableau_rows(Tableau, Rows), 471 ( find_row(Variable, Rows, Row) -> Row = row(_, _, Value) 472 ; Value = 0 473 ) 474 ; F == clpr_solved -> clpr_variable_value(State, Variable, Value) 475 ). 476 477var_zero(State, _Coeff*Var) :- variable_value(State, Var, 0). 478 479list_first(Ls, F, Index) :- once(nth0(Index, Ls, F)).

487shadow_price(State, Name, Value) :- 488 functor(State, F, _), 489 ( F == solved -> 490 solved_tableau(State, Tableau), 491 tableau_objective(Tableau, row(_,Left,_)), 492 tableau_variables(Tableau, Variables), 493 solved_names(State, Names), 494 memberchk(user(Name)-Var, Names), 495 list_first(Variables, Var, Nth0), 496 nth0(Nth0, Left, Value) 497 ; F == clpr_solved -> clpr_shadow_price(State, Name, Value) 498 ).

505objective(State, Obj) :- 506 functor(State, F, _), 507 ( F == solved -> 508 solved_tableau(State, Tableau), 509 tableau_objective(Tableau, Objective), 510 Objective = row(_, _, Obj) 511 ; clpr_objective(State, Obj) 512 ). 513 514 % interface functions that access tableau components 515 516tableau_objective(tableau(Obj, _, _, _), Obj). 517tableau_rows(tableau(_, _, _, Rows), Rows). 518tableau_indicators(tableau(_, _, Inds, _), Inds). 519tableau_variables(tableau(_, Vars, _, _), Vars). 520 521/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 522 interface functions that access and modify state components 523 524 state is a structure of the form 525 state(Num, Names, Cs, Is) 526 Num: used to obtain new unique names for slack variables in a side-effect 527 free way (increased by one and threaded through) 528 Names: list of Name-Var, correspondence between constraint-names and 529 names of slack/artificial variables to obtain shadow prices later 530 Cs: list of constraints 531 Is: list of integer variables 532 533 constraints are initially represented as c(Name, Left, Op, Right), 534 and after normalizing as c(Var, Left, Right). Name of unnamed constraints 535 is 0. The distinction is important for merging constraints (mainly in 536 branch and bound) with existing ones. 537- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 538 539 540constraint_name(c(Name, _, _, _), Name). 541constraint_op(c(_, _, Op, _), Op). 542constraint_left(c(_, Left, _, _), Left). 543constraint_right(c(_, _, _, Right), Right).

549gen_state(state(0,[],[],[])). 550 551state_add_constraint(C, S0, S) :- 552 ( constraint_name(C, 0), constraint_left(C, [_Coeff*_Var]) -> 553 state_merge_constraint(C, S0, S) 554 ; state_add_constraint_(C, S0, S) 555 ). 556 557state_add_constraint_(C, state(VID,Ns,Cs,Is), state(VID,Ns,[C|Cs],Is)). 558 559state_merge_constraint(C, S0, S) :- 560 constraint_left(C, [Coeff0*Var0]), 561 constraint_right(C, Right0), 562 constraint_op(C, Op), 563 ( Coeff0 =:= 0 -> 564 ( Op == (=) -> Right0 =:= 0, S0 = S 565 ; Op == (=<) -> S0 = S 566 ; Op == (>=) -> Right0 =:= 0, S0 = S 567 ) 568 ; Coeff0 < 0 -> state_add_constraint_(C, S0, S) 569 ; Right is Right0 rdiv Coeff0, 570 state_constraints(S0, Cs), 571 ( select(c(0, [1*Var0], Op, CRight), Cs, RestCs) -> 572 ( Op == (=) -> CRight =:= Right, S0 = S 573 ; Op == (=<) -> 574 NewRight is min(Right, CRight), 575 NewCs = [c(0, [1*Var0], Op, NewRight)|RestCs], 576 state_set_constraints(NewCs, S0, S) 577 ; Op == (>=) -> 578 NewRight is max(Right, CRight), 579 NewCs = [c(0, [1*Var0], Op, NewRight)|RestCs], 580 state_set_constraints(NewCs, S0, S) 581 ) 582 ; state_add_constraint_(c(0, [1*Var0], Op, Right), S0, S) 583 ) 584 ). 585 586 587state_add_name(Name, Var), [state(VID,[Name-Var|Ns],Cs,Is)] --> 588 [state(VID,Ns,Cs,Is)]. 589 590state_add_integral(Var, state(VID,Ns,Cs,Is), state(VID,Ns,Cs,[Var|Is])). 591 592state_constraints(state(_, _, Cs, _), Cs). 593state_names(state(_,Names,_,_), Names). 594state_integrals(state(_,_,_,Is), Is). 595state_set_constraints(Cs, state(VID,Ns,_,Is), state(VID,Ns,Cs,Is)). 596state_set_integrals(Is, state(VID,Ns,Cs,_), state(VID,Ns,Cs,Is)). 597 598 599state_next_var(VarID0), [state(VarID1,Names,Cs,Is)] --> 600 [state(VarID0,Names,Cs,Is)], 601 { VarID1 is VarID0 + 1 }. 602 603solved_tableau(solved(Tableau, _, _), Tableau). 604solved_names(solved(_, Names,_), Names). 605solved_integrals(solved(_,_,Is), Is). 606 607% User-named constraints are wrapped with user/1 to also allow "0" in 608% constraint names.

623constraint(C, S0, S) :- 624 functor(S0, F, _), 625 ( F == state -> 626 ( C = integral(Var) -> state_add_integral(Var, S0, S) 627 ; constraint_(0, C, S0, S) 628 ) 629 ; F == clpr_state -> clpr_constraint(C, S0, S) 630 ).

637constraint(Name, C, S0, S) :- constraint_(user(Name), C, S0, S). 638 639constraint_(Name, C, S0, S) :- 640 functor(S0, F, _), 641 ( F == state -> 642 ( C = integral(Var) -> state_add_integral(Var, S0, S) 643 ; C =.. [Op, Left0, Right0], 644 coeff_one(Left0, Left), 645 Right0 >= 0, 646 Right is rationalize(Right0), 647 state_add_constraint(c(Name, Left, Op, Right), S0, S) 648 ) 649 ; F == clpr_state -> clpr_constraint(Name, C, S0, S) 650 ).

658constraint_add(Name, A, S0, S) :- 659 functor(S0, F, _), 660 ( F == state -> 661 state_constraints(S0, Cs), 662 add_left(Cs, user(Name), A, Cs1), 663 state_set_constraints(Cs1, S0, S) 664 ; F == clpr_state -> 665 clpr_state_constraints(S0, Cs), 666 add_left(Cs, Name, A, Cs1), 667 clpr_state_set_constraints(Cs1, S0, S) 668 ). 669 670 671add_left([c(Name,Left0,Op,Right)|Cs], V, A, [c(Name,Left,Op,Right)|Rest]) :- 672 ( Name == V -> append(A, Left0, Left), Rest = Cs 673 ; Left0 = Left, add_left(Cs, V, A, Rest) 674 ). 675 676branching_variable(State, Variable) :- 677 solved_integrals(State, Integrals), 678 member(Variable, Integrals), 679 variable_value(State, Variable, Value), 680 \+ integer(Value). 681 682 683worth_investigating(ZStar0, _, _) :- var(ZStar0). 684worth_investigating(ZStar0, AllInt, Z) :- 685 nonvar(ZStar0), 686 ( AllInt =:= 1 -> Z1 is floor(Z) 687 ; Z1 = Z 688 ), 689 Z1 > ZStar0. 690 691 692branch_and_bound(Objective, Solved, AllInt, ZStar0, ZStar, S0, S, Found) :- 693 objective(Solved, Z), 694 ( worth_investigating(ZStar0, AllInt, Z) -> 695 ( branching_variable(Solved, BrVar) -> 696 variable_value(Solved, BrVar, Value), 697 Value1 is floor(Value), 698 Value2 is Value1 + 1, 699 constraint([BrVar] =< Value1, S0, SubProb1), 700 ( maximize_(Objective, SubProb1, SubSolved1) -> 701 Sub1Feasible = 1, 702 objective(SubSolved1, Obj1) 703 ; Sub1Feasible = 0 704 ), 705 constraint([BrVar] >= Value2, S0, SubProb2), 706 ( maximize_(Objective, SubProb2, SubSolved2) -> 707 Sub2Feasible = 1, 708 objective(SubSolved2, Obj2) 709 ; Sub2Feasible = 0 710 ), 711 ( Sub1Feasible =:= 1, Sub2Feasible =:= 1 -> 712 ( Obj1 >= Obj2 -> 713 First = SubProb1, 714 Second = SubProb2, 715 FirstSolved = SubSolved1, 716 SecondSolved = SubSolved2 717 ; First = SubProb2, 718 Second = SubProb1, 719 FirstSolved = SubSolved2, 720 SecondSolved = SubSolved1 721 ), 722 branch_and_bound(Objective, FirstSolved, AllInt, ZStar0, ZStar1, First, Solved1, Found1), 723 branch_and_bound(Objective, SecondSolved, AllInt, ZStar1, ZStar2, Second, Solved2, Found2) 724 ; Sub1Feasible =:= 1 -> 725 branch_and_bound(Objective, SubSolved1, AllInt, ZStar0, ZStar1, SubProb1, Solved1, Found1), 726 Found2 = 0 727 ; Sub2Feasible =:= 1 -> 728 Found1 = 0, 729 branch_and_bound(Objective, SubSolved2, AllInt, ZStar0, ZStar2, SubProb2, Solved2, Found2) 730 ; Found1 = 0, Found2 = 0 731 ), 732 ( Found1 =:= 1, Found2 =:= 1 -> S = Solved2, ZStar = ZStar2 733 ; Found1 =:= 1 -> S = Solved1, ZStar = ZStar1 734 ; Found2 =:= 1 -> S = Solved2, ZStar = ZStar2 735 ; S = S0, ZStar = ZStar0 736 ), 737 Found is max(Found1, Found2) 738 ; S = Solved, ZStar = Z, Found = 1 739 ) 740 ; ZStar = ZStar0, S = S0, Found = 0 741 ).

751maximize(Z0, S0, S) :- 752 coeff_one(Z0, Z1), 753 functor(S0, F, _), 754 ( F == state -> maximize_mip(Z1, S0, S) 755 ; F == clpr_state -> clpr_maximize(Z1, S0, S) 756 ). 757 758maximize_mip(Z, S0, S) :- 759 maximize_(Z, S0, Solved), 760 state_integrals(S0, Is), 761 ( Is == [] -> S = Solved 762 ; % arrange it so that branch and bound branches on variables 763 % in the same order the integrality constraints were stated in 764 reverse(Is, Is1), 765 state_set_integrals(Is1, S0, S1), 766 ( all_integers(Z, Is1) -> AllInt = 1 767 ; AllInt = 0 768 ), 769 branch_and_bound(Z, Solved, AllInt, _, _, S1, S, 1) 770 ). 771 772all_integers([], _). 773all_integers([Coeff*V|Rest], Is) :- 774 integer(Coeff), 775 memberchk(V, Is), 776 all_integers(Rest, Is).

782minimize(Z0, S0, S) :- 783 coeff_one(Z0, Z1), 784 functor(S0, F, _), 785 ( F == state -> 786 maplist(linsum_negate, Z1, Z2), 787 maximize_mip(Z2, S0, S1), 788 solved_tableau(S1, tableau(Obj, Vars, Inds, Rows)), 789 solved_names(S1, Names), 790 Obj = row(z, Left0, Right0), 791 all_times(Left0, -1, Left), 792 Right is -Right0, 793 Obj1 = row(z, Left, Right), 794 state_integrals(S0, Is), 795 S = solved(tableau(Obj1, Vars, Inds, Rows), Names, Is) 796 ; F == clpr_state -> clpr_minimize(Z1, S0, S) 797 ). 798 799op_pendant(>=, =<). 800op_pendant(=<, >=). 801 802constraints_collapse([]) --> []. 803constraints_collapse([C|Cs]) --> 804 { C = c(Name, Left, Op, Right) }, 805 ( { Name == 0, Left = [1*Var], op_pendant(Op, P) } -> 806 { Pendant = c(0, [1*Var], P, Right) }, 807 ( { select(Pendant, Cs, Rest) } -> 808 [c(0, Left, (=), Right)], 809 { CsLeft = Rest } 810 ; [C], 811 { CsLeft = Cs } 812 ) 813 ; [C], 814 { CsLeft = Cs } 815 ), 816 constraints_collapse(CsLeft). 817 818% solve a (relaxed) LP in standard form 819 820maximize_(Z, S0, S) :- 821 state_constraints(S0, Cs0), 822 phrase(constraints_collapse(Cs0), Cs1), 823 phrase(constraints_normalize(Cs1, Cs, As0), [S0], [S1]), 824 flatten(As0, As1), 825 ( As1 == [] -> 826 make_tableau(Z, Cs, Tableau0), 827 simplex(Tableau0, Tableau), 828 state_names(S1, Names), 829 state_integrals(S1, Is), 830 S = solved(Tableau, Names, Is) 831 ; state_names(S1, Names), 832 state_integrals(S1, Is), 833 two_phase_simplex(Z, Cs, As1, Names, Is, S) 834 ). 835 836make_tableau(Z, Cs, Tableau) :- 837 ZC = c(_, Z, _), 838 phrase(constraints_variables([ZC|Cs]), Variables0), 839 sort(Variables0, Variables), 840 constraints_rows(Cs, Variables, Rows), 841 linsum_row(Variables, Z, Objective1), 842 all_times(Objective1, -1, Obj), 843 length(Variables, LVs), 844 length(Ones, LVs), 845 all_one(Ones), 846 Tableau = tableau(row(z, Obj, 0), Variables, Ones, Rows). 847 848all_one(Ones) :- maplist(=(1), Ones). 849 850proper_form(Variables, Rows, _Coeff*A, Obj0, Obj) :- 851 ( find_row(A, Rows, PivotRow) -> 852 list_first(Variables, A, Col), 853 row_eliminate(Obj0, PivotRow, Col, Obj) 854 ; Obj = Obj0 855 ). 856 857 858two_phase_simplex(Z, Cs, As, Names, Is, S) :- 859 % phase 1: minimize sum of articifial variables 860 make_tableau(As, Cs, Tableau0), 861 Tableau0 = tableau(Obj0, Variables, Inds, Rows), 862 foldl(proper_form(Variables, Rows), As, Obj0, Obj), 863 simplex(tableau(Obj, Variables, Inds, Rows), Tableau1), 864 maplist(var_zero(solved(Tableau1, _, _)), As), 865 % phase 2: remove artificial variables and solve actual LP. 866 tableau_rows(Tableau1, Rows2), 867 eliminate_artificial(As, As, Variables, Rows2, Rows3), 868 list_nths(As, Variables, Nths0), 869 nths_to_zero(Nths0, Inds, Inds1), 870 linsum_row(Variables, Z, Objective), 871 all_times(Objective, -1, Objective1), 872 foldl(proper_form(Variables, Rows3), Z, row(z, Objective1, 0), ObjRow), 873 simplex(tableau(ObjRow, Variables, Inds1, Rows3), Tableau), 874 S = solved(Tableau, Names, Is). 875 876/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 877 If artificial variables are still in the basis, replace them with 878 non-artificial variables if possible. If that is not possible, the 879 constraint is ignored because it is redundant. 880- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 881 882eliminate_artificial([], _, _, Rows, Rows). 883eliminate_artificial([_Coeff*A|Rest], As, Variables, Rows0, Rows) :- 884 ( select(row(A, Left, 0), Rows0, Others) -> 885 ( nth0(Col, Left, Coeff), 886 Coeff =\= 0, 887 nth0(Col, Variables, Var), 888 \+ memberchk(_*Var, As) -> 889 row_divide(row(A, Left, 0), Coeff, Row), 890 gauss_elimination(Rows0, Row, Col, Rows1), 891 swap_basic(Rows1, A, Var, Rows2) 892 ; Rows2 = Others 893 ) 894 ; Rows2 = Rows0 895 ), 896 eliminate_artificial(Rest, As, Variables, Rows2, Rows). 897 898nths_to_zero([], Inds, Inds). 899nths_to_zero([Nth|Nths], Inds0, Inds) :- 900 nth_to_zero(Inds0, 0, Nth, Inds1), 901 nths_to_zero(Nths, Inds1, Inds). 902 903nth_to_zero([], _, _, []). 904nth_to_zero([I|Is], Curr, Nth, [Z|Zs]) :- 905 ( Curr =:= Nth -> [Z|Zs] = [0|Is] 906 ; Z = I, 907 Next is Curr + 1, 908 nth_to_zero(Is, Next, Nth, Zs) 909 ). 910 911 912list_nths([], _, []). 913list_nths([_Coeff*A|As], Variables, [Nth|Nths]) :- 914 list_first(Variables, A, Nth), 915 list_nths(As, Variables, Nths). 916 917 918linsum_negate(Coeff0*Var, Coeff*Var) :- Coeff is -Coeff0. 919 920linsum_row([], _, []). 921linsum_row([V|Vs], Ls, [C|Cs]) :- 922 ( memberchk(Coeff*V, Ls) -> C is rationalize(Coeff) 923 ; C = 0 924 ), 925 linsum_row(Vs, Ls, Cs). 926 927constraints_rows([], _, []). 928constraints_rows([C|Cs], Vars, [R|Rs]) :- 929 C = c(Var, Left0, Right), 930 linsum_row(Vars, Left0, Left), 931 R = row(Var, Left, Right), 932 constraints_rows(Cs, Vars, Rs). 933 934constraints_normalize([], [], []) --> []. 935constraints_normalize([C0|Cs0], [C|Cs], [A|As]) --> 936 { constraint_op(C0, Op), 937 constraint_left(C0, Left), 938 constraint_right(C0, Right), 939 constraint_name(C0, Name), 940 Con =.. [Op, Left, Right] }, 941 constraint_normalize(Con, Name, C, A), 942 constraints_normalize(Cs0, Cs, As). 943 944constraint_normalize(As0 =< B0, Name, c(Slack, [1*Slack|As0], B0), []) --> 945 state_next_var(Slack), 946 state_add_name(Name, Slack). 947constraint_normalize(As0 = B0, Name, c(AID, [1*AID|As0], B0), [-1*AID]) --> 948 state_next_var(AID), 949 state_add_name(Name, AID). 950constraint_normalize(As0 >= B0, Name, c(AID, [-1*Slack,1*AID|As0], B0), [-1*AID]) --> 951 state_next_var(Slack), 952 state_next_var(AID), 953 state_add_name(Name, AID). 954 955coeff_one([], []). 956coeff_one([L|Ls], [Coeff*Var|Rest]) :- 957 ( L = A*B -> Coeff = A, Var = B 958 ; Coeff = 1, Var = L 959 ), 960 coeff_one(Ls, Rest). 961 962 963tableau_optimal(Tableau) :- 964 tableau_objective(Tableau, Objective), 965 tableau_indicators(Tableau, Indicators), 966 Objective = row(_, Left, _), 967 all_nonnegative(Left, Indicators). 968 969all_nonnegative([], []). 970all_nonnegative([Coeff|As], [I|Is]) :- 971 ( I =:= 0 -> true 972 ; Coeff >= 0 973 ), 974 all_nonnegative(As, Is). 975 976pivot_column(Tableau, PCol) :- 977 tableau_objective(Tableau, row(_, Left, _)), 978 tableau_indicators(Tableau, Indicators), 979 first_negative(Left, Indicators, 0, Index0, Val, RestL, RestI), 980 Index1 is Index0 + 1, 981 pivot_column(RestL, RestI, Val, Index1, Index0, PCol). 982 983first_negative([L|Ls], [I|Is], Index0, N, Val, RestL, RestI) :- 984 Index1 is Index0 + 1, 985 ( I =:= 0 -> first_negative(Ls, Is, Index1, N, Val, RestL, RestI) 986 ; ( L < 0 -> N = Index0, Val = L, RestL = Ls, RestI = Is 987 ; first_negative(Ls, Is, Index1, N, Val, RestL, RestI) 988 ) 989 ). 990 991 992pivot_column([], _, _, _, N, N). 993pivot_column([L|Ls], [I|Is], Coeff0, Index0, N0, N) :- 994 ( I =:= 0 -> Coeff1 = Coeff0, N1 = N0 995 ; ( L < Coeff0 -> Coeff1 = L, N1 = Index0 996 ; Coeff1 = Coeff0, N1 = N0 997 ) 998 ), 999 Index1 is Index0 + 1, 1000 pivot_column(Ls, Is, Coeff1, Index1, N1, N). 1001 1002 1003pivot_row(Tableau, PCol, PRow) :- 1004 tableau_rows(Tableau, Rows), 1005 pivot_row(Rows, PCol, false, _, 0, 0, PRow). 1006 1007pivot_row([], _, true, _, _, Row, Row). 1008pivot_row([Row|Rows], PCol, Bounded0, Min0, Index0, PRow0, PRow) :- 1009 Row = row(_Var, Left, B), 1010 nth0(PCol, Left, Ae), 1011 ( Ae > 0 -> 1012 Bounded1 = true, 1013 Bound is B rdiv Ae, 1014 ( Bounded0 == true -> 1015 ( Bound < Min0 -> Min1 = Bound, PRow1 = Index0 1016 ; Min1 = Min0, PRow1 = PRow0 1017 ) 1018 ; Min1 = Bound, PRow1 = Index0 1019 ) 1020 ; Bounded1 = Bounded0, Min1 = Min0, PRow1 = PRow0 1021 ), 1022 Index1 is Index0 + 1, 1023 pivot_row(Rows, PCol, Bounded1, Min1, Index1, PRow1, PRow). 1024 1025 1026row_divide(row(Var, Left0, Right0), Div, row(Var, Left, Right)) :- 1027 all_divide(Left0, Div, Left), 1028 Right is Right0 rdiv Div. 1029 1030 1031all_divide([], _, []). 1032all_divide([R|Rs], Div, [DR|DRs]) :- 1033 DR is R rdiv Div, 1034 all_divide(Rs, Div, DRs). 1035 1036gauss_elimination([], _, _, []). 1037gauss_elimination([Row0|Rows0], PivotRow, Col, [Row|Rows]) :- 1038 PivotRow = row(PVar, _, _), 1039 Row0 = row(Var, _, _), 1040 ( PVar == Var -> Row = PivotRow 1041 ; row_eliminate(Row0, PivotRow, Col, Row) 1042 ), 1043 gauss_elimination(Rows0, PivotRow, Col, Rows). 1044 1045row_eliminate(row(Var, Ls0, R0), row(_, PLs, PR), Col, row(Var, Ls, R)) :- 1046 nth0(Col, Ls0, Coeff), 1047 ( Coeff =:= 0 -> Ls = Ls0, R = R0 1048 ; all_times_minus([PR|PLs], Coeff, [R0|Ls0], [R|Ls]) 1049 ). 1050 1051all_times_minus([], _, _, []). 1052all_times_minus([A|As], T, [B|Bs], [AT|ATs]) :- 1053 ( A == 0 -> AT = B 1054 ; AT is B - A * T 1055 ), 1056 all_times_minus(As, T, Bs, ATs). 1057 1058all_times([], _, []). 1059all_times([A|As], T, [AT|ATs]) :- 1060 AT is A * T, 1061 all_times(As, T, ATs). 1062 1063simplex(Tableau0, Tableau) :- 1064 ( tableau_optimal(Tableau0) -> Tableau0 = Tableau 1065 ; pivot_column(Tableau0, PCol), 1066 pivot_row(Tableau0, PCol, PRow), 1067 Tableau0 = tableau(Obj0,Variables,Inds,Matrix0), 1068 nth0(PRow, Matrix0, Row0), 1069 Row0 = row(Leaving, Left0, _Right0), 1070 nth0(PCol, Left0, PivotElement), 1071 row_divide(Row0, PivotElement, Row1), 1072 gauss_elimination([Obj0|Matrix0], Row1, PCol, [Obj|Matrix1]), 1073 nth0(PCol, Variables, Entering), 1074 swap_basic(Matrix1, Leaving, Entering, Matrix), 1075 simplex(tableau(Obj,Variables,Inds,Matrix), Tableau) 1076 ). 1077 1078swap_basic([Row0|Rows0], Old, New, Matrix) :- 1079 Row0 = row(Var, Left, Right), 1080 ( Var == Old -> Matrix = [row(New, Left, Right)|Rows0] 1081 ; Matrix = [Row0|Rest], 1082 swap_basic(Rows0, Old, New, Rest) 1083 ). 1084 1085constraints_variables([]) --> []. 1086constraints_variables([c(_,Left,_)|Cs]) --> 1087 variables(Left), 1088 constraints_variables(Cs). 1089 1090variables([]) --> []. 1091variables([_Coeff*Var|Rest]) --> [Var], variables(Rest). 1092 1093 1094%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1095%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1096 1097/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1098 A dual algorithm ("algorithm alpha-beta" in Papadimitriou and 1099 Steiglitz) is used for transportation and assignment problems. The 1100 arising max-flow problem is solved with Edmonds-Karp, itself a dual 1101 algorithm. 1102- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1103 1104/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1105 An attributed variable is introduced for each node. Attributes: 1106 node: Original name of the node. 1107 edges: arc_to(To,F,Capacity) (F has an attribute "flow") or 1108 arc_from(From,F,Capacity) 1109 parent: used in breadth-first search 1110- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1111 1112arcs([], Assoc, Assoc). 1113arcs([arc(From0,To0,C)|As], Assoc0, Assoc) :- 1114 ( get_assoc(From0, Assoc0, From) -> Assoc1 = Assoc0 1115 ; put_assoc(From0, Assoc0, From, Assoc1), 1116 put_attr(From, node, From0) 1117 ), 1118 ( get_attr(From, edges, Es) -> true 1119 ; Es = [] 1120 ), 1121 put_attr(F, flow, 0), 1122 put_attr(From, edges, [arc_to(To,F,C)|Es]), 1123 ( get_assoc(To0, Assoc1, To) -> Assoc2 = Assoc1 1124 ; put_assoc(To0, Assoc1, To, Assoc2), 1125 put_attr(To, node, To0) 1126 ), 1127 ( get_attr(To, edges, Es1) -> true 1128 ; Es1 = [] 1129 ), 1130 put_attr(To, edges, [arc_from(From,F,C)|Es1]), 1131 arcs(As, Assoc2, Assoc). 1132 1133 1134edmonds_karp(Arcs0, Arcs) :- 1135 empty_assoc(E), 1136 arcs(Arcs0, E, Assoc), 1137 get_assoc(s, Assoc, S), 1138 get_assoc(t, Assoc, T), 1139 maximum_flow(S, T), 1140 % fetch attvars before deleting visited edges 1141 term_attvars(S, AttVars), 1142 phrase(flow_to_arcs(S), Ls), 1143 arcs_assoc(Ls, Arcs), 1144 maplist(del_attrs, AttVars). 1145 1146flow_to_arcs(V) --> 1147 ( { get_attr(V, edges, Es) } -> 1148 { del_attr(V, edges), 1149 get_attr(V, node, Name) }, 1150 flow_to_arcs_(Es, Name) 1151 ; [] 1152 ). 1153 1154flow_to_arcs_([], _) --> []. 1155flow_to_arcs_([E|Es], Name) --> 1156 edge_to_arc(E, Name), 1157 flow_to_arcs_(Es, Name). 1158 1159edge_to_arc(arc_from(_,_,_), _) --> []. 1160edge_to_arc(arc_to(To,F,C), Name) --> 1161 { get_attr(To, node, NTo), 1162 get_attr(F, flow, Flow) }, 1163 [arc(Name,NTo,Flow,C)], 1164 flow_to_arcs(To). 1165 1166arcs_assoc(Arcs, Hash) :- 1167 empty_assoc(E), 1168 arcs_assoc(Arcs, E, Hash). 1169 1170arcs_assoc([], Hs, Hs). 1171arcs_assoc([arc(From,To,F,C)|Rest], Hs0, Hs) :- 1172 ( get_assoc(From, Hs0, As) -> Hs1 = Hs0 1173 ; put_assoc(From, Hs0, [], Hs1), 1174 empty_assoc(As) 1175 ), 1176 put_assoc(To, As, arc(From,To,F,C), As1), 1177 put_assoc(From, Hs1, As1, Hs2), 1178 arcs_assoc(Rest, Hs2, Hs). 1179 1180 1181/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1182 Strategy: Breadth-first search until we find a free right vertex in 1183 the value graph, then find an augmenting path in reverse. 1184- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1185 1186maximum_flow(S, T) :- 1187 ( augmenting_path([[S]], Levels, T) -> 1188 phrase(augmenting_path(S, T), Path), 1189 Path = [augment(_,First,_)|Rest], 1190 path_minimum(Rest, First, Min), 1191 % format("augmenting path: ~w\n", [Min]), 1192 maplist(augment(Min), Path), 1193 maplist(maplist(clear_parent), Levels), 1194 maximum_flow(S, T) 1195 ; true 1196 ). 1197 1198clear_parent(V) :- del_attr(V, parent). 1199 1200augmenting_path(Levels0, Levels, T) :- 1201 Levels0 = [Vs|_], 1202 Levels1 = [Tos|Levels0], 1203 phrase(reachables(Vs), Tos), 1204 Tos = [_|_], 1205 ( member(To, Tos), To == T -> Levels = Levels1 1206 ; augmenting_path(Levels1, Levels, T) 1207 ). 1208 1209reachables([]) --> []. 1210reachables([V|Vs]) --> 1211 { get_attr(V, edges, Es) }, 1212 reachables_(Es, V), 1213 reachables(Vs). 1214 1215reachables_([], _) --> []. 1216reachables_([E|Es], V) --> 1217 reachable(E, V), 1218 reachables_(Es, V). 1219 1220reachable(arc_from(V,F,_), P) --> 1221 ( { \+ get_attr(V, parent, _), 1222 get_attr(F, flow, Flow), 1223 Flow > 0 } -> 1224 { put_attr(V, parent, P-augment(F,Flow,-)) }, 1225 [V] 1226 ; [] 1227 ). 1228reachable(arc_to(V,F,C), P) --> 1229 ( { \+ get_attr(V, parent, _), 1230 get_attr(F, flow, Flow), 1231 ( C == inf ; Flow < C )} -> 1232 { ( C == inf -> Diff = inf 1233 ; Diff is C - Flow 1234 ), 1235 put_attr(V, parent, P-augment(F,Diff,+)) }, 1236 [V] 1237 ; [] 1238 ). 1239 1240 1241path_minimum([], Min, Min). 1242path_minimum([augment(_,A,_)|As], Min0, Min) :- 1243 ( A == inf -> Min1 = Min0 1244 ; Min1 is min(Min0,A) 1245 ), 1246 path_minimum(As, Min1, Min). 1247 1248augment(Min, augment(F,_,Sign)) :- 1249 get_attr(F, flow, Flow0), 1250 flow_(Sign, Flow0, Min, Flow), 1251 put_attr(F, flow, Flow). 1252 1253flow_(+, F0, A, F) :- F is F0 + A. 1254flow_(-, F0, A, F) :- F is F0 - A. 1255 1256augmenting_path(S, V) --> 1257 ( { V == S } -> [] 1258 ; { get_attr(V, parent, V1-Augment) }, 1259 [Augment], 1260 augmenting_path(S, V1) 1261 ). 1262 1263%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1264 1265naive_init(Supplies, _, Costs, Alphas, Betas) :- 1266 same_length(Supplies, Alphas), 1267 maplist(=(0), Alphas), 1268 lists_transpose(Costs, TCs), 1269 maplist(min_list, TCs, Betas). 1270 1271%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1272 1273lists_transpose([], []). 1274lists_transpose([L|Ls], Ts) :- foldl(transpose_, L, Ts, [L|Ls], _). 1275 1276transpose_(_, Fs, Lists0, Lists) :- 1277 maplist(list_first_rest, Lists0, Fs, Lists). 1278 1279list_first_rest([L|Ls], L, Ls). 1280 1281%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1282/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1283 TODO: use attributed variables throughout 1284- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

1298transportation(Supplies, Demands, Costs, Transport) :- 1299 length(Supplies, LAs), 1300 length(Demands, LBs), 1301 naive_init(Supplies, Demands, Costs, Alphas, Betas), 1302 network_head(Supplies, 1, SArcs, []), 1303 network_tail(Demands, 1, DArcs, []), 1304 numlist(1, LAs, Sources), 1305 numlist(1, LBs, Sinks0), 1306 maplist(make_sink, Sinks0, Sinks), 1307 append(SArcs, DArcs, Torso), 1308 alpha_beta(Torso, Sources, Sinks, Costs, Alphas, Betas, Flow), 1309 flow_transport(Supplies, 1, Demands, Flow, Transport). 1310 1311flow_transport([], _, _, _, []). 1312flow_transport([_|Rest], N, Demands, Flow, [Line|Lines]) :- 1313 transport_line(Demands, N, 1, Flow, Line), 1314 N1 is N + 1, 1315 flow_transport(Rest, N1, Demands, Flow, Lines). 1316 1317transport_line([], _, _, _, []). 1318transport_line([_|Rest], I, J, Flow, [L|Ls]) :- 1319 ( get_assoc(I, Flow, As), get_assoc(p(J), As, arc(I,p(J),F,_)) -> L = F 1320 ; L = 0 1321 ), 1322 J1 is J + 1, 1323 transport_line(Rest, I, J1, Flow, Ls). 1324 1325 1326make_sink(N, p(N)). 1327 1328network_head([], _) --> []. 1329network_head([S|Ss], N) --> 1330 [arc(s,N,S)], 1331 { N1 is N + 1 }, 1332 network_head(Ss, N1). 1333 1334network_tail([], _) --> []. 1335network_tail([D|Ds], N) --> 1336 [arc(p(N),t,D)], 1337 { N1 is N + 1 }, 1338 network_tail(Ds, N1). 1339 1340network_connections([], _, _, _) --> []. 1341network_connections([A|As], Betas, [Cs|Css], N) --> 1342 network_connections(Betas, Cs, A, N, 1), 1343 { N1 is N + 1 }, 1344 network_connections(As, Betas, Css, N1). 1345 1346network_connections([], _, _, _, _) --> []. 1347network_connections([B|Bs], [C|Cs], A, N, PN) --> 1348 ( { C =:= A + B } -> [arc(N,p(PN),inf)] 1349 ; [] 1350 ), 1351 { PN1 is PN + 1 }, 1352 network_connections(Bs, Cs, A, N, PN1). 1353 1354alpha_beta(Torso, Sources, Sinks, Costs, Alphas, Betas, Flow) :- 1355 network_connections(Alphas, Betas, Costs, 1, Cons, []), 1356 append(Torso, Cons, Arcs), 1357 edmonds_karp(Arcs, MaxFlow), 1358 mark_hashes(MaxFlow, MArcs, MRevArcs), 1359 all_markable(MArcs, MRevArcs, Markable), 1360 mark_unmark(Sources, Markable, MarkSources, UnmarkSources), 1361 ( MarkSources == [] -> Flow = MaxFlow 1362 ; mark_unmark(Sinks, Markable, MarkSinks0, UnmarkSinks0), 1363 maplist(un_p, MarkSinks0, MarkSinks), 1364 maplist(un_p, UnmarkSinks0, UnmarkSinks), 1365 MarkSources = [FirstSource|_], 1366 UnmarkSinks = [FirstSink|_], 1367 theta(FirstSource, FirstSink, Costs, Alphas, Betas, TInit), 1368 theta(MarkSources, UnmarkSinks, Costs, Alphas, Betas, TInit, Theta), 1369 duals_add(MarkSources, Alphas, Theta, Alphas1), 1370 duals_add(UnmarkSinks, Betas, Theta, Betas1), 1371 Theta1 is -Theta, 1372 duals_add(UnmarkSources, Alphas1, Theta1, Alphas2), 1373 duals_add(MarkSinks, Betas1, Theta1, Betas2), 1374 alpha_beta(Torso, Sources, Sinks, Costs, Alphas2, Betas2, Flow) 1375 ). 1376 1377mark_hashes(MaxFlow, Arcs, RevArcs) :- 1378 assoc_to_list(MaxFlow, FlowList), 1379 maplist(un_arc, FlowList, FlowList1), 1380 flatten(FlowList1, FlowList2), 1381 empty_assoc(E), 1382 mark_arcs(FlowList2, E, Arcs), 1383 mark_revarcs(FlowList2, E, RevArcs). 1384 1385un_arc(_-Ls0, Ls) :- 1386 assoc_to_list(Ls0, Ls1), 1387 maplist(un_arc_, Ls1, Ls). 1388 1389un_arc_(_-Ls, Ls). 1390 1391mark_arcs([], Arcs, Arcs). 1392mark_arcs([arc(From,To,F,C)|Rest], Arcs0, Arcs) :- 1393 ( get_assoc(From, Arcs0, As) -> true 1394 ; As = [] 1395 ), 1396 ( C == inf -> As1 = [To|As] 1397 ; F < C -> As1 = [To|As] 1398 ; As1 = As 1399 ), 1400 put_assoc(From, Arcs0, As1, Arcs1), 1401 mark_arcs(Rest, Arcs1, Arcs). 1402 1403mark_revarcs([], Arcs, Arcs). 1404mark_revarcs([arc(From,To,F,_)|Rest], Arcs0, Arcs) :- 1405 ( get_assoc(To, Arcs0, Fs) -> true 1406 ; Fs = [] 1407 ), 1408 ( F > 0 -> Fs1 = [From|Fs] 1409 ; Fs1 = Fs 1410 ), 1411 put_assoc(To, Arcs0, Fs1, Arcs1), 1412 mark_revarcs(Rest, Arcs1, Arcs). 1413 1414 1415un_p(p(N), N). 1416 1417duals_add([], Alphas, _, Alphas). 1418duals_add([S|Ss], Alphas0, Theta, Alphas) :- 1419 add_to_nth(1, S, Alphas0, Theta, Alphas1), 1420 duals_add(Ss, Alphas1, Theta, Alphas). 1421 1422add_to_nth(N, N, [A0|As], Theta, [A|As]) :- !, 1423 A is A0 + Theta. 1424add_to_nth(N0, N, [A|As0], Theta, [A|As]) :- 1425 N1 is N0 + 1, 1426 add_to_nth(N1, N, As0, Theta, As). 1427 1428 1429theta(Source, Sink, Costs, Alphas, Betas, Theta) :- 1430 nth1(Source, Costs, Row), 1431 nth1(Sink, Row, C), 1432 nth1(Source, Alphas, A), 1433 nth1(Sink, Betas, B), 1434 Theta is (C - A - B) rdiv 2. 1435 1436theta([], _, _, _, _, Theta, Theta). 1437theta([Source|Sources], Sinks, Costs, Alphas, Betas, Theta0, Theta) :- 1438 theta_(Sinks, Source, Costs, Alphas, Betas, Theta0, Theta1), 1439 theta(Sources, Sinks, Costs, Alphas, Betas, Theta1, Theta). 1440 1441theta_([], _, _, _, _, Theta, Theta). 1442theta_([Sink|Sinks], Source, Costs, Alphas, Betas, Theta0, Theta) :- 1443 theta(Source, Sink, Costs, Alphas, Betas, Theta1), 1444 Theta2 is min(Theta0, Theta1), 1445 theta_(Sinks, Source, Costs, Alphas, Betas, Theta2, Theta). 1446 1447 1448mark_unmark(Nodes, Hash, Mark, Unmark) :- 1449 mark_unmark(Nodes, Hash, Mark, [], Unmark, []). 1450 1451mark_unmark([], _, Mark, Mark, Unmark, Unmark). 1452mark_unmark([Node|Nodes], Markable, Mark0, Mark, Unmark0, Unmark) :- 1453 ( memberchk(Node, Markable) -> 1454 Mark0 = [Node|Mark1], 1455 Unmark0 = Unmark1 1456 ; Mark0 = Mark1, 1457 Unmark0 = [Node|Unmark1] 1458 ), 1459 mark_unmark(Nodes, Markable, Mark1, Mark, Unmark1, Unmark). 1460 1461all_markable(Flow, RevArcs, Markable) :- 1462 phrase(markable(s, [], _, Flow, RevArcs), Markable). 1463 1464all_markable([], Visited, Visited, _, _) --> []. 1465all_markable([To|Tos], Visited0, Visited, Arcs, RevArcs) --> 1466 ( { memberchk(To, Visited0) } -> { Visited0 = Visited1 } 1467 ; markable(To, [To|Visited0], Visited1, Arcs, RevArcs) 1468 ), 1469 all_markable(Tos, Visited1, Visited, Arcs, RevArcs). 1470 1471markable(Current, Visited0, Visited, Arcs, RevArcs) --> 1472 { ( Current = p(_) -> 1473 ( get_assoc(Current, RevArcs, Fs) -> true 1474 ; Fs = [] 1475 ) 1476 ; ( get_assoc(Current, Arcs, Fs) -> true 1477 ; Fs = [] 1478 ) 1479 ) }, 1480 [Current], 1481 all_markable(Fs, [Current|Visited0], Visited, Arcs, RevArcs). 1482 1483/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1484 solve(File) -- read input from File. 1485 1486 Format (NS = number of sources, ND = number of demands): 1487 1488 NS 1489 ND 1490 S1 S2 S3 ... 1491 D1 D2 D3 ... 1492 C11 C12 C13 ... 1493 C21 C22 C23 ... 1494 ... ... ... ... 1495- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1496 1497input(Ss, Ds, Costs) --> 1498 integer(NS), 1499 integer(ND), 1500 n_integers(NS, Ss), 1501 n_integers(ND, Ds), 1502 n_kvectors(NS, ND, Costs). 1503 1504n_kvectors(0, _, []) --> !. 1505n_kvectors(N, K, [V|Vs]) --> 1506 n_integers(K, V), 1507 { N1 is N - 1 }, 1508 n_kvectors(N1, K, Vs). 1509 1510n_integers(0, []) --> !. 1511n_integers(N, [I|Is]) --> integer(I), { N1 is N - 1 }, n_integers(N1, Is). 1512 1513 1514number([D|Ds]) --> digit(D), number(Ds). 1515number([D]) --> digit(D). 1516 1517digit(D) --> [D], { between(0'0, 0'9, D) }. 1518 1519integer(N) --> number(Ds), !, ws, { name(N, Ds) }. 1520 1521ws --> [W], { W =< 0' }, !, ws. % closing quote for syntax highlighting: ' 1522ws --> []. 1523 1524solve(File) :- 1525 time((phrase_from_file(input(Supplies, Demands, Costs), File), 1526 transportation(Supplies, Demands, Costs, Matrix), 1527 maplist(print_row, Matrix))), 1528 halt. 1529 1530print_row(R) :- maplist(print_row_, R), nl. 1531 1532print_row_(N) :- format("~w ", [N]). 1533 1534 1535% ?- call_residue_vars(transportation([12,7,14], [3,15,9,6], [[20,50,10,60],[70,40,60,30],[40,80,70,40]], Ms), Vs). 1536%@ Ms = [[0, 3, 9, 0], [0, 7, 0, 0], [3, 5, 0, 6]], 1537%@ Vs = []. 1538 1539 1540%?- call_residue_vars(simplex:solve('instance_80_80.txt'), Vs). 1541 1542%?- call_residue_vars(simplex:solve('instance_3_4.txt'), Vs). 1543 1544%?- call_residue_vars(simplex:solve('instance_100_100.txt'), Vs). 1545 1546 1547%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

1558% Assignment problem - for now, reduce to transportation problem 1559assignment(Costs, Assignment) :- 1560 length(Costs, LC), 1561 length(Supply, LC), 1562 all_one(Supply), 1563 transportation(Supply, Supply, Costs, Assignment). 1564 1565/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1566 Messages 1567- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1568 1569:- multifile prolog:message//1. 1570 1571prolog:message(simplex(bounded)) --> 1572 ['Using library(simplex) with bounded arithmetic may yield wrong results.'-[]]. 1573 1574warn_if_bounded_arithmetic :- 1575 ( current_prolog_flag(bounded, true) -> 1576 print_message(warning, simplex(bounded)) 1577 ; true 1578 ). 1579 1580:- initialization(warn_if_bounded_arithmetic).