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)).

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 */
  235/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  236   CLP(R) bindings
  237   the (unsolved) state is stored as a structure of the form
  238      clpr_state(Options, Cs, Is)
  239   Options: list of Option=Value pairs, currently only eps=Eps
  240   Cs: list of constraints, i.e., structures of the form
  241      c(Name, Left, Op, Right)
  242      anonymous constraints have Name == 0
  243   Is: list of integral variables
  244- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  245
  246gen_state_clpr(State) :- gen_state_clpr([], State).
  247
  248gen_state_clpr(Options, State) :-
  249        ( memberchk(eps=_, Options) -> Options1 = Options
  250        ;   Options1 = [eps=1.0e-6|Options]
  251        ),
  252        State = clpr_state(Options1, [], []).
  253
  254clpr_state_options(clpr_state(Os, _, _), Os).
  255clpr_state_constraints(clpr_state(_, Cs, _), Cs).
  256clpr_state_integrals(clpr_state(_, _, Is), Is).
  257clpr_state_add_integral(I, clpr_state(Os, Cs, Is), clpr_state(Os, Cs, [I|Is])).
  258clpr_state_add_constraint(C, clpr_state(Os, Cs, Is), clpr_state(Os, [C|Cs], Is)).
  259clpr_state_set_constraints(Cs, clpr_state(Os,_,Is), clpr_state(Os, Cs, Is)).
  260
  261clpr_constraint(Name, Constraint, S0, S) :-
  262        (   Constraint = integral(Var) -> clpr_state_add_integral(Var, S0, S)
  263        ;   Constraint =.. [Op,Left,Right],
  264            coeff_one(Left, Left1),
  265            clpr_state_add_constraint(c(Name, Left1, Op, Right), S0, S)
  266        ).
  267
  268clpr_constraint(Constraint, S0, S) :-
  269        clpr_constraint(0, Constraint, S0, S).
  270
  271clpr_shadow_price(clpr_solved(_,_,Duals,_), Name, Value) :-
  272        memberchk(Name-Value0, Duals),
  273        Value is Value0.
  274        %( var(Value0) ->
  275        %       Value = 0
  276        %;
  277        %       Value is Value0
  278        %).
  279
  280
  281
  282clpr_make_variables(Cs, Aliases) :-
  283        clpr_constraints_variables(Cs, Variables0, []),
  284        sort(Variables0, Variables1),
  285        pairs_keys(Aliases, Variables1).
  286
  287clpr_constraints_variables([]) --> [].
  288clpr_constraints_variables([c(_, Left, _, _)|Cs]) -->
  289        variables(Left),
  290        clpr_constraints_variables(Cs).
  291
  292clpr_set_up(Aliases, c(_Name, Left, Op, Right)) :-
  293        clpr_translate_linsum(Left, Aliases, LinSum),
  294        CLPRConstraint =.. [Op, LinSum, Right],
  295        clpr:{ CLPRConstraint }.
  296
  297clpr_set_up_noneg(Aliases, Var) :-
  298        memberchk(Var-CLPVar, Aliases),
  299        { CLPVar >= 0 }.
  300
  301clpr_translate_linsum([], _, 0).
  302clpr_translate_linsum([Coeff*Var|Ls], Aliases, LinSum) :-
  303        memberchk(Var-CLPVar, Aliases),
  304        LinSum = Coeff*CLPVar + LinRest,
  305        clpr_translate_linsum(Ls, Aliases, LinRest).
  306
  307clpr_dual(Objective0, S0, DualValues) :-
  308        clpr_state_constraints(S0, Cs0),
  309        clpr_constraints_variables(Cs0, Variables0, []),
  310        sort(Variables0, Variables1),
  311        maplist(clpr_standard_form, Cs0, Cs1),
  312        clpr_include_all_vars(Cs1, Variables1, Cs2),
  313        clpr_merge_into(Variables1, Objective0, Objective, []),
  314        clpr_unique_names(Cs2, 0, Names),
  315        maplist(clpr_constraint_coefficient, Cs2, Coefficients),
  316        lists_transpose(Coefficients, TCs),
  317        maplist(clpr_dual_constraints(Names), TCs, Objective, DualConstraints),
  318        phrase(clpr_nonneg_constraints(Cs2, Names), DualNonNeg),
  319        append(DualConstraints, DualNonNeg, DualConstraints1),
  320        maplist(clpr_dual_objective, Cs2, Names, DualObjective),
  321        clpr_make_variables(DualConstraints1, Aliases),
  322        maplist(clpr_set_up(Aliases), DualConstraints1),
  323        clpr_translate_linsum(DualObjective, Aliases, LinExpr),
  324        minimize(LinExpr),
  325        Aliases = DualValues.
  326
  327
  328
  329clpr_dual_objective(c(_, _, _, Right), Name, Right*Name).
  330
  331clpr_nonneg_constraints([], _) --> [].
  332clpr_nonneg_constraints([C|Cs], [Name|Names]) -->
  333        { C = c(_, _, Op, _) },
  334        (   { Op == (=<) } -> [c(0, [1*Name], (>=), 0)]
  335        ;   []
  336        ),
  337        clpr_nonneg_constraints(Cs, Names).
  338
  339
  340clpr_dual_constraints(Names, Coeffs, O*_, Constraint) :-
  341        maplist(clpr_dual_linsum, Coeffs, Names, Linsum),
  342        Constraint = c(0, Linsum, (>=), O).
  343
  344clpr_dual_linsum(Coeff, Name, Coeff*Name).
  345
  346clpr_constraint_coefficient(c(_, Left, _, _), Coeff) :-
  347        maplist(coeff_, Left, Coeff).
  348
  349all_coeffs(Coeff*_, Coeff).
  350
  351clpr_unique_names([], _, []).
  352clpr_unique_names([C0|Cs0], Num, [N|Ns]) :-
  353        C0 = c(Name, _, _, _),
  354        (   Name == 0 -> N = Num, Num1 is Num + 1
  355        ;   N = Name, Num1 = Num
  356        ),
  357        clpr_unique_names(Cs0, Num1, Ns).
  358
  359clpr_include_all_vars([], _, []).
  360clpr_include_all_vars([C0|Cs0], Variables, [C|Cs]) :-
  361        C0 = c(Name, Left0, Op, Right),
  362        clpr_merge_into(Variables, Left0, Left, []),
  363        C = c(Name, Left, Op, Right),
  364        clpr_include_all_vars(Cs0, Variables, Cs).
  365
  366clpr_merge_into([], _, Ls, Ls).
  367clpr_merge_into([V|Vs], Left, Ls0, Ls) :-
  368        ( member(Coeff*V, Left) ->
  369                Ls0 = [Coeff*V|Rest]
  370        ;
  371                Ls0 = [0*V|Rest]
  372        ),
  373        clpr_merge_into(Vs, Left, Rest, Ls).
  374
  375
  376
  377
  378clpr_maximize(Expr0, S0, S) :-
  379        coeff_one(Expr0, Expr),
  380        clpr_state_constraints(S0, Cs),
  381        clpr_make_variables(Cs, Aliases),
  382        maplist(clpr_set_up(Aliases), Cs),
  383        clpr_constraints_variables(Cs, Variables0, []),
  384        sort(Variables0, Variables1),
  385        maplist(clpr_set_up_noneg(Aliases), Variables1),
  386        clpr_translate_linsum(Expr, Aliases, LinExpr),
  387        clpr_state_integrals(S0, Is),
  388        ( Is == [] ->
  389                maximize(LinExpr),
  390                Sup is LinExpr,
  391                clpr_dual(Expr, S0, DualValues),
  392                S = clpr_solved(Sup, Aliases, DualValues, S0)
  393        ;
  394                clpr_state_options(S0, Options),
  395                memberchk(eps=Eps, Options),
  396                maplist(clpr_fetch_var(Aliases), Is, Vars),
  397                bb_inf(Vars, -LinExpr, Sup, Vertex, Eps),
  398                pairs_keys_values(Values, Is, Vertex),
  399                % what about the dual in MIPs?
  400                Sup1 is -Sup,
  401                S = clpr_solved(Sup1, Values, [], S0)
  402        ).
  403
  404clpr_minimize(Expr0, S0, S) :-
  405        coeff_one(Expr0, Expr1),
  406        maplist(linsum_negate, Expr1, Expr2),
  407        clpr_maximize(Expr2, S0, S1),
  408        S1 = clpr_solved(Sup, Values, Duals, S0),
  409        Inf is -Sup,
  410        S = clpr_solved(Inf, Values, Duals, S0).
  411
  412clpr_fetch_var(Aliases, Var, X) :- memberchk(Var-X, Aliases).
  413
  414clpr_variable_value(clpr_solved(_, Aliases, _, _), Variable, Value) :-
  415        memberchk(Variable-Value0, Aliases),
  416        Value is Value0.
  417        %( var(Value0) ->
  418        %       Value = 0
  419        %;
  420        %       Value is Value0
  421        %).
  422
  423clpr_objective(clpr_solved(Obj, _, _, _), Obj).
  424
  425clpr_standard_form(c(Name, Left, Op, Right), S) :-
  426        clpr_standard_form_(Op, Name, Left, Right, S).
  427
  428clpr_standard_form_((=), Name, Left, Right, c(Name, Left, (=), Right)).
  429clpr_standard_form_((>=), Name, Left, Right, c(Name, Left1, (=<), Right1)) :-
  430        Right1 is -Right,
  431        maplist(linsum_negate, Left, Left1).
  432clpr_standard_form_((=<), Name, Left, Right, c(Name, Left, (=<), Right)).
  433
  434
  435/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  436   General Simplex Algorithm
  437   Structures used:
  438
  439   tableau(Objective, Variables, Indicators, Constraints)
  440     *) objective function, represented as row
  441     *) list of variables corresponding to columns
  442     *) indicators denoting which variables are still active
  443     *) constraints as rows
  444
  445   row(Var, Left, Right)
  446     *) the basic variable corresponding to this row
  447     *) coefficients of the left-hand side of the constraint
  448     *) right-hand side of the constraint
  449- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  450
  451
  452find_row(Variable, [Row|Rows], R) :-
  453        Row = row(V, _, _),
  454        (   V == Variable -> R = Row
  455        ;   find_row(Variable, Rows, R)
  456        ).
 variable_value(+State, +Variable, -Value)
