View source with raw 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]).

Solve linear programming problems

Introduction

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

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:

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        ).
 variable_value(+State, +Variable, -Value)
Value is unified with the value obtained for Variable. State must correspond to a solved instance.
  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)).
 shadow_price(+State, +Name, -Value)
Unifies Value with the shadow price corresponding to the linear constraint whose name is Name. State must correspond to a solved instance.
  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        ).
 objective(+State, -Objective)
Unifies Objective with the result of the objective function at the obtained extremum. State must correspond to a solved instance.
  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).
 gen_state(-State)
Generates an initial state corresponding to an empty linear program.
  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.
 constraint(+Constraint, +S0, -S)
Adds a linear or integrality constraint to the linear program corresponding to state S0. A linear constraint is of the form Left Op C, where Left is a list of Coefficient*Variable terms (variables in the context of linear programs can be atoms or compound terms) and C is a non-negative numeric constant. The list represents the sum of its elements. Op can be =, =< or >=. The coefficient 1 can be omitted. An integrality constraint is of the form integral(Variable) and constrains Variable to an integral value.
  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        ).
 constraint(+Name, +Constraint, +S0, -S)
Like constraint/3, and attaches the name Name (an atom or compound term) to the new constraint.
  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        ).
 constraint_add(+Name, +Left, +S0, -S)
Left is a list of Coefficient*Variable terms. The terms are added to the left-hand side of the constraint named Name. S is unified with the resulting state.
  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        ).
 maximize(+Objective, +S0, -S)
Maximizes the objective function, stated as a list of Coefficient*Variable terms that represents the sum of its elements, with respect to the linear program corresponding to state S0. \arg{S} is unified with an internal representation of the solved instance.
  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).
 minimize(+Objective, +S0, -S)
