simplex.pl

View source with formatted comments or as raw
    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]).   61
   62/** <module> Solve linear programming problems
   63
   64## Introduction			{#simplex-intro}
   65
   66A *linear programming problem* or simply *linear program* (LP)
   67consists of:
   68
   69  - a set of _linear_ **constraints**
   70  - a set of **variables**
   71  - a _linear_ **objective function**.
   72
   73The goal is to assign values to the variables so as to _maximize_ (or
   74minimize) the value of the objective function while satisfying all
   75constraints.
   76
   77Many optimization problems can be modeled in this way. As one basic
   78example, consider a knapsack with fixed capacity C, and a number of
   79items with sizes `s(i)` and values `v(i)`. The goal is to put as many
   80items as possible in the knapsack (not exceeding its capacity) while
   81maximizing the sum of their values.
   82
   83As another example, suppose you are given a set of _coins_ with
   84certain values, and you are to find the minimum number of coins such
   85that their values sum up to a fixed amount. Instances of these
   86problems are solved below.
   87
   88All numeric quantities are converted to rationals via `rationalize/1`,
   89and rational arithmetic is used throughout solving linear programs. In
   90the current implementation, all variables are implicitly constrained
   91to be _non-negative_. This may change in future versions, and
   92non-negativity constraints should therefore be stated explicitly.
   93
   94
   95## Example 1     {#simplex-ex-1}
   96
   97This is the "radiation therapy" example, taken from _Introduction to
   98Operations Research_ by Hillier and Lieberman.
   99
  100[**Prolog DCG notation**](https://www.metalevel.at/prolog/dcg) is
  101used to _implicitly_ thread the state through posting the constraints:
  102
  103==
  104:- use_module(library(simplex)).
  105
  106radiation(S) :-
  107	gen_state(S0),
  108	post_constraints(S0, S1),
  109	minimize([0.4*x1, 0.5*x2], S1, S).
  110
  111post_constraints -->
  112	constraint([0.3*x1, 0.1*x2] =< 2.7),
  113	constraint([0.5*x1, 0.5*x2] = 6),
  114	constraint([0.6*x1, 0.4*x2] >= 6),
  115	constraint([x1] >= 0),
  116	constraint([x2] >= 0).
  117==
  118
  119An example query:
  120
  121==
  122?- radiation(S), variable_value(S, x1, Val1),
  123		 variable_value(S, x2, Val2).
  124Val1 = 15 rdiv 2,
  125Val2 = 9 rdiv 2.
  126==
  127
  128## Example 2     {#simplex-ex-2}
  129
  130Here is an instance of the knapsack problem described above, where `C
  131= 8`, and we have two types of items: One item with value 7 and size
  1326, and 2 items each having size 4 and value 4. We introduce two
  133variables, `x(1)` and `x(2)` that denote how many items to take of
  134each type.
  135
  136==
  137:- use_module(library(simplex)).
  138
  139knapsack(S) :-
  140	knapsack_constraints(S0),
  141	maximize([7*x(1), 4*x(2)], S0, S).
  142
  143knapsack_constraints(S) :-
  144	gen_state(S0),
  145	constraint([6*x(1), 4*x(2)] =< 8, S0, S1),
  146	constraint([x(1)] =< 1, S1, S2),
  147	constraint([x(2)] =< 2, S2, S).
  148==
  149
  150An example query yields:
  151
  152==
  153?- knapsack(S), variable_value(S, x(1), X1),
  154	        variable_value(S, x(2), X2).
  155X1 = 1
  156X2 = 1 rdiv 2.
  157==
  158
  159That is, we are to take the one item of the first type, and half of one of
  160the items of the other type to maximize the total value of items in the
  161knapsack.
  162
  163If items can not be split, integrality constraints have to be imposed:
  164
  165==
  166knapsack_integral(S) :-
  167	knapsack_constraints(S0),
  168	constraint(integral(x(1)), S0, S1),
  169	constraint(integral(x(2)), S1, S2),
  170	maximize([7*x(1), 4*x(2)], S2, S).
  171==
  172
  173Now the result is different:
  174
  175==
  176?- knapsack_integral(S), variable_value(S, x(1), X1),
  177			 variable_value(S, x(2), X2).
  178
  179X1 = 0
  180X2 = 2
  181==
  182
  183That is, we are to take only the _two_ items of the second type.
  184Notice in particular that always choosing the remaining item with best
  185performance (ratio of value to size) that still fits in the knapsack
  186does not necessarily yield an optimal solution in the presence of
  187integrality constraints.
  188
  189## Example 3      {#simplex-ex-3}
  190
  191We are given:
  192
  193  - 3 coins each worth 1 unit
  194  - 20 coins each worth 5 units and
  195  - 10 coins each worth 20 units.
  196
  197The task is to find a _minimal_ number of these coins that amount to
  198111 units in total. We introduce variables `c(1)`, `c(5)` and `c(20)`
  199denoting how many coins to take of the respective type:
  200
  201==
  202:- use_module(library(simplex)).
  203
  204coins(S) :-
  205	gen_state(S0),
  206	coins(S0, S).
  207
  208coins -->
  209	constraint([c(1), 5*c(5), 20*c(20)] = 111),
  210	constraint([c(1)] =< 3),
  211	constraint([c(5)] =< 20),
  212	constraint([c(20)] =< 10),
  213	constraint([c(1)] >= 0),
  214	constraint([c(5)] >= 0),
  215	constraint([c(20)] >= 0),
  216	constraint(integral(c(1))),
  217	constraint(integral(c(5))),
  218	constraint(integral(c(20))),
  219	minimize([c(1), c(5), c(20)]).
  220==
  221
  222An example query:
  223
  224==
  225?- coins(S), variable_value(S, c(1), C1),
  226	     variable_value(S, c(5), C5),
  227	     variable_value(S, c(20), C20).
  228
  229C1 = 1,
  230C5 = 2,
  231C20 = 5.
  232==
  233
  234@author [Markus Triska](https://www.metalevel.at)
  235*/
  236
  237
  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        ).
  460
  461%% variable_value(+State, +Variable, -Value)
  462%
  463% Value is unified with the value obtained for Variable. State must
  464% correspond to a solved instance.
  465
  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)).
  480
  481%% shadow_price(+State, +Name, -Value)
  482%
  483% Unifies Value with the shadow price corresponding to the linear
  484% constraint whose name is Name. State must correspond to a solved
  485% instance.
  486
  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        ).
  499
  500%% objective(+State, -Objective)
  501%
  502% Unifies Objective with the result of the objective function at the
  503% obtained extremum. State must correspond to a solved instance.
  504
  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).
  544
  545%% gen_state(-State)
  546%
  547% Generates an initial state corresponding to an empty linear program.
  548
  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.
  609
  610%% constraint(+Constraint, +S0, -S)
  611%
  612% Adds a linear or integrality constraint to the linear program
  613% corresponding to state S0. A linear constraint is of the form =|Left
  614% Op C|=, where `Left` is a list of `Coefficient*Variable` terms
  615% (variables in the context of linear programs can be atoms or
  616% compound terms) and `C` is a non-negative numeric constant. The list
  617% represents the sum of its elements. `Op` can be `=`, `=<` or `>=`.
  618% The coefficient `1` can be omitted. An integrality constraint is of
  619% the form integral(Variable) and constrains Variable to an integral
  620% value.
  621
  622
  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        ).
  631
  632%% constraint(+Name, +Constraint, +S0, -S)
  633%
  634% Like constraint/3, and attaches the name Name (an atom or compound
  635% term) to the new constraint.
  636
  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        ).
  651
  652%% constraint_add(+Name, +Left, +S0, -S)
  653%
  654% Left is a list of `Coefficient*Variable` terms. The terms are added
  655% to the left-hand side of the constraint named Name. S is unified
  656% with the resulting state.
  657
  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        ).
  742
  743%% maximize(+Objective, +S0, -S)
  744%
  745% Maximizes the objective function, stated as a list of
  746% `Coefficient*Variable` terms that represents the sum of its
  747% elements, with respect to the linear program corresponding to state
  748% S0. \arg{S} is unified with an internal representation of the solved
  749% instance.
  750
  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).
  777
  778%% minimize(+Objective, +S0, -S)
  779%
  780% Analogous to maximize/3.
  781
  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- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1285
 1286
 1287%% transportation(+Supplies, +Demands, +Costs, -Transport)
 1288%
 1289% Solves a transportation problem. Supplies and Demands must be lists
 1290% of non-negative integers. Their respective sums must be equal. Costs
 1291% is a list of lists representing the cost matrix, where an entry
 1292% (_i_,_j_) denotes the integer cost of transporting one unit from _i_
 1293% to _j_. A transportation plan having minimum cost is computed and
 1294% unified with Transport in the form of a list of lists that
 1295% represents the transportation matrix, where element (_i_,_j_)
 1296% denotes how many units to ship from _i_ to _j_.
 1297
 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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 1548
 1549%% assignment(+Cost, -Assignment)
 1550%
 1551%  Solves a linear assignment problem. Cost is a list of lists
 1552%  representing the quadratic cost matrix, where element (i,j) denotes
 1553%  the integer cost of assigning entity $i$ to entity $j$. An
 1554%  assignment with minimal cost is computed and unified with
 1555%  Assignment as a list of lists, representing an adjacency matrix.
 1556
 1557
 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).