Value is unified with the value obtained for Variable. State must correspond to a solved instance.
  463variable_value(State, Variable, Value) :-
  464        functor(State, F, _),
  465        (   F == solved ->
  466            solved_tableau(State, Tableau),
  467            tableau_rows(Tableau, Rows),
  468            (   find_row(Variable, Rows, Row) -> Row = row(_, _, Value)
  469            ;   Value = 0
  470            )
  471        ;   F == clpr_solved -> clpr_variable_value(State, Variable, Value)
  472        ).
  473
  474var_zero(State, _Coeff*Var) :- variable_value(State, Var, 0).
  475
  476list_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.
  484shadow_price(State, Name, Value) :-
  485        functor(State, F, _),
  486        (   F == solved ->
  487            solved_tableau(State, Tableau),
  488            tableau_objective(Tableau, row(_,Left,_)),
  489            tableau_variables(Tableau, Variables),
  490            solved_names(State, Names),
  491            memberchk(user(Name)-Var, Names),
  492            list_first(Variables, Var, Nth0),
  493            nth0(Nth0, Left, Value)
  494        ;   F == clpr_solved -> clpr_shadow_price(State, Name, Value)
  495        ).
 objective(+State, -Objective)
Unifies Objective with the result of the objective function at the obtained extremum. State must correspond to a solved instance.
  502objective(State, Obj) :-
  503        functor(State, F, _),
  504        (   F == solved ->
  505            solved_tableau(State, Tableau),
  506            tableau_objective(Tableau, Objective),
  507            Objective = row(_, _, Obj)
  508        ;   clpr_objective(State, Obj)
  509        ).
  510
  511   % interface functions that access tableau components
  512
  513tableau_objective(tableau(Obj, _, _, _), Obj).
  514tableau_rows(tableau(_, _, _, Rows), Rows).
  515tableau_indicators(tableau(_, _, Inds, _), Inds).
  516tableau_variables(tableau(_, Vars, _, _), Vars).
  517
  518/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  519   interface functions that access and modify state components
  520
  521   state is a structure of the form
  522       state(Num, Names, Cs, Is)
  523   Num: used to obtain new unique names for slack variables in a side-effect
  524        free way (increased by one and threaded through)
  525   Names: list of Name-Var, correspondence between constraint-names and
  526        names of slack/artificial variables to obtain shadow prices later
  527   Cs: list of constraints
  528   Is: list of integer variables
  529
  530   constraints are initially represented as c(Name, Left, Op, Right),
  531   and after normalizing as c(Var, Left, Right). Name of unnamed constraints
  532   is 0. The distinction is important for merging constraints (mainly in
  533   branch and bound) with existing ones.
  534- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  535
  536
  537constraint_name(c(Name, _, _, _), Name).
  538constraint_op(c(_, _, Op, _), Op).
  539constraint_left(c(_, Left, _, _), Left).
  540constraint_right(c(_, _, _, Right), Right).
 gen_state(-State)