Analogous to maximize/3.
  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        (   member(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([], _, Bounded, _, _, Row, Row) :- Bounded.
 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 ->
 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        ;   MCoeff is -Coeff,
 1049            all_times_plus([PR|PLs], MCoeff, [R0|Ls0], [R|Ls])
 1050        ).
 1051
 1052all_times_plus([], _, _, []).
 1053all_times_plus([A|As], T, [B|Bs], [AT|ATs]) :-
 1054        AT is A * T + B,
 1055        all_times_plus(As, T, Bs, ATs).
 1056
 1057all_times([], _, []).
 1058all_times([A|As], T, [AT|ATs]) :-
 1059        AT is A * T,
 1060        all_times(As, T, ATs).
 1061
 1062simplex(Tableau0, Tableau) :-
 1063        (   tableau_optimal(Tableau0) -> Tableau0 = Tableau
 1064        ;   pivot_column(Tableau0, PCol),
 1065            pivot_row(Tableau0, PCol, PRow),
 1066            Tableau0 = tableau(Obj0,Variables,Inds,Matrix0),
 1067            nth0(PRow, Matrix0, Row0),
 1068            Row0 = row(Leaving, Left0, _Right0),
 1069            nth0(PCol, Left0, PivotElement),
 1070            row_divide(Row0, PivotElement, Row1),
 1071            gauss_elimination([Obj0|Matrix0], Row1, PCol, [Obj|Matrix1]),
 1072            nth0(PCol, Variables, Entering),
 1073            swap_basic(Matrix1, Leaving, Entering, Matrix),
 1074            simplex(tableau(Obj,Variables,Inds,Matrix), Tableau)
 1075        ).
 1076
 1077swap_basic([Row0|Rows0], Old, New, Matrix) :-
 1078        Row0 = row(Var, Left, Right),
 1079        (   Var == Old -> Matrix = [row(New, Left, Right)|Rows0]
 1080        ;   Matrix = [Row0|Rest],
 1081            swap_basic(Rows0, Old, New, Rest)
 1082        ).
 1083
 1084constraints_variables([]) --> [].
 1085constraints_variables([c(_,Left,_)|Cs]) -->
 1086        variables(Left),
 1087        constraints_variables(Cs).
 1088
 1089variables([]) --> [].
 1090variables([_Coeff*Var|Rest]) --> [Var], variables(Rest).
 1091
 1092
 1093%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 1094%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 1095
 1096/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1097   A dual algorithm ("algorithm alpha-beta" in Papadimitriou and
 1098   Steiglitz) is used for transportation and assignment problems. The
 1099   arising max-flow problem is solved with Edmonds-Karp, itself a dual
 1100   algorithm.
 1101- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1102
 1103/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1104   An attributed variable is introduced for each node. Attributes:
 1105   node: Original name of the node.
 1106   edges: arc_to(To,F,Capacity) (F has an attribute "flow") or
 1107          arc_from(From,F,Capacity)
 1108   parent: used in breadth-first search
 1109- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1110
 1111arcs([], Assoc, Assoc).
 1112arcs([arc(From0,To0,C)|As], Assoc0, Assoc) :-
 1113        (   get_assoc(From0, Assoc0, From) -> Assoc1 = Assoc0
 1114        ;   put_assoc(From0, Assoc0, From, Assoc1),
 1115            put_attr(From, node, From0)
 1116        ),
 1117        (   get_attr(From, edges, Es) -> true
 1118        ;   Es = []
 1119        ),
 1120        put_attr(F, flow, 0),
 1121        put_attr(From, edges, [arc_to(To,F,C)|Es]),
 1122        (   get_assoc(To0, Assoc1, To) -> Assoc2 = Assoc1
 1123        ;   put_assoc(To0, Assoc1, To, Assoc2),
 1124            put_attr(To, node, To0)
 1125        ),
 1126        (   get_attr(To, edges, Es1) -> true
 1127        ;   Es1 = []
 1128        ),
 1129        put_attr(To, edges, [arc_from(From,F,C)|Es1]),
 1130        arcs(As, Assoc2, Assoc).
 1131
 1132
 1133edmonds_karp(Arcs0, Arcs) :-
 1134        empty_assoc(E),
 1135        arcs(Arcs0, E, Assoc),
 1136        get_assoc(s, Assoc, S),
 1137        get_assoc(t, Assoc, T),
 1138        maximum_flow(S, T),
 1139        % fetch attvars before deleting visited edges
 1140        term_attvars(S, AttVars),
 1141        phrase(flow_to_arcs(S), Ls),
 1142        arcs_assoc(Ls, Arcs),
 1143        maplist(del_attrs, AttVars).
 1144
 1145flow_to_arcs(V) -->
 1146        (   { get_attr(V, edges, Es) } ->
 1147            { del_attr(V, edges),
 1148              get_attr(V, node, Name) },
 1149            flow_to_arcs_(Es, Name)
 1150        ;   []
 1151        ).
 1152
 1153flow_to_arcs_([], _) --> [].
 1154flow_to_arcs_([E|Es], Name) -->
 1155        edge_to_arc(E, Name),
 1156        flow_to_arcs_(Es, Name).
 1157
 1158edge_to_arc(arc_from(_,_,_), _) --> [].
 1159edge_to_arc(arc_to(To,F,C), Name) -->
 1160        { get_attr(To, node, NTo),
 1161          get_attr(F, flow, Flow) },
 1162        [arc(Name,NTo,Flow,C)],
 1163        flow_to_arcs(To).
 1164
 1165arcs_assoc(Arcs, Hash) :-
 1166        empty_assoc(E),
 1167        arcs_assoc(Arcs, E, Hash).
 1168
 1169arcs_assoc([], Hs, Hs).
 1170arcs_assoc([arc(From,To,F,C)|Rest], Hs0, Hs) :-
 1171        (   get_assoc(From, Hs0, As) -> Hs1 = Hs0
 1172        ;   put_assoc(From, Hs0, [], Hs1),
 1173            empty_assoc(As)
 1174        ),
 1175        put_assoc(To, As, arc(From,To,F,C), As1),
 1176        put_assoc(From, Hs1, As1, Hs2),
 1177        arcs_assoc(Rest, Hs2, Hs).
 1178
 1179
 1180/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1181   Strategy: Breadth-first search until we find a free right vertex in
 1182   the value graph, then find an augmenting path in reverse.
 1183- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1184
 1185maximum_flow(S, T) :-
 1186        (   augmenting_path([[S]], Levels, T) ->
 1187            phrase(augmenting_path(S, T), Path),
 1188            Path = [augment(_,First,_)|Rest],
 1189            path_minimum(Rest, First, Min),
 1190            % format("augmenting path: ~w\n", [Min]),
 1191            maplist(augment(Min), Path),
 1192            maplist(maplist(clear_parent), Levels),
 1193            maximum_flow(S, T)
 1194        ;   true
 1195        ).
 1196
 1197clear_parent(V) :- del_attr(V, parent).
 1198
 1199augmenting_path(Levels0, Levels, T) :-
 1200        Levels0 = [Vs|_],
 1201        Levels1 = [Tos|Levels0],
 1202        phrase(reachables(Vs), Tos),
 1203        Tos = [_|_],
 1204        (   member(To, Tos), To == T -> Levels = Levels1
 1205        ;   augmenting_path(Levels1, Levels, T)
 1206        ).
 1207
 1208reachables([])     --> [].
 1209reachables([V|Vs]) -->
 1210        { get_attr(V, edges, Es) },
 1211        reachables_(Es, V),
 1212        reachables(Vs).
 1213
 1214reachables_([], _)     --> [].
 1215reachables_([E|Es], V) -->
 1216        reachable(E, V),
 1217        reachables_(Es, V).
 1218
 1219reachable(arc_from(V,F,_), P) -->
 1220        (   { \+ get_attr(V, parent, _),
 1221              get_attr(F, flow, Flow),
 1222              Flow > 0 } ->
 1223            { put_attr(V, parent, P-augment(F,Flow,-)) },
 1224            [V]
 1225        ;   []
 1226        ).
 1227reachable(arc_to(V,F,C), P) -->
 1228        (   { \+ get_attr(V, parent, _),
 1229              get_attr(F, flow, Flow),
 1230              (   C == inf ; Flow < C )} ->
 1231            { ( C == inf -> Diff = inf
 1232              ;   Diff is C - Flow
 1233              ),
 1234              put_attr(V, parent, P-augment(F,Diff,+)) },
 1235            [V]
 1236        ;   []
 1237        ).
 1238
 1239
 1240path_minimum([], Min, Min).
 1241path_minimum([augment(_,A,_)|As], Min0, Min) :-
 1242        (   A == inf -> Min1 = Min0
 1243        ;   Min1 is min(Min0,A)
 1244        ),
 1245        path_minimum(As, Min1, Min).
 1246
 1247augment(Min, augment(F,_,Sign)) :-
 1248        get_attr(F, flow, Flow0),
 1249        flow_(Sign, Flow0, Min, Flow),
 1250        put_attr(F, flow, Flow).
 1251
 1252flow_(+, F0, A, F) :- F is F0 + A.
 1253flow_(-, F0, A, F) :- F is F0 - A.
 1254
 1255augmenting_path(S, V) -->
 1256        (   { V == S } -> []
 1257        ;   { get_attr(V, parent, V1-Augment) },
 1258            [Augment],
 1259            augmenting_path(S, V1)
 1260        ).
 1261
 1262%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 1263
 1264naive_init(Supplies, _, Costs, Alphas, Betas) :-
 1265        same_length(Supplies, Alphas),
 1266        maplist(=(0), Alphas),
 1267        lists_transpose(Costs, TCs),
 1268        maplist(min_list, TCs, Betas).
 1269
 1270%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 1271
 1272lists_transpose([], []).
 1273lists_transpose([L|Ls], Ts) :- foldl(transpose_, L, Ts, [L|Ls], _).
 1274
 1275transpose_(_, Fs, Lists0, Lists) :-
 1276        maplist(list_first_rest, Lists0, Fs, Lists).
 1277
 1278list_first_rest([L|Ls], L, Ls).
 1279
 1280%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 1281/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1282   TODO: use attributed variables throughout
 1283- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 transportation(+Supplies, +Demands, +Costs, -Transport)
Solves a transportation problem. Supplies and Demands must be lists of non-negative integers. Their respective sums must be equal. Costs is a list of lists representing the cost matrix, where an entry (i,j) denotes the integer cost of transporting one unit from i to j. A transportation plan having minimum cost is computed and unified with Transport in the form of a list of lists that represents the transportation matrix, where element (i,j) denotes how many units to ship from i to j.
 1297transportation(Supplies, Demands, Costs, Transport) :-
 1298        length(Supplies, LAs),
 1299        length(Demands, LBs),
 1300        naive_init(Supplies, Demands, Costs, Alphas, Betas),
 1301        network_head(Supplies, 1, SArcs, []),
 1302        network_tail(Demands, 1, DArcs, []),
 1303        numlist(1, LAs, Sources),
 1304        numlist(1, LBs, Sinks0),
 1305        maplist(make_sink, Sinks0, Sinks),
 1306        append(SArcs, DArcs, Torso),
 1307        alpha_beta(Torso, Sources, Sinks, Costs, Alphas, Betas, Flow),
 1308        flow_transport(Supplies, 1, Demands, Flow, Transport).
 1309
 1310flow_transport([], _, _, _, []).
 1311flow_transport([_|Rest], N, Demands, Flow, [Line|Lines]) :-
 1312        transport_line(Demands, N, 1, Flow, Line),
 1313        N1 is N + 1,
 1314        flow_transport(Rest, N1, Demands, Flow, Lines).
 1315
 1316transport_line([], _, _, _, []).
 1317transport_line([_|Rest], I, J, Flow, [L|Ls]) :-
 1318        (   get_assoc(I, Flow, As), get_assoc(p(J), As, arc(I,p(J),F,_)) -> L = F
 1319        ;   L = 0
 1320        ),
 1321        J1 is J + 1,
 1322        transport_line(Rest, I, J1, Flow, Ls).
 1323
 1324
 1325make_sink(N, p(N)).
 1326
 1327network_head([], _) --> [].
 1328network_head([S|Ss], N) -->
 1329        [arc(s,N,S)],
 1330        { N1 is N + 1 },
 1331        network_head(Ss, N1).
 1332
 1333network_tail([], _) --> [].
 1334network_tail([D|Ds], N) -->
 1335        [arc(p(N),t,D)],
 1336        { N1 is N + 1 },
 1337        network_tail(Ds, N1).
 1338
 1339network_connections([], _, _, _) --> [].
 1340network_connections([A|As], Betas, [Cs|Css], N) -->
 1341        network_connections(Betas, Cs, A, N, 1),
 1342        { N1 is N + 1 },
 1343        network_connections(As, Betas, Css, N1).
 1344
 1345network_connections([], _, _, _, _) --> [].
 1346network_connections([B|Bs], [C|Cs], A, N, PN) -->
 1347        (   { C =:= A + B } -> [arc(N,p(PN),inf)]
 1348        ;   []
 1349        ),
 1350        { PN1 is PN + 1 },
 1351        network_connections(Bs, Cs, A, N, PN1).
 1352
 1353alpha_beta(Torso, Sources, Sinks, Costs, Alphas, Betas, Flow) :-
 1354        network_connections(Alphas, Betas, Costs, 1, Cons, []),
 1355        append(Torso, Cons, Arcs),
 1356        edmonds_karp(Arcs, MaxFlow),
 1357        mark_hashes(MaxFlow, MArcs, MRevArcs),
 1358        all_markable(MArcs, MRevArcs, Markable),
 1359        mark_unmark(Sources, Markable, MarkSources, UnmarkSources),
 1360        (   MarkSources == [] -> Flow = MaxFlow
 1361        ;   mark_unmark(Sinks, Markable, MarkSinks0, UnmarkSinks0),
 1362            maplist(un_p, MarkSinks0, MarkSinks),
 1363            maplist(un_p, UnmarkSinks0, UnmarkSinks),
 1364            MarkSources = [FirstSource|_],
 1365            UnmarkSinks = [FirstSink|_],
 1366            theta(FirstSource, FirstSink, Costs, Alphas, Betas, TInit),
 1367            theta(MarkSources, UnmarkSinks, Costs, Alphas, Betas, TInit, Theta),
 1368            duals_add(MarkSources, Alphas, Theta, Alphas1),
 1369            duals_add(UnmarkSinks, Betas, Theta, Betas1),
 1370            Theta1 is -Theta,
 1371            duals_add(UnmarkSources, Alphas1, Theta1, Alphas2),
 1372            duals_add(MarkSinks, Betas1, Theta1, Betas2),
 1373            alpha_beta(Torso, Sources, Sinks, Costs, Alphas2, Betas2, Flow)
 1374        ).
 1375
 1376mark_hashes(MaxFlow, Arcs, RevArcs) :-
 1377        assoc_to_list(MaxFlow, FlowList),
 1378        maplist(un_arc, FlowList, FlowList1),
 1379        flatten(FlowList1, FlowList2),
 1380        empty_assoc(E),
 1381        mark_arcs(FlowList2, E, Arcs),
 1382        mark_revarcs(FlowList2, E, RevArcs).
 1383
 1384un_arc(_-Ls0, Ls) :-
 1385        assoc_to_list(Ls0, Ls1),
 1386        maplist(un_arc_, Ls1, Ls).
 1387
 1388un_arc_(_-Ls, Ls).
 1389
 1390mark_arcs([], Arcs, Arcs).
 1391mark_arcs([arc(From,To,F,C)|Rest], Arcs0, Arcs) :-
 1392        (   get_assoc(From, Arcs0, As) -> true
 1393        ;   As = []
 1394        ),
 1395        (   C == inf -> As1 = [To|As]
 1396        ;   F < C -> As1 = [To|As]
 1397        ;   As1 = As
 1398        ),
 1399        put_assoc(From, Arcs0, As1, Arcs1),
 1400        mark_arcs(Rest, Arcs1, Arcs).
 1401
 1402mark_revarcs([], Arcs, Arcs).
 1403mark_revarcs([arc(From,To,F,_)|Rest], Arcs0, Arcs) :-
 1404        (   get_assoc(To, Arcs0, Fs) -> true
 1405        ;   Fs = []
 1406        ),
 1407        (   F > 0 -> Fs1 = [From|Fs]
 1408        ;   Fs1 = Fs
 1409        ),
 1410        put_assoc(To, Arcs0, Fs1, Arcs1),
 1411        mark_revarcs(Rest, Arcs1, Arcs).
 1412
 1413
 1414un_p(p(N), N).
 1415
 1416duals_add([], Alphas, _, Alphas).
 1417duals_add([S|Ss], Alphas0, Theta, Alphas) :-
 1418        add_to_nth(1, S, Alphas0, Theta, Alphas1),
 1419        duals_add(Ss, Alphas1, Theta, Alphas).
 1420
 1421add_to_nth(N, N, [A0|As], Theta, [A|As]) :- !,
 1422        A is A0 + Theta.
 1423add_to_nth(N0, N, [A|As0], Theta, [A|As]) :-
 1424        N1 is N0 + 1,
 1425        add_to_nth(N1, N, As0, Theta, As).
 1426
 1427
 1428theta(Source, Sink, Costs, Alphas, Betas, Theta) :-
 1429        nth1(Source, Costs, Row),
 1430        nth1(Sink, Row, C),
 1431        nth1(Source, Alphas, A),
 1432        nth1(Sink, Betas, B),
 1433        Theta is (C - A - B) rdiv 2.
 1434
 1435theta([], _, _, _, _, Theta, Theta).
 1436theta([Source|Sources], Sinks, Costs, Alphas, Betas, Theta0, Theta) :-
 1437        theta_(Sinks, Source, Costs, Alphas, Betas, Theta0, Theta1),
 1438        theta(Sources, Sinks, Costs, Alphas, Betas, Theta1, Theta).
 1439
 1440theta_([], _, _, _, _, Theta, Theta).
 1441theta_([Sink|Sinks], Source, Costs, Alphas, Betas, Theta0, Theta) :-
 1442        theta(Source, Sink, Costs, Alphas, Betas, Theta1),
 1443        Theta2 is min(Theta0, Theta1),
 1444        theta_(Sinks, Source, Costs, Alphas, Betas, Theta2, Theta).
 1445
 1446
 1447mark_unmark(Nodes, Hash, Mark, Unmark) :-
 1448        mark_unmark(Nodes, Hash, Mark, [], Unmark, []).
 1449
 1450mark_unmark([], _, Mark, Mark, Unmark, Unmark).
 1451mark_unmark([Node|Nodes], Markable, Mark0, Mark, Unmark0, Unmark) :-
 1452        (   memberchk(Node, Markable) ->
 1453            Mark0 = [Node|Mark1],
 1454            Unmark0 = Unmark1
 1455        ;   Mark0 = Mark1,
 1456            Unmark0 = [Node|Unmark1]
 1457        ),
 1458        mark_unmark(Nodes, Markable, Mark1, Mark, Unmark1, Unmark).
 1459
 1460all_markable(Flow, RevArcs, Markable) :-
 1461        phrase(markable(s, [], _, Flow, RevArcs), Markable).
 1462
 1463all_markable([], Visited, Visited, _, _) --> [].
 1464all_markable([To|Tos], Visited0, Visited, Arcs, RevArcs) -->
 1465        (   { memberchk(To, Visited0) } -> { Visited0 = Visited1 }
 1466        ;   markable(To, [To|Visited0], Visited1, Arcs, RevArcs)
 1467        ),
 1468        all_markable(Tos, Visited1, Visited, Arcs, RevArcs).
 1469
 1470markable(Current, Visited0, Visited, Arcs, RevArcs) -->
 1471        { (   Current = p(_) ->
 1472              (   get_assoc(Current, RevArcs, Fs) -> true
 1473              ;   Fs = []
 1474              )
 1475          ;   (   get_assoc(Current, Arcs, Fs) -> true
 1476              ;   Fs = []
 1477              )
 1478          ) },
 1479        [Current],
 1480        all_markable(Fs, [Current|Visited0], Visited, Arcs, RevArcs).
 1481
 1482/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1483   solve(File) -- read input from File.
 1484
 1485   Format (NS = number of sources, ND = number of demands):
 1486
 1487   NS
 1488   ND
 1489   S1 S2 S3 ...
 1490   D1 D2 D3 ...
 1491   C11 C12 C13 ...
 1492   C21 C22 C23 ...
 1493   ... ... ... ...
 1494- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1495
 1496input(Ss, Ds, Costs) -->
 1497        integer(NS),
 1498        integer(ND),
 1499        n_integers(NS, Ss),
 1500        n_integers(ND, Ds),
 1501        n_kvectors(NS, ND, Costs).
 1502
 1503n_kvectors(0, _, []) --> !.
 1504n_kvectors(N, K, [V|Vs]) -->
 1505        n_integers(K, V),
 1506        { N1 is N - 1 },
 1507        n_kvectors(N1, K, Vs).
 1508
 1509n_integers(0, []) --> !.
 1510n_integers(N, [I|Is]) --> integer(I), { N1 is N - 1 }, n_integers(N1, Is).
 1511
 1512
 1513number([D|Ds]) --> digit(D), number(Ds).
 1514number([D])    --> digit(D).
 1515
 1516digit(D) --> [D], { between(0'0, 0'9, D) }.
 1517
 1518integer(N) --> number(Ds), !, ws, { name(N, Ds) }.
 1519
 1520ws --> [W], { W =< 0' }, !, ws.  % closing quote for syntax highlighting: '
 1521ws --> [].
 1522
 1523solve(File) :-
 1524        time((phrase_from_file(input(Supplies, Demands, Costs), File),
 1525              transportation(Supplies, Demands, Costs, Matrix),
 1526              maplist(print_row, Matrix))),
 1527        halt.
 1528
 1529print_row(R) :- maplist(print_row_, R), nl.
 1530
 1531print_row_(N) :- format("~w ", [N]).
 1532
 1533
 1534% ?- 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).
 1535%@ Ms = [[0, 3, 9, 0], [0, 7, 0, 0], [3, 5, 0, 6]],
 1536%@ Vs = [].
 1537
 1538
 1539%?- call_residue_vars(simplex:solve('instance_80_80.txt'), Vs).
 1540
 1541%?- call_residue_vars(simplex:solve('instance_3_4.txt'), Vs).
 1542
 1543%?- call_residue_vars(simplex:solve('instance_100_100.txt'), Vs).
 1544
 1545
 1546%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 assignment(+Cost, -Assignment)
Solves a linear assignment problem. Cost is a list of lists representing the quadratic cost matrix, where element (i,j) denotes the integer cost of assigning entity $i$ to entity $j$. An assignment with minimal cost is computed and unified with Assignment as a list of lists, representing an adjacency matrix.
 1557% Assignment problem - for now, reduce to transportation problem
 1558assignment(Costs, Assignment) :-
 1559        length(Costs, LC),
 1560        length(Supply, LC),
 1561        all_one(Supply),
 1562        transportation(Supply, Supply, Costs, Assignment).
 1563
 1564/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1565   Messages
 1566- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1567
 1568:- multifile prolog:message//1. 1569
 1570prolog:message(simplex(bounded)) -->
 1571        ['Using library(simplex) with bounded arithmetic may yield wrong results.'-[]].
 1572
 1573warn_if_bounded_arithmetic :-
 1574        (   current_prolog_flag(bounded, true) ->
 1575            print_message(warning, simplex(bounded))
 1576        ;   true
 1577        ).
 1578
 1579:- initialization(warn_if_bounded_arithmetic).