Generates an initial state corresponding to an empty linear program.
  546gen_state(state(0,[],[],[])).
  547
  548state_add_constraint(C, S0, S) :-
  549        (   constraint_name(C, 0), constraint_left(C, [_Coeff*_Var]) ->
  550            state_merge_constraint(C, S0, S)
  551        ;   state_add_constraint_(C, S0, S)
  552        ).
  553
  554state_add_constraint_(C, state(VID,Ns,Cs,Is), state(VID,Ns,[C|Cs],Is)).
  555
  556state_merge_constraint(C, S0, S) :-
  557        constraint_left(C, [Coeff0*Var0]),
  558        constraint_right(C, Right0),
  559        constraint_op(C, Op),
  560        (   Coeff0 =:= 0 ->
  561            (   Op == (=) -> Right0 =:= 0, S0 = S
  562            ;   Op == (=<) -> S0 = S
  563            ;   Op == (>=) -> Right0 =:= 0, S0 = S
  564            )
  565        ;   Coeff0 < 0 -> state_add_constraint_(C, S0, S)
  566        ;   Right is Right0 rdiv Coeff0,
  567            state_constraints(S0, Cs),
  568            (   select(c(0, [1*Var0], Op, CRight), Cs, RestCs) ->
  569                (   Op == (=) -> CRight =:= Right, S0 = S
  570                ;   Op == (=<) ->
  571                    NewRight is min(Right, CRight),
  572                    NewCs = [c(0, [1*Var0], Op, NewRight)|RestCs],
  573                    state_set_constraints(NewCs, S0, S)
  574                ;   Op == (>=) ->
  575                    NewRight is max(Right, CRight),
  576                    NewCs = [c(0, [1*Var0], Op, NewRight)|RestCs],
  577                    state_set_constraints(NewCs, S0, S)
  578                )
  579            ;   state_add_constraint_(c(0, [1*Var0], Op, Right), S0, S)
  580            )
  581        ).
  582
  583
  584state_add_name(Name, Var), [state(VID,[Name-Var|Ns],Cs,Is)] -->
  585        [state(VID,Ns,Cs,Is)].
  586
  587state_add_integral(Var, state(VID,Ns,Cs,Is), state(VID,Ns,Cs,[Var|Is])).
  588
  589state_constraints(state(_, _, Cs, _), Cs).
  590state_names(state(_,Names,_,_), Names).
  591state_integrals(state(_,_,_,Is), Is).
  592state_set_constraints(Cs, state(VID,Ns,_,Is), state(VID,Ns,Cs,Is)).
  593state_set_integrals(Is, state(VID,Ns,Cs,_), state(VID,Ns,Cs,Is)).
  594
  595
  596state_next_var(VarID0), [state(VarID1,Names,Cs,Is)] -->
  597        [state(VarID0,Names,Cs,Is)],
  598        { VarID1 is VarID0 + 1 }.
  599
  600solved_tableau(solved(Tableau, _, _), Tableau).
  601solved_names(solved(_, Names,_), Names).
  602solved_integrals(solved(_,_,Is), Is).
  603
  604% User-named constraints are wrapped with user/1 to also allow "0" in
  605% 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.
  620constraint(C, S0, S) :-
  621        functor(S0, F, _),
  622        (   F == state ->
  623            (   C = integral(Var) -> state_add_integral(Var, S0, S)
  624            ;   constraint_(0, C, S0, S)
  625            )
  626        ;   F == clpr_state -> clpr_constraint(C, S0, S)
  627        ).
 constraint(+Name, +Constraint, +S0, -S)
Like constraint/3, and attaches the name Name (an atom or compound term) to the new constraint.
  634constraint(Name, C, S0, S) :- constraint_(user(Name), C, S0, S).
  635
  636constraint_(Name, C, S0, S) :-
  637        functor(S0, F, _),
  638        (   F == state ->
  639            (   C = integral(Var) -> state_add_integral(Var, S0, S)
  640            ;   C =.. [Op, Left0, Right0],
  641                coeff_one(Left0, Left),
  642                Right0 >= 0,
  643                Right is rationalize(Right0),
  644                state_add_constraint(c(Name, Left, Op, Right), S0, S)
  645            )
  646        ;   F == clpr_state -> clpr_constraint(Name, C, S0, S)
  647        ).
 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.
  655constraint_add(Name, A, S0, S) :-
  656        functor(S0, F, _),
  657        (   F == state ->
  658            state_constraints(S0, Cs),
  659            add_left(Cs, user(Name), A, Cs1),
  660            state_set_constraints(Cs1, S0, S)
  661        ;   F == clpr_state ->
  662            clpr_state_constraints(S0, Cs),
  663            add_left(Cs, Name, A, Cs1),
  664            clpr_state_set_constraints(Cs1, S0, S)
  665        ).
  666
  667
  668add_left([c(Name,Left0,Op,Right)|Cs], V, A, [c(Name,Left,Op,Right)|Rest]) :-
  669        (   Name == V -> append(A, Left0, Left), Rest = Cs
  670        ;   Left0 = Left, add_left(Cs, V, A, Rest)
  671        ).
  672
  673branching_variable(State, Variable) :-
  674        solved_integrals(State, Integrals),
  675        member(Variable, Integrals),
  676        variable_value(State, Variable, Value),
  677        \+ integer(Value).
  678
  679
  680worth_investigating(ZStar0, _, _) :- var(ZStar0).
  681worth_investigating(ZStar0, AllInt, Z) :-
  682        nonvar(ZStar0),
  683        (   AllInt =:= 1 -> Z1 is floor(Z)
  684        ;   Z1 = Z
  685        ),
  686        Z1 > ZStar0.
  687
  688
  689branch_and_bound(Objective, Solved, AllInt, ZStar0, ZStar, S0, S, Found) :-
  690        objective(Solved, Z),
  691        (   worth_investigating(ZStar0, AllInt, Z) ->
  692            (   branching_variable(Solved, BrVar) ->
  693                variable_value(Solved, BrVar, Value),
  694                Value1 is floor(Value),
  695                Value2 is Value1 + 1,
  696                constraint([BrVar] =< Value1, S0, SubProb1),
  697                (   maximize_(Objective, SubProb1, SubSolved1) ->
  698                    Sub1Feasible = 1,
  699                    objective(SubSolved1, Obj1)
  700                ;   Sub1Feasible = 0
  701                ),
  702                constraint([BrVar] >= Value2, S0, SubProb2),
  703                (   maximize_(Objective, SubProb2, SubSolved2) ->
  704                    Sub2Feasible = 1,
  705                    objective(SubSolved2, Obj2)
  706                ;   Sub2Feasible = 0
  707                ),
  708                (   Sub1Feasible =:= 1, Sub2Feasible =:= 1 ->
  709                    (   Obj1 >= Obj2 ->
  710                        First = SubProb1,
  711                        Second = SubProb2,
  712                        FirstSolved = SubSolved1,
  713                        SecondSolved = SubSolved2
  714                    ;   First = SubProb2,
  715                        Second = SubProb1,
  716                        FirstSolved = SubSolved2,
  717                        SecondSolved = SubSolved1
  718                    ),
  719                    branch_and_bound(Objective, FirstSolved, AllInt, ZStar0, ZStar1, First, Solved1, Found1),
  720                    branch_and_bound(Objective, SecondSolved, AllInt, ZStar1, ZStar2, Second, Solved2, Found2)
  721                ;   Sub1Feasible =:= 1 ->
  722                    branch_and_bound(Objective, SubSolved1, AllInt, ZStar0, ZStar1, SubProb1, Solved1, Found1),
  723                    Found2 = 0
  724                ;   Sub2Feasible =:= 1 ->
  725                    Found1 = 0,
  726                    branch_and_bound(Objective, SubSolved2, AllInt, ZStar0, ZStar2, SubProb2, Solved2, Found2)
  727                ;   Found1 = 0, Found2 = 0
  728                ),
  729                (   Found1 =:= 1, Found2 =:= 1 -> S = Solved2, ZStar = ZStar2
  730                ;   Found1 =:= 1 -> S = Solved1, ZStar = ZStar1
  731                ;   Found2 =:= 1 -> S = Solved2, ZStar = ZStar2
  732                ;   S = S0, ZStar = ZStar0
  733                ),
  734                Found is max(Found1, Found2)
  735            ;   S = Solved, ZStar = Z, Found = 1
  736            )
  737        ;   ZStar = ZStar0, S = S0, Found = 0
  738        ).
 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.
  748maximize(Z0, S0, S) :-
  749        coeff_one(Z0, Z1),
  750        functor(S0, F, _),
  751        (   F == state -> maximize_mip(Z1, S0, S)
  752        ;   F == clpr_state -> clpr_maximize(Z1, S0, S)
  753        ).
  754
  755maximize_mip(Z, S0, S) :-
  756        maximize_(Z, S0, Solved),
  757        state_integrals(S0, Is),
  758        (   Is == [] -> S = Solved
  759        ;   % arrange it so that branch and bound branches on variables
  760            % in the same order the integrality constraints were stated in
  761            reverse(Is, Is1),
  762            state_set_integrals(Is1, S0, S1),
  763            (   all_integers(Z, Is1) -> AllInt = 1
  764            ;   AllInt = 0
  765            ),
  766            branch_and_bound(Z, Solved, AllInt, _, _, S1, S, 1)
  767        ).
  768
  769all_integers([], _).
  770all_integers([Coeff*V|Rest], Is) :-
  771        integer(Coeff),
  772        memberchk(V, Is),
  773        all_integers(Rest, Is).
 minimize(+Objective, +S0, -S)
Analogous to maximize/3.
  779minimize(Z0, S0, S) :-
  780        coeff_one(Z0, Z1),
  781        functor(S0, F, _),
  782        (   F == state ->
  783            maplist(linsum_negate, Z1, Z2),
  784            maximize_mip(Z2, S0, S1),
  785            solved_tableau(S1, tableau(Obj, Vars, Inds, Rows)),
  786            solved_names(S1, Names),
  787            Obj = row(z, Left0, Right0),
  788            all_times(Left0, -1, Left),
  789            Right is -Right0,
  790            Obj1 = row(z, Left, Right),
  791            state_integrals(S0, Is),
  792            S = solved(tableau(Obj1, Vars, Inds, Rows), Names, Is)
  793        ;   F == clpr_state -> clpr_minimize(Z1, S0, S)
  794        ).
  795
  796op_pendant(>=, =<).
  797op_pendant(=<, >=).
  798
  799constraints_collapse([]) --> [].
  800constraints_collapse([C|Cs]) -->
  801        { C = c(Name, Left, Op, Right) },
  802        (   { Name == 0, Left = [1*Var], op_pendant(Op, P) } ->
  803            { Pendant = c(0, [1*Var], P, Right) },
  804            (   { select(Pendant, Cs, Rest) } ->
  805                [c(0, Left, (=), Right)],
  806                { CsLeft = Rest }
  807            ;   [C],
  808                { CsLeft = Cs }
  809            )
  810        ;   [C],
  811            { CsLeft = Cs }
  812        ),
  813        constraints_collapse(CsLeft).
  814
  815% solve a (relaxed) LP in standard form
  816
  817maximize_(Z, S0, S) :-
  818        state_constraints(S0, Cs0),
  819        phrase(constraints_collapse(Cs0), Cs1),
  820        phrase(constraints_normalize(Cs1, Cs, As0), [S0], [S1]),
  821        flatten(As0, As1),
  822        (   As1 == [] ->
  823            make_tableau(Z, Cs, Tableau0),
  824            simplex(Tableau0, Tableau),
  825            state_names(S1, Names),
  826            state_integrals(S1, Is),
  827            S = solved(Tableau, Names, Is)
  828        ;   state_names(S1, Names),
  829            state_integrals(S1, Is),
  830            two_phase_simplex(Z, Cs, As1, Names, Is, S)
  831        ).
  832
  833make_tableau(Z, Cs, Tableau) :-
  834        ZC = c(_, Z, _),
  835        phrase(constraints_variables([ZC|Cs]), Variables0),
  836        sort(Variables0, Variables),
  837        constraints_rows(Cs, Variables, Rows),
  838        linsum_row(Variables, Z, Objective1),
  839        all_times(Objective1, -1, Obj),
  840        length(Variables, LVs),
  841        length(Ones, LVs),
  842        all_one(Ones),
  843        Tableau = tableau(row(z, Obj, 0), Variables, Ones, Rows).
  844
  845all_one(Ones) :- maplist(=(1), Ones).
  846
  847proper_form(Variables, Rows, _Coeff*A, Obj0, Obj) :-
  848        (   find_row(A, Rows, PivotRow) ->
  849            list_first(Variables, A, Col),
  850            row_eliminate(Obj0, PivotRow, Col, Obj)
  851        ;   Obj = Obj0
  852        ).
  853
  854
  855two_phase_simplex(Z, Cs, As, Names, Is, S) :-
  856        % phase 1: minimize sum of articifial variables
  857        make_tableau(As, Cs, Tableau0),
  858        Tableau0 = tableau(Obj0, Variables, Inds, Rows),
  859        foldl(proper_form(Variables, Rows), As, Obj0, Obj),
  860        simplex(tableau(Obj, Variables, Inds, Rows), Tableau1),
  861        maplist(var_zero(solved(Tableau1, _, _)), As),
  862        % phase 2: remove artificial variables and solve actual LP.
  863        tableau_rows(Tableau1, Rows2),
  864        eliminate_artificial(As, As, Variables, Rows2, Rows3),
  865        list_nths(As, Variables, Nths0),
  866        nths_to_zero(Nths0, Inds, Inds1),
  867        linsum_row(Variables, Z, Objective),
  868        all_times(Objective, -1, Objective1),
  869        foldl(proper_form(Variables, Rows3), Z, row(z, Objective1, 0), ObjRow),
  870        simplex(tableau(ObjRow, Variables, Inds1, Rows3), Tableau),
  871        S = solved(Tableau, Names, Is).
  872
  873/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  874   If artificial variables are still in the basis, replace them with
  875   non-artificial variables if possible. If that is not possible, the
  876   constraint is ignored because it is redundant.
  877- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  878
  879eliminate_artificial([], _, _, Rows, Rows).
  880eliminate_artificial([_Coeff*A|Rest], As, Variables, Rows0, Rows) :-
  881        (   select(row(A, Left, 0), Rows0, Others) ->
  882            (   nth0(Col, Left, Coeff),
  883                Coeff =\= 0,
  884                nth0(Col, Variables, Var),
  885                \+ memberchk(_*Var, As) ->
  886                row_divide(row(A, Left, 0), Coeff, Row),
  887                gauss_elimination(Rows0, Row, Col, Rows1),
  888                swap_basic(Rows1, A, Var, Rows2)
  889            ;   Rows2 = Others
  890            )
  891        ;   Rows2 = Rows0
  892        ),
  893        eliminate_artificial(Rest, As, Variables, Rows2, Rows).
  894
  895nths_to_zero([], Inds, Inds).
  896nths_to_zero([Nth|Nths], Inds0, Inds) :-
  897        nth_to_zero(Inds0, 0, Nth, Inds1),
  898        nths_to_zero(Nths, Inds1, Inds).
  899
  900nth_to_zero([], _, _, []).
  901nth_to_zero([I|Is], Curr, Nth, [Z|Zs]) :-
  902        (   Curr =:= Nth -> [Z|Zs] = [0|Is]
  903        ;   Z = I,
  904            Next is Curr + 1,
  905            nth_to_zero(Is, Next, Nth, Zs)
  906        ).
  907
  908
  909list_nths([], _, []).
  910list_nths([_Coeff*A|As], Variables, [Nth|Nths]) :-
  911        list_first(Variables, A, Nth),
  912        list_nths(As, Variables, Nths).
  913
  914
  915linsum_negate(Coeff0*Var, Coeff*Var) :- Coeff is -Coeff0.
  916
  917linsum_row([], _, []).
  918linsum_row([V|Vs], Ls, [C|Cs]) :-
  919        (   member(Coeff*V, Ls) -> C is rationalize(Coeff)
  920        ;   C = 0
  921        ),
  922        linsum_row(Vs, Ls, Cs).
  923
  924constraints_rows([], _, []).
  925constraints_rows([C|Cs], Vars, [R|Rs]) :-
  926        C = c(Var, Left0, Right),
  927        linsum_row(Vars, Left0, Left),
  928        R = row(Var, Left, Right),
  929        constraints_rows(Cs, Vars, Rs).
  930
  931constraints_normalize([], [], []) --> [].
  932constraints_normalize([C0|Cs0], [C|Cs], [A|As]) -->
  933        { constraint_op(C0, Op),
  934          constraint_left(C0, Left),
  935          constraint_right(C0, Right),
  936          constraint_name(C0, Name),
  937          Con =.. [Op, Left, Right] },
  938        constraint_normalize(Con, Name, C, A),
  939        constraints_normalize(Cs0, Cs, As).
  940
  941constraint_normalize(As0 =< B0, Name, c(Slack, [1*Slack|As0], B0), []) -->
  942        state_next_var(Slack),
  943        state_add_name(Name, Slack).
  944constraint_normalize(As0 = B0, Name, c(AID, [1*AID|As0], B0), [-1*AID]) -->
  945        state_next_var(AID),
  946        state_add_name(Name, AID).
  947constraint_normalize(As0 >= B0, Name, c(AID, [-1*Slack,1*AID|As0], B0), [-1*AID]) -->
  948        state_next_var(Slack),
  949        state_next_var(AID),
  950        state_add_name(Name, AID).
  951
  952coeff_one([], []).
  953coeff_one([L|Ls], [Coeff*Var|Rest]) :-
  954        (   L = A*B -> Coeff = A, Var = B
  955        ;   Coeff = 1, Var = L
  956        ),
  957        coeff_one(Ls, Rest).
  958
  959
  960tableau_optimal(Tableau) :-
  961        tableau_objective(Tableau, Objective),
  962        tableau_indicators(Tableau, Indicators),
  963        Objective = row(_, Left, _),
  964        all_nonnegative(Left, Indicators).
  965
  966all_nonnegative([], []).
  967all_nonnegative([Coeff|As], [I|Is]) :-
  968        (   I =:= 0 -> true
  969        ;   Coeff >= 0
  970        ),
  971        all_nonnegative(As, Is).
  972
  973pivot_column(Tableau, PCol) :-
  974        tableau_objective(Tableau, row(_, Left, _)),
  975        tableau_indicators(Tableau, Indicators),
  976        first_negative(Left, Indicators, 0, Index0, Val, RestL, RestI),
  977        Index1 is Index0 + 1,
  978        pivot_column(RestL, RestI, Val, Index1, Index0, PCol).
  979
  980first_negative([L|Ls], [I|Is], Index0, N, Val, RestL, RestI) :-
  981        Index1 is Index0 + 1,
  982        (   I =:= 0 -> first_negative(Ls, Is, Index1, N, Val, RestL, RestI)
  983        ;   (   L < 0 -> N = Index0, Val = L, RestL = Ls, RestI = Is
  984            ;   first_negative(Ls, Is, Index1, N, Val, RestL, RestI)
  985            )
  986        ).
  987
  988
  989pivot_column([], _, _, _, N, N).
  990pivot_column([L|Ls], [I|Is], Coeff0, Index0, N0, N) :-
  991        (   I =:= 0 -> Coeff1 = Coeff0, N1 = N0
  992        ;   (   L < Coeff0 -> Coeff1 = L, N1 = Index0
  993            ;   Coeff1 = Coeff0, N1 = N0
  994            )
  995        ),
  996        Index1 is Index0 + 1,
  997        pivot_column(Ls, Is, Coeff1, Index1, N1, N).
  998
  999
 1000pivot_row(Tableau, PCol, PRow) :-
 1001        tableau_rows(Tableau, Rows),
 1002        pivot_row(Rows, PCol, false, _, 0, 0, PRow).
 1003
 1004pivot_row([], _, Bounded, _, _, Row, Row) :- Bounded.
 1005pivot_row([Row|Rows], PCol, Bounded0, Min0, Index0, PRow0, PRow) :-
 1006        Row = row(_Var, Left, B),
 1007        nth0(PCol, Left, Ae),
 1008        (   Ae > 0 ->
 1009            Bounded1 = true,
 1010            Bound is B rdiv Ae,
 1011            (   Bounded0 ->
 1012                (   Bound < Min0 -> Min1 = Bound, PRow1 = Index0
 1013                ;   Min1 = Min0, PRow1 = PRow0
 1014                )
 1015            ;   Min1 = Bound, PRow1 = Index0
 1016            )
 1017        ;   Bounded1 = Bounded0, Min1 = Min0, PRow1 = PRow0
 1018        ),
 1019        Index1 is Index0 + 1,
 1020        pivot_row(Rows, PCol, Bounded1, Min1, Index1, PRow1, PRow).
 1021
 1022
 1023row_divide(row(Var, Left0, Right0), Div, row(Var, Left, Right)) :-
 1024        all_divide(Left0, Div, Left),
 1025        Right is Right0 rdiv Div.
 1026
 1027
 1028all_divide([], _, []).
 1029all_divide([R|Rs], Div, [DR|DRs]) :-
 1030        DR is R rdiv Div,
 1031        all_divide(Rs, Div, DRs).
 1032
 1033gauss_elimination([], _, _, []).
 1034gauss_elimination([Row0|Rows0], PivotRow, Col, [Row|Rows]) :-
 1035        PivotRow = row(PVar, _, _),
 1036        Row0 = row(Var, _, _),
 1037        (   PVar == Var -> Row = PivotRow
 1038        ;   row_eliminate(Row0, PivotRow, Col, Row)
 1039        ),
 1040        gauss_elimination(Rows0, PivotRow, Col, Rows).
 1041
 1042row_eliminate(row(Var, Ls0, R0), row(_, PLs, PR), Col, row(Var, Ls, R)) :-
 1043        nth0(Col, Ls0, Coeff),
 1044        (   Coeff =:= 0 -> Ls = Ls0, R = R0
 1045        ;   MCoeff is -Coeff,
 1046            all_times_plus([PR|PLs], MCoeff, [R0|Ls0], [R|Ls])
 1047        ).
 1048
 1049all_times_plus([], _, _, []).
 1050all_times_plus([A|As], T, [B|Bs], [AT|ATs]) :-
 1051        AT is A * T + B,
 1052        all_times_plus(As, T, Bs, ATs).
 1053
 1054all_times([], _, []).
 1055all_times([A|As], T, [AT|ATs]) :-
 1056        AT is A * T,
 1057        all_times(As, T, ATs).
 1058
 1059simplex(Tableau0, Tableau) :-
 1060        (   tableau_optimal(Tableau0) -> Tableau0 = Tableau
 1061        ;   pivot_column(Tableau0, PCol),
 1062            pivot_row(Tableau0, PCol, PRow),
 1063            Tableau0 = tableau(Obj0,Variables,Inds,Matrix0),
 1064            nth0(PRow, Matrix0, Row0),
 1065            Row0 = row(Leaving, Left0, _Right0),
 1066            nth0(PCol, Left0, PivotElement),
 1067            row_divide(Row0, PivotElement, Row1),
 1068            gauss_elimination([Obj0|Matrix0], Row1, PCol, [Obj|Matrix1]),
 1069            nth0(PCol, Variables, Entering),
 1070            swap_basic(Matrix1, Leaving, Entering, Matrix),
 1071            simplex(tableau(Obj,Variables,Inds,Matrix), Tableau)
 1072        ).
 1073
 1074swap_basic([Row0|Rows0], Old, New, Matrix) :-
 1075        Row0 = row(Var, Left, Right),
 1076        (   Var == Old -> Matrix = [row(New, Left, Right)|Rows0]
 1077        ;   Matrix = [Row0|Rest],
 1078            swap_basic(Rows0, Old, New, Rest)
 1079        ).
 1080
 1081constraints_variables([]) --> [].
 1082constraints_variables([c(_,Left,_)|Cs]) -->
 1083        variables(Left),
 1084        constraints_variables(Cs).
 1085
 1086variables([]) --> [].
 1087variables([_Coeff*Var|Rest]) --> [Var], variables(Rest).
 1088
 1089
 1090%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 1091%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 1092
 1093/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1094   A dual algorithm ("algorithm alpha-beta" in Papadimitriou and
 1095   Steiglitz) is used for transportation and assignment problems. The
 1096   arising max-flow problem is solved with Edmonds-Karp, itself a dual
 1097   algorithm.
 1098- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1099
 1100/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1101   An attributed variable is introduced for each node. Attributes:
 1102   node: Original name of the node.
 1103   edges: arc_to(To,F,Capacity) (F has an attribute "flow") or
 1104          arc_from(From,F,Capacity)
 1105   parent: used in breadth-first search
 1106- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1107
 1108arcs([], Assoc, Assoc).
 1109arcs([arc(From0,To0,C)|As], Assoc0, Assoc) :-
 1110        (   get_assoc(From0, Assoc0, From) -> Assoc1 = Assoc0
 1111        ;   put_assoc(From0, Assoc0, From, Assoc1),
 1112            put_attr(From, node, From0)
 1113        ),
 1114        (   get_attr(From, edges, Es) -> true
 1115        ;   Es = []
 1116        ),
 1117        put_attr(F, flow, 0),
 1118        put_attr(From, edges, [arc_to(To,F,C)|Es]),
 1119        (   get_assoc(To0, Assoc1, To) -> Assoc2 = Assoc1
 1120        ;   put_assoc(To0, Assoc1, To, Assoc2),
 1121            put_attr(To, node, To0)
 1122        ),
 1123        (   get_attr(To, edges, Es1) -> true
 1124        ;   Es1 = []
 1125        ),
 1126        put_attr(To, edges, [arc_from(From,F,C)|Es1]),
 1127        arcs(As, Assoc2, Assoc).
 1128
 1129
 1130edmonds_karp(Arcs0, Arcs) :-
 1131        empty_assoc(E),
 1132        arcs(Arcs0, E, Assoc),
 1133        get_assoc(s, Assoc, S),
 1134        get_assoc(t, Assoc, T),
 1135        maximum_flow(S, T),
 1136        % fetch attvars before deleting visited edges
 1137        term_attvars(S, AttVars),
 1138        phrase(flow_to_arcs(S), Ls),
 1139        arcs_assoc(Ls, Arcs),
 1140        maplist(del_attrs, AttVars).
 1141
 1142flow_to_arcs(V) -->
 1143        (   { get_attr(V, edges, Es) } ->
 1144            { del_attr(V, edges),
 1145              get_attr(V, node, Name) },
 1146            flow_to_arcs_(Es, Name)
 1147        ;   []
 1148        ).
 1149
 1150flow_to_arcs_([], _) --> [].
 1151flow_to_arcs_([E|Es], Name) -->
 1152        edge_to_arc(E, Name),
 1153        flow_to_arcs_(Es, Name).
 1154
 1155edge_to_arc(arc_from(_,_,_), _) --> [].
 1156edge_to_arc(arc_to(To,F,C), Name) -->
 1157        { get_attr(To, node, NTo),
 1158          get_attr(F, flow, Flow) },
 1159        [arc(Name,NTo,Flow,C)],
 1160        flow_to_arcs(To).
 1161
 1162arcs_assoc(Arcs, Hash) :-
 1163        empty_assoc(E),
 1164        arcs_assoc(Arcs, E, Hash).
 1165
 1166arcs_assoc([], Hs, Hs).
 1167arcs_assoc([arc(From,To,F,C)|Rest], Hs0, Hs) :-
 1168        (   get_assoc(From, Hs0, As) -> Hs1 = Hs0
 1169        ;   put_assoc(From, Hs0, [], Hs1),
 1170            empty_assoc(As)
 1171        ),
 1172        put_assoc(To, As, arc(From,To,F,C), As1),
 1173        put_assoc(From, Hs1, As1, Hs2),
 1174        arcs_assoc(Rest, Hs2, Hs).
 1175
 1176
 1177/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1178   Strategy: Breadth-first search until we find a free right vertex in
 1179   the value graph, then find an augmenting path in reverse.
 1180- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1181
 1182maximum_flow(S, T) :-
 1183        (   augmenting_path([[S]], Levels, T) ->
 1184            phrase(augmenting_path(S, T), Path),
 1185            Path = [augment(_,First,_)|Rest],
 1186            path_minimum(Rest, First, Min),
 1187            % format("augmenting path: ~w\n", [Min]),
 1188            maplist(augment(Min), Path),
 1189            maplist(maplist(clear_parent), Levels),
 1190            maximum_flow(S, T)
 1191        ;   true
 1192        ).
 1193
 1194clear_parent(V) :- del_attr(V, parent).
 1195
 1196augmenting_path(Levels0, Levels, T) :-
 1197        Levels0 = [Vs|_],
 1198        Levels1 = [Tos|Levels0],
 1199        phrase(reachables(Vs), Tos),
 1200        Tos = [_|_],
 1201        (   member(To, Tos), To == T -> Levels = Levels1
 1202        ;   augmenting_path(Levels1, Levels, T)
 1203        ).
 1204
 1205reachables([])     --> [].
 1206reachables([V|Vs]) -->
 1207        { get_attr(V, edges, Es) },
 1208        reachables_(Es, V),
 1209        reachables(Vs).
 1210
 1211reachables_([], _)     --> [].
 1212reachables_([E|Es], V) -->
 1213        reachable(E, V),
 1214        reachables_(Es, V).
 1215
 1216reachable(arc_from(V,F,_), P) -->
 1217        (   { \+ get_attr(V, parent, _),
 1218              get_attr(F, flow, Flow),
 1219              Flow > 0 } ->
 1220            { put_attr(V, parent, P-augment(F,Flow,-)) },
 1221            [V]
 1222        ;   []
 1223        ).
 1224reachable(arc_to(V,F,C), P) -->
 1225        (   { \+ get_attr(V, parent, _),
 1226              get_attr(F, flow, Flow),
 1227              (   C == inf ; Flow < C )} ->
 1228            { ( C == inf -> Diff = inf
 1229              ;   Diff is C - Flow
 1230              ),
 1231              put_attr(V, parent, P-augment(F,Diff,+)) },
 1232            [V]
 1233        ;   []
 1234        ).
 1235
 1236
 1237path_minimum([], Min, Min).
 1238path_minimum([augment(_,A,_)|As], Min0, Min) :-
 1239        (   A == inf -> Min1 = Min0
 1240        ;   Min1 is min(Min0,A)
 1241        ),
 1242        path_minimum(As, Min1, Min).
 1243
 1244augment(Min, augment(F,_,Sign)) :-
 1245        get_attr(F, flow, Flow0),
 1246        flow_(Sign, Flow0, Min, Flow),
 1247        put_attr(F, flow, Flow).
 1248
 1249flow_(+, F0, A, F) :- F is F0 + A.
 1250flow_(-, F0, A, F) :- F is F0 - A.
 1251
 1252augmenting_path(S, V) -->
 1253        (   { V == S } -> []
 1254        ;   { get_attr(V, parent, V1-Augment) },
 1255            [Augment],
 1256            augmenting_path(S, V1)
 1257        ).
 1258
 1259%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 1260
 1261naive_init(Supplies, _, Costs, Alphas, Betas) :-
 1262        same_length(Supplies, Alphas),
 1263        maplist(=(0), Alphas),
 1264        lists_transpose(Costs, TCs),
 1265        maplist(min_list, TCs, Betas).
 1266
 1267%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 1268
 1269lists_transpose([], []).
 1270lists_transpose([L|Ls], Ts) :- foldl(transpose_, L, Ts, [L|Ls], _).
 1271
 1272transpose_(_, Fs, Lists0, Lists) :-
 1273        maplist(list_first_rest, Lists0, Fs, Lists).
 1274
 1275list_first_rest([L|Ls], L, Ls).
 1276
 1277%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 1278/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1279   TODO: use attributed variables throughout
 1280- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 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.
 1294transportation(Supplies, Demands, Costs, Transport) :-
 1295        length(Supplies, LAs),
 1296        length(Demands, LBs),
 1297        naive_init(Supplies, Demands, Costs, Alphas, Betas),
 1298        network_head(Supplies, 1, SArcs, []),
 1299        network_tail(Demands, 1, DArcs, []),
 1300        numlist(1, LAs, Sources),
 1301        numlist(1, LBs, Sinks0),
 1302        maplist(make_sink, Sinks0, Sinks),
 1303        append(SArcs, DArcs, Torso),
 1304        alpha_beta(Torso, Sources, Sinks, Costs, Alphas, Betas, Flow),
 1305        flow_transport(Supplies, 1, Demands, Flow, Transport).
 1306
 1307flow_transport([], _, _, _, []).
 1308flow_transport([_|Rest], N, Demands, Flow, [Line|Lines]) :-
 1309        transport_line(Demands, N, 1, Flow, Line),
 1310        N1 is N + 1,
 1311        flow_transport(Rest, N1, Demands, Flow, Lines).
 1312
 1313transport_line([], _, _, _, []).
 1314transport_line([_|Rest], I, J, Flow, [L|Ls]) :-
 1315        (   get_assoc(I, Flow, As), get_assoc(p(J), As, arc(I,p(J),F,_)) -> L = F
 1316        ;   L = 0
 1317        ),
 1318        J1 is J + 1,
 1319        transport_line(Rest, I, J1, Flow, Ls).
 1320
 1321
 1322make_sink(N, p(N)).
 1323
 1324network_head([], _) --> [].
 1325network_head([S|Ss], N) -->
 1326        [arc(s,N,S)],
 1327        { N1 is N + 1 },
 1328        network_head(Ss, N1).
 1329
 1330network_tail([], _) --> [].
 1331network_tail([D|Ds], N) -->
 1332        [arc(p(N),t,D)],
 1333        { N1 is N + 1 },
 1334        network_tail(Ds, N1).
 1335
 1336network_connections([], _, _, _) --> [].
 1337network_connections([A|As], Betas, [Cs|Css], N) -->
 1338        network_connections(Betas, Cs, A, N, 1),
 1339        { N1 is N + 1 },
 1340        network_connections(As, Betas, Css, N1).
 1341
 1342network_connections([], _, _, _, _) --> [].
 1343network_connections([B|Bs], [C|Cs], A, N, PN) -->
 1344        (   { C =:= A + B } -> [arc(N,p(PN),inf)]
 1345        ;   []
 1346        ),
 1347        { PN1 is PN + 1 },
 1348        network_connections(Bs, Cs, A, N, PN1).
 1349
 1350alpha_beta(Torso, Sources, Sinks, Costs, Alphas, Betas, Flow) :-
 1351        network_connections(Alphas, Betas, Costs, 1, Cons, []),
 1352        append(Torso, Cons, Arcs),
 1353        edmonds_karp(Arcs, MaxFlow),
 1354        mark_hashes(MaxFlow, MArcs, MRevArcs),
 1355        all_markable(MArcs, MRevArcs, Markable),
 1356        mark_unmark(Sources, Markable, MarkSources, UnmarkSources),
 1357        (   MarkSources == [] -> Flow = MaxFlow
 1358        ;   mark_unmark(Sinks, Markable, MarkSinks0, UnmarkSinks0),
 1359            maplist(un_p, MarkSinks0, MarkSinks),
 1360            maplist(un_p, UnmarkSinks0, UnmarkSinks),
 1361            MarkSources = [FirstSource|_],
 1362            UnmarkSinks = [FirstSink|_],
 1363            theta(FirstSource, FirstSink, Costs, Alphas, Betas, TInit),
 1364            theta(MarkSources, UnmarkSinks, Costs, Alphas, Betas, TInit, Theta),
 1365            duals_add(MarkSources, Alphas, Theta, Alphas1),
 1366            duals_add(UnmarkSinks, Betas, Theta, Betas1),
 1367            Theta1 is -Theta,
 1368            duals_add(UnmarkSources, Alphas1, Theta1, Alphas2),
 1369            duals_add(MarkSinks, Betas1, Theta1, Betas2),
 1370            alpha_beta(Torso, Sources, Sinks, Costs, Alphas2, Betas2, Flow)
 1371        ).
 1372
 1373mark_hashes(MaxFlow, Arcs, RevArcs) :-
 1374        assoc_to_list(MaxFlow, FlowList),
 1375        maplist(un_arc, FlowList, FlowList1),
 1376        flatten(FlowList1, FlowList2),
 1377        empty_assoc(E),
 1378        mark_arcs(FlowList2, E, Arcs),
 1379        mark_revarcs(FlowList2, E, RevArcs).
 1380
 1381un_arc(_-Ls0, Ls) :-
 1382        assoc_to_list(Ls0, Ls1),
 1383        maplist(un_arc_, Ls1, Ls).
 1384
 1385un_arc_(_-Ls, Ls).
 1386
 1387mark_arcs([], Arcs, Arcs).
 1388mark_arcs([arc(From,To,F,C)|Rest], Arcs0, Arcs) :-
 1389        (   get_assoc(From, Arcs0, As) -> true
 1390        ;   As = []
 1391        ),
 1392        (   C == inf -> As1 = [To|As]
 1393        ;   F < C -> As1 = [To|As]
 1394        ;   As1 = As
 1395        ),
 1396        put_assoc(From, Arcs0, As1, Arcs1),
 1397        mark_arcs(Rest, Arcs1, Arcs).
 1398
 1399mark_revarcs([], Arcs, Arcs).
 1400mark_revarcs([arc(From,To,F,_)|Rest], Arcs0, Arcs) :-
 1401        (   get_assoc(To, Arcs0, Fs) -> true
 1402        ;   Fs = []
 1403        ),
 1404        (   F > 0 -> Fs1 = [From|Fs]
 1405        ;   Fs1 = Fs
 1406        ),
 1407        put_assoc(To, Arcs0, Fs1, Arcs1),
 1408        mark_revarcs(Rest, Arcs1, Arcs).
 1409
 1410
 1411un_p(p(N), N).
 1412
 1413duals_add([], Alphas, _, Alphas).
 1414duals_add([S|Ss], Alphas0, Theta, Alphas) :-
 1415        add_to_nth(1, S, Alphas0, Theta, Alphas1),
 1416        duals_add(Ss, Alphas1, Theta, Alphas).
 1417
 1418add_to_nth(N, N, [A0|As], Theta, [A|As]) :- !,
 1419        A is A0 + Theta.
 1420add_to_nth(N0, N, [A|As0], Theta, [A|As]) :-
 1421        N1 is N0 + 1,
 1422        add_to_nth(N1, N, As0, Theta, As).
 1423
 1424
 1425theta(Source, Sink, Costs, Alphas, Betas, Theta) :-
 1426        nth1(Source, Costs, Row),
 1427        nth1(Sink, Row, C),
 1428        nth1(Source, Alphas, A),
 1429        nth1(Sink, Betas, B),
 1430        Theta is (C - A - B) rdiv 2.
 1431
 1432theta([], _, _, _, _, Theta, Theta).
 1433theta([Source|Sources], Sinks, Costs, Alphas, Betas, Theta0, Theta) :-
 1434        theta_(Sinks, Source, Costs, Alphas, Betas, Theta0, Theta1),
 1435        theta(Sources, Sinks, Costs, Alphas, Betas, Theta1, Theta).
 1436
 1437theta_([], _, _, _, _, Theta, Theta).
 1438theta_([Sink|Sinks], Source, Costs, Alphas, Betas, Theta0, Theta) :-
 1439        theta(Source, Sink, Costs, Alphas, Betas, Theta1),
 1440        Theta2 is min(Theta0, Theta1),
 1441        theta_(Sinks, Source, Costs, Alphas, Betas, Theta2, Theta).
 1442
 1443
 1444mark_unmark(Nodes, Hash, Mark, Unmark) :-
 1445        mark_unmark(Nodes, Hash, Mark, [], Unmark, []).
 1446
 1447mark_unmark([], _, Mark, Mark, Unmark, Unmark).
 1448mark_unmark([Node|Nodes], Markable, Mark0, Mark, Unmark0, Unmark) :-
 1449        (   memberchk(Node, Markable) ->
 1450            Mark0 = [Node|Mark1],
 1451            Unmark0 = Unmark1
 1452        ;   Mark0 = Mark1,
 1453            Unmark0 = [Node|Unmark1]
 1454        ),
 1455        mark_unmark(Nodes, Markable, Mark1, Mark, Unmark1, Unmark).
 1456
 1457all_markable(Flow, RevArcs, Markable) :-
 1458        phrase(markable(s, [], _, Flow, RevArcs), Markable).
 1459
 1460all_markable([], Visited, Visited, _, _) --> [].
 1461all_markable([To|Tos], Visited0, Visited, Arcs, RevArcs) -->
 1462        (   { memberchk(To, Visited0) } -> { Visited0 = Visited1 }
 1463        ;   markable(To, [To|Visited0], Visited1, Arcs, RevArcs)
 1464        ),
 1465        all_markable(Tos, Visited1, Visited, Arcs, RevArcs).
 1466
 1467markable(Current, Visited0, Visited, Arcs, RevArcs) -->
 1468        { (   Current = p(_) ->
 1469              (   get_assoc(Current, RevArcs, Fs) -> true
 1470              ;   Fs = []
 1471              )
 1472          ;   (   get_assoc(Current, Arcs, Fs) -> true
 1473              ;   Fs = []
 1474              )
 1475          ) },
 1476        [Current],
 1477        all_markable(Fs, [Current|Visited0], Visited, Arcs, RevArcs).
 1478
 1479/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1480   solve(File) -- read input from File.
 1481
 1482   Format (NS = number of sources, ND = number of demands):
 1483
 1484   NS
 1485   ND
 1486   S1 S2 S3 ...
 1487   D1 D2 D3 ...
 1488   C11 C12 C13 ...
 1489   C21 C22 C23 ...
 1490   ... ... ... ...
 1491- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1492
 1493input(Ss, Ds, Costs) -->
 1494        integer(NS),
 1495        integer(ND),
 1496        n_integers(NS, Ss),
 1497        n_integers(ND, Ds),
 1498        n_kvectors(NS, ND, Costs).
 1499
 1500n_kvectors(0, _, []) --> !.
 1501n_kvectors(N, K, [V|Vs]) -->
 1502        n_integers(K, V),
 1503        { N1 is N - 1 },
 1504        n_kvectors(N1, K, Vs).
 1505
 1506n_integers(0, []) --> !.
 1507n_integers(N, [I|Is]) --> integer(I), { N1 is N - 1 }, n_integers(N1, Is).
 1508
 1509
 1510number([D|Ds]) --> digit(D), number(Ds).
 1511number([D])    --> digit(D).
 1512
 1513digit(D) --> [D], { between(0'0, 0'9, D) }.
 1514
 1515integer(N) --> number(Ds), !, ws, { name(N, Ds) }.
 1516
 1517ws --> [W], { W =< 0' }, !, ws.  % closing quote for syntax highlighting: '
 1518ws --> [].
 1519
 1520solve(File) :-
 1521        time((phrase_from_file(input(Supplies, Demands, Costs), File),
 1522              transportation(Supplies, Demands, Costs, Matrix),
 1523              maplist(print_row, Matrix))),
 1524        halt.
 1525
 1526print_row(R) :- maplist(print_row_, R), nl.
 1527
 1528print_row_(N) :- format("~w ", [N]).
 1529
 1530
 1531% ?- 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).
 1532%@ Ms = [[0, 3, 9, 0], [0, 7, 0, 0], [3, 5, 0, 6]],
 1533%@ Vs = [].
 1534
 1535
 1536%?- call_residue_vars(simplex:solve('instance_80_80.txt'), Vs).
 1537
 1538%?- call_residue_vars(simplex:solve('instance_3_4.txt'), Vs).
 1539
 1540%?- call_residue_vars(simplex:solve('instance_100_100.txt'), Vs).
 1541
 1542
 1543%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 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.
 1554% Assignment problem - for now, reduce to transportation problem
 1555assignment(Costs, Assignment) :-
 1556        length(Costs, LC),
 1557        length(Supply, LC),
 1558        all_one(Supply),
 1559        transportation(Supply, Supply, Costs, Assignment).
 1560
 1561/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1562   Messages
 1563- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1564
 1565:- multifile prolog:message//1. 1566
 1567prolog:message(simplex(bounded)) -->
 1568        ['Using library(simplex) with bounded arithmetic may yield wrong results.'-[]].
 1569
 1570warn_if_bounded_arithmetic :-
 1571        (   current_prolog_flag(bounded, true) ->
 1572            print_message(warning, simplex(bounded))
 1573        ;   true
 1574        ).
 1575
 1576:- initialization(warn_if_bounded_arithmetic).