From c6264cf098a8200b922def8fb68e549e8782624c Mon Sep 17 00:00:00 2001 From: Markus Triska Date: Sun, 6 Mar 2022 21:36:55 +0100 Subject: [PATCH] ADDED: library(simplex), for solving linear programming problems. --- README.md | 3 + src/lib/simplex.pl | 1383 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 1386 insertions(+) create mode 100644 src/lib/simplex.pl diff --git a/README.md b/README.md index 7999dcb4..cb551e3a 100644 --- a/README.md +++ b/README.md @@ -563,6 +563,9 @@ The modules that ship with Scryer Prolog are also called * [`tls`](src/lib/tls.pl) Predicates for negotiating TLS connections explicitly. * [`ugraphs`](src/lib/ugraphs.pl) Graph manipulation library +* [`simplex`](src/lib/simplex.pl) Providing `assignment/2`, + `transportation/4` and other predicates for solving linear + programming problems. To use predicates provided by the `lists` library, write: diff --git a/src/lib/simplex.pl b/src/lib/simplex.pl new file mode 100644 index 00000000..9522dab6 --- /dev/null +++ b/src/lib/simplex.pl @@ -0,0 +1,1383 @@ +/* + Author: Markus Triska + E-mail: triska@metalevel.at + WWW: http://www.metalevel.at + Copyright (C): 2005-2022, Markus Triska + + Part of Scryer Prolog. All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the + distribution. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS + FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE + COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, + INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, + BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN + ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + POSSIBILITY OF SUCH DAMAGE. +*/ + + +:- module(simplex, + [ + assignment/2, + constraint/3, + constraint/4, + constraint_add/4, + gen_state/1, + maximize/3, + minimize/3, + objective/2, + shadow_price/3, + transportation/4, + variable_value/3 + ]). + +:- use_module(library(assoc)). +:- use_module(library(pio)). +:- use_module(library(lists)). +:- use_module(library(dcgs)). +:- use_module(library(charsio)). +:- use_module(library(format)). +:- use_module(library(between)). +:- use_module(library(atts)). +:- use_module(library(arithmetic)). + +/** library(simplex): Solve linear programming problems + +This library provides several predicates for solving linear +programming problems with the simplex algorithm, and also includes +efficient algorithms for transportation and assignment problems. + +Efficiency could be improved significantly by changing the +implementation to the revised simplex method, benefiting from sparse +matrices. + +If you are interested in cooperating on such improvements, please +contact me! Enhancing the performance of this library would be a great +thesis project, for example. + +## Introduction {#simplex-intro} + +A *linear programming problem* or simply *linear program* (LP) +consists of: + + - a set of _linear_ **constraints** + - a set of **variables** + - a _linear_ **objective function**. + +The goal is to assign values to the variables so as to _maximize_ (or +minimize) the value of the objective function while satisfying all +constraints. + +Many optimization problems can be modeled in this way. As one basic +example, consider a knapsack with fixed capacity C, and a number of +items with sizes `s(i)` and values `v(i)`. The goal is to put as many +items as possible in the knapsack (not exceeding its capacity) while +maximizing the sum of their values. + +As another example, suppose you are given a set of _coins_ with +certain values, and you are to find the minimum number of coins such +that their values sum up to a fixed amount. Instances of these +problems are solved below. + +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 {#simplex-ex-1} + +This is the "radiation therapy" example, taken from _Introduction to +Operations Research_ by Hillier and Lieberman. + +[**Prolog DCG notation**](https://www.metalevel.at/prolog/dcg) is +used to _implicitly_ thread the state through posting the constraints: + +== +:- use_module(library(simplex)). +:- use_module(library(dcgs)). + +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). + S = solved(...), Val1 = 15 rdiv 2, Val2 = 9 rdiv 2. +== + +## Example 2 {#simplex-ex-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). + S = solved(...), X1 = 1 rdiv 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 {#simplex-ex-3} + +We are given: + + - 3 coins each worth 1 unit + - 20 coins each worth 5 units and + - 10 coins each worth 20 units. + +The task is to find a _minimal_ number of these coins that amount to +111 units in total. We introduce variables `c(1)`, `c(5)` and `c(20)` +denoting how many coins to take of the respective type: + +== +:- use_module(library(simplex)). + +coins(S) :- + gen_state(S0), + coins(S0, S). + +coins --> + constraint([c(1), 5*c(5), 20*c(20)] = 111), + constraint([c(1)] =< 3), + constraint([c(5)] =< 20), + constraint([c(20)] =< 10), + constraint([c(1)] >= 0), + constraint([c(5)] >= 0), + constraint([c(20)] >= 0), + constraint(integral(c(1))), + constraint(integral(c(5))), + constraint(integral(c(20))), + minimize([c(1), c(5), c(20)]). +== + +An example query: + +== +?- coins(S), variable_value(S, c(1), C1), + variable_value(S, c(5), C5), + variable_value(S, c(20), C20). + S = solved(...), C1 = 1 rdiv 1, C5 = 2 rdiv 1, C20 = 5 rdiv 1. +== + +@author [Markus Triska](https://www.metalevel.at) +*/ + + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + General Simplex Algorithm + Structures used: + + tableau(Objective, Variables, Indicators, Constraints) + *) objective function, represented as row + *) list of variables corresponding to columns + *) indicators denoting which variables are still active + *) constraints as rows + + row(Var, Left, Right) + *) the basic variable corresponding to this row + *) coefficients of the left-hand side of the constraint + *) right-hand side of the constraint +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + +find_row(Variable, [Row|Rows], R) :- + Row = row(V, _, _), + ( V == Variable -> R = Row + ; find_row(Variable, Rows, R) + ). + +%% variable_value(+State, +Variable, -Value) +% +% Value is unified with the value obtained for Variable. State must +% correspond to a solved instance. + +variable_value(State, Variable, Value) :- + functor(State, F, _), + ( F == solved -> + solved_tableau(State, Tableau), + tableau_rows(Tableau, Rows), + ( find_row(Variable, Rows, Row) -> Row = row(_, _, Value) + ; Value = 0 + ) + ; F == clpr_solved -> no_clpr + ). + +no_clpr :- throw(clpr_not_supported). + +var_zero(State, _Coeff*Var) :- variable_value(State, Var, 0). + +list_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. + +shadow_price(State, Name, Value) :- + functor(State, F, _), + ( F == solved -> + solved_tableau(State, Tableau), + tableau_objective(Tableau, row(_,Left,_)), + tableau_variables(Tableau, Variables), + solved_names(State, Names), + memberchk(user(Name)-Var, Names), + list_first(Variables, Var, Nth0), + nth0(Nth0, Left, Value) + ; F == clpr_solved -> no_clpr + ). + +%% objective(+State, -Objective) +% +% Unifies Objective with the result of the objective function at the +% obtained extremum. State must correspond to a solved instance. + +objective(State, Obj) :- + functor(State, F, _), + ( F == solved -> + solved_tableau(State, Tableau), + tableau_objective(Tableau, Objective), + Objective = row(_, _, Obj) + ; no_clpr + ). + + % interface functions that access tableau components + +tableau_objective(tableau(Obj, _, _, _), Obj). +tableau_rows(tableau(_, _, _, Rows), Rows). +tableau_indicators(tableau(_, _, Inds, _), Inds). +tableau_variables(tableau(_, Vars, _, _), Vars). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + interface functions that access and modify state components + + state is a structure of the form + state(Num, Names, Cs, Is) + Num: used to obtain new unique names for slack variables in a side-effect + free way (increased by one and threaded through) + Names: list of Name-Var, correspondence between constraint-names and + names of slack/artificial variables to obtain shadow prices later + Cs: list of constraints + Is: list of integer variables + + constraints are initially represented as c(Name, Left, Op, Right), + and after normalizing as c(Var, Left, Right). Name of unnamed constraints + is 0. The distinction is important for merging constraints (mainly in + branch and bound) with existing ones. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + +constraint_name(c(Name, _, _, _), Name). +constraint_op(c(_, _, Op, _), Op). +constraint_left(c(_, Left, _, _), Left). +constraint_right(c(_, _, _, Right), Right). + +%% gen_state(-State) +% +% Generates an initial state corresponding to an empty linear program. + +gen_state(state(0,[],[],[])). + +state_add_constraint(C, S0, S) :- + ( constraint_name(C, 0), constraint_left(C, [_Coeff*_Var]) -> + state_merge_constraint(C, S0, S) + ; state_add_constraint_(C, S0, S) + ). + +state_add_constraint_(C, state(VID,Ns,Cs,Is), state(VID,Ns,[C|Cs],Is)). + +state_merge_constraint(C, S0, S) :- + constraint_left(C, [Coeff0*Var0]), + constraint_right(C, Right0), + constraint_op(C, Op), + ( Coeff0 =:= 0 -> + ( Op == (=) -> Right0 =:= 0, S0 = S + ; Op == (=<) -> S0 = S + ; Op == (>=) -> Right0 =:= 0, S0 = S + ) + ; Coeff0 < 0 -> state_add_constraint_(C, S0, S) + ; Right is Right0 rdiv Coeff0, + state_constraints(S0, Cs), + ( select(c(0, [1*Var0], Op, CRight), Cs, RestCs) -> + ( Op == (=) -> CRight =:= Right, S0 = S + ; Op == (=<) -> + NewRight is min(Right, CRight), + NewCs = [c(0, [1*Var0], Op, NewRight)|RestCs], + state_set_constraints(NewCs, S0, S) + ; Op == (>=) -> + NewRight is max(Right, CRight), + NewCs = [c(0, [1*Var0], Op, NewRight)|RestCs], + state_set_constraints(NewCs, S0, S) + ) + ; state_add_constraint_(c(0, [1*Var0], Op, Right), S0, S) + ) + ). + + +state_add_name(Name, Var), [state(VID,[Name-Var|Ns],Cs,Is)] --> + [state(VID,Ns,Cs,Is)]. + +state_add_integral(Var, state(VID,Ns,Cs,Is), state(VID,Ns,Cs,[Var|Is])). + +state_constraints(state(_, _, Cs, _), Cs). +state_names(state(_,Names,_,_), Names). +state_integrals(state(_,_,_,Is), Is). +state_set_constraints(Cs, state(VID,Ns,_,Is), state(VID,Ns,Cs,Is)). +state_set_integrals(Is, state(VID,Ns,Cs,_), state(VID,Ns,Cs,Is)). + + +state_next_var(VarID0), [state(VarID1,Names,Cs,Is)] --> + [state(VarID0,Names,Cs,Is)], + { VarID1 is VarID0 + 1 }. + +solved_tableau(solved(Tableau, _, _), Tableau). +solved_names(solved(_, Names,_), Names). +solved_integrals(solved(_,_,Is), Is). + +% User-named constraints are wrapped with user/1 to also allow "0" in +% 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. + + +constraint(C, S0, S) :- + functor(S0, F, _), + ( F == state -> + ( C = integral(Var) -> state_add_integral(Var, S0, S) + ; constraint_(0, C, S0, S) + ) + ; F == clpr_state -> no_clpr + ). + +%% constraint(+Name, +Constraint, +S0, -S) +% +% Like constraint/3, and attaches the name Name (an atom or compound +% term) to the new constraint. + +constraint(Name, C, S0, S) :- constraint_(user(Name), C, S0, S). + +constraint_(Name, C, S0, S) :- + functor(S0, F, _), + ( F == state -> + ( C = integral(Var) -> state_add_integral(Var, S0, S) + ; C =.. [Op, Left0, Right0], + coeff_one(Left0, Left), + Right0 >= 0, + number_to_rational(Right0, Right), + state_add_constraint(c(Name, Left, Op, Right), S0, S) + ) + ; F == clpr_state -> no_clpr + ). + +%% 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. + +constraint_add(Name, A, S0, S) :- + functor(S0, F, _), + ( F == state -> + state_constraints(S0, Cs), + add_left(Cs, user(Name), A, Cs1), + state_set_constraints(Cs1, S0, S) + ; F == clpr_state -> no_clpr + ). + + +add_left([c(Name,Left0,Op,Right)|Cs], V, A, [c(Name,Left,Op,Right)|Rest]) :- + ( Name == V -> append(A, Left0, Left), Rest = Cs + ; Left0 = Left, add_left(Cs, V, A, Rest) + ). + +branching_variable(State, Variable) :- + solved_integrals(State, Integrals), + member(Variable, Integrals), + variable_value(State, Variable, Value), + \+ integer(Value). + + +worth_investigating(ZStar0, _, _) :- var(ZStar0). +worth_investigating(ZStar0, AllInt, Z) :- + nonvar(ZStar0), + ( AllInt =:= 1 -> Z1 is floor(Z) + ; Z1 = Z + ), + Z1 > ZStar0. + + +branch_and_bound(Objective, Solved, AllInt, ZStar0, ZStar, S0, S, Found) :- + objective(Solved, Z), + ( worth_investigating(ZStar0, AllInt, Z) -> + ( branching_variable(Solved, BrVar) -> + variable_value(Solved, BrVar, Value), + Value1 is floor(Value), + Value2 is Value1 + 1, + constraint([BrVar] =< Value1, S0, SubProb1), + ( maximize_(Objective, SubProb1, SubSolved1) -> + Sub1Feasible = 1, + objective(SubSolved1, Obj1) + ; Sub1Feasible = 0 + ), + constraint([BrVar] >= Value2, S0, SubProb2), + ( maximize_(Objective, SubProb2, SubSolved2) -> + Sub2Feasible = 1, + objective(SubSolved2, Obj2) + ; Sub2Feasible = 0 + ), + ( Sub1Feasible =:= 1, Sub2Feasible =:= 1 -> + ( Obj1 >= Obj2 -> + First = SubProb1, + Second = SubProb2, + FirstSolved = SubSolved1, + SecondSolved = SubSolved2 + ; First = SubProb2, + Second = SubProb1, + FirstSolved = SubSolved2, + SecondSolved = SubSolved1 + ), + branch_and_bound(Objective, FirstSolved, AllInt, ZStar0, ZStar1, First, Solved1, Found1), + branch_and_bound(Objective, SecondSolved, AllInt, ZStar1, ZStar2, Second, Solved2, Found2) + ; Sub1Feasible =:= 1 -> + branch_and_bound(Objective, SubSolved1, AllInt, ZStar0, ZStar1, SubProb1, Solved1, Found1), + Found2 = 0 + ; Sub2Feasible =:= 1 -> + Found1 = 0, + branch_and_bound(Objective, SubSolved2, AllInt, ZStar0, ZStar2, SubProb2, Solved2, Found2) + ; Found1 = 0, Found2 = 0 + ), + ( Found1 =:= 1, Found2 =:= 1 -> S = Solved2, ZStar = ZStar2 + ; Found1 =:= 1 -> S = Solved1, ZStar = ZStar1 + ; Found2 =:= 1 -> S = Solved2, ZStar = ZStar2 + ; S = S0, ZStar = ZStar0 + ), + Found is max(Found1, Found2) + ; S = Solved, ZStar = Z, Found = 1 + ) + ; ZStar = ZStar0, S = S0, Found = 0 + ). + +%% 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. + +maximize(Z0, S0, S) :- + coeff_one(Z0, Z1), + functor(S0, F, _), + ( F == state -> maximize_mip(Z1, S0, S) + ; F == clpr_state -> no_clpr + ). + +maximize_mip(Z, S0, S) :- + maximize_(Z, S0, Solved), + state_integrals(S0, Is), + ( Is == [] -> S = Solved + ; % arrange it so that branch and bound branches on variables + % in the same order the integrality constraints were stated in + reverse(Is, Is1), + state_set_integrals(Is1, S0, S1), + ( all_integers(Z, Is1) -> AllInt = 1 + ; AllInt = 0 + ), + branch_and_bound(Z, Solved, AllInt, _, _, S1, S, 1) + ). + +all_integers([], _). +all_integers([Coeff*V|Rest], Is) :- + integer(Coeff), + memberchk(V, Is), + all_integers(Rest, Is). + +%% minimize(+Objective, +S0, -S) +% +% Analogous to maximize/3. + +minimize(Z0, S0, S) :- + coeff_one(Z0, Z1), + functor(S0, F, _), + ( F == state -> + maplist(linsum_negate, Z1, Z2), + maximize_mip(Z2, S0, S1), + solved_tableau(S1, tableau(Obj, Vars, Inds, Rows)), + solved_names(S1, Names), + Obj = row(z, Left0, Right0), + all_times(Left0, -1, Left), + Right is -Right0, + Obj1 = row(z, Left, Right), + state_integrals(S0, Is), + S = solved(tableau(Obj1, Vars, Inds, Rows), Names, Is) + ; F == clpr_state -> no_clpr + ). + +op_pendant(>=, =<). +op_pendant(=<, >=). + +constraints_collapse([]) --> []. +constraints_collapse([C|Cs]) --> + { C = c(Name, Left, Op, Right) }, + ( { Name == 0, Left = [1*Var], op_pendant(Op, P) } -> + { Pendant = c(0, [1*Var], P, Right) }, + ( { select(Pendant, Cs, Rest) } -> + [c(0, Left, (=), Right)], + { CsLeft = Rest } + ; [C], + { CsLeft = Cs } + ) + ; [C], + { CsLeft = Cs } + ), + constraints_collapse(CsLeft). + +% solve a (relaxed) LP in standard form + +maximize_(Z, S0, S) :- + state_constraints(S0, Cs0), + phrase(constraints_collapse(Cs0), Cs1), + phrase(constraints_normalize(Cs1, Cs, As0), [S0], [S1]), + flatten(As0, As1), + ( As1 == [] -> + make_tableau(Z, Cs, Tableau0), + simplex(Tableau0, Tableau), + state_names(S1, Names), + state_integrals(S1, Is), + S = solved(Tableau, Names, Is) + ; state_names(S1, Names), + state_integrals(S1, Is), + two_phase_simplex(Z, Cs, As1, Names, Is, S) + ). + +flatten(Lss, Ls) :- + phrase(seqq(Lss), Ls). + +make_tableau(Z, Cs, Tableau) :- + ZC = c(_, Z, _), + phrase(constraints_variables([ZC|Cs]), Variables0), + sort(Variables0, Variables), + constraints_rows(Cs, Variables, Rows), + linsum_row(Variables, Z, Objective1), + all_times(Objective1, -1, Obj), + length(Variables, LVs), + length(Ones, LVs), + all_one(Ones), + Tableau = tableau(row(z, Obj, 0), Variables, Ones, Rows). + +all_one(Ones) :- maplist(=(1), Ones). + +proper_form(Variables, Rows, _Coeff*A, Obj0, Obj) :- + ( find_row(A, Rows, PivotRow) -> + list_first(Variables, A, Col), + row_eliminate(Obj0, PivotRow, Col, Obj) + ; Obj = Obj0 + ). + + +two_phase_simplex(Z, Cs, As, Names, Is, S) :- + % phase 1: minimize sum of articifial variables + make_tableau(As, Cs, Tableau0), + Tableau0 = tableau(Obj0, Variables, Inds, Rows), + foldl(proper_form(Variables, Rows), As, Obj0, Obj), + simplex(tableau(Obj, Variables, Inds, Rows), Tableau1), + maplist(var_zero(solved(Tableau1, _, _)), As), + % phase 2: remove artificial variables and solve actual LP. + tableau_rows(Tableau1, Rows2), + eliminate_artificial(As, As, Variables, Rows2, Rows3), + list_nths(As, Variables, Nths0), + nths_to_zero(Nths0, Inds, Inds1), + linsum_row(Variables, Z, Objective), + all_times(Objective, -1, Objective1), + foldl(proper_form(Variables, Rows3), Z, row(z, Objective1, 0), ObjRow), + simplex(tableau(ObjRow, Variables, Inds1, Rows3), Tableau), + S = solved(Tableau, Names, Is). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + If artificial variables are still in the basis, replace them with + non-artificial variables if possible. If that is not possible, the + constraint is ignored because it is redundant. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +eliminate_artificial([], _, _, Rows, Rows). +eliminate_artificial([_Coeff*A|Rest], As, Variables, Rows0, Rows) :- + ( select(row(A, Left, 0), Rows0, Others) -> + ( nth0(Col, Left, Coeff), + Coeff =\= 0, + nth0(Col, Variables, Var), + \+ memberchk(_*Var, As) -> + row_divide(row(A, Left, 0), Coeff, Row), + gauss_elimination(Rows0, Row, Col, Rows1), + swap_basic(Rows1, A, Var, Rows2) + ; Rows2 = Others + ) + ; Rows2 = Rows0 + ), + eliminate_artificial(Rest, As, Variables, Rows2, Rows). + +nths_to_zero([], Inds, Inds). +nths_to_zero([Nth|Nths], Inds0, Inds) :- + nth_to_zero(Inds0, 0, Nth, Inds1), + nths_to_zero(Nths, Inds1, Inds). + +nth_to_zero([], _, _, []). +nth_to_zero([I|Is], Curr, Nth, [Z|Zs]) :- + ( Curr =:= Nth -> [Z|Zs] = [0|Is] + ; Z = I, + Next is Curr + 1, + nth_to_zero(Is, Next, Nth, Zs) + ). + + +list_nths([], _, []). +list_nths([_Coeff*A|As], Variables, [Nth|Nths]) :- + list_first(Variables, A, Nth), + list_nths(As, Variables, Nths). + + +linsum_negate(Coeff0*Var, Coeff*Var) :- Coeff is -Coeff0. + +linsum_row([], _, []). +linsum_row([V|Vs], Ls, [C|Cs]) :- + ( member(Coeff*V, Ls) -> C = Coeff + ; C = 0 + ), + linsum_row(Vs, Ls, Cs). + +constraints_rows([], _, []). +constraints_rows([C|Cs], Vars, [R|Rs]) :- + C = c(Var, Left0, Right), + linsum_row(Vars, Left0, Left), + R = row(Var, Left, Right), + constraints_rows(Cs, Vars, Rs). + +constraints_normalize([], [], []) --> []. +constraints_normalize([C0|Cs0], [C|Cs], [A|As]) --> + { constraint_op(C0, Op), + constraint_left(C0, Left), + constraint_right(C0, Right), + constraint_name(C0, Name), + Con =.. [Op, Left, Right] }, + constraint_normalize(Con, Name, C, A), + constraints_normalize(Cs0, Cs, As). + +constraint_normalize(As0 =< B0, Name, c(Slack, [1*Slack|As0], B0), []) --> + state_next_var(Slack), + state_add_name(Name, Slack). +constraint_normalize(As0 = B0, Name, c(AID, [1*AID|As0], B0), [-1*AID]) --> + state_next_var(AID), + state_add_name(Name, AID). +constraint_normalize(As0 >= B0, Name, c(AID, [-1*Slack,1*AID|As0], B0), [-1*AID]) --> + state_next_var(Slack), + state_next_var(AID), + state_add_name(Name, AID). + +coeff_one([], []). +coeff_one([L|Ls], [Coeff*Var|Rest]) :- + ( L = A*B -> + number_to_rational(A, Coeff), + Var = B + ; Coeff = 1, Var = L + ), + coeff_one(Ls, Rest). + + +tableau_optimal(Tableau) :- + tableau_objective(Tableau, Objective), + tableau_indicators(Tableau, Indicators), + Objective = row(_, Left, _), + all_nonnegative(Left, Indicators). + +all_nonnegative([], []). +all_nonnegative([Coeff|As], [I|Is]) :- + ( I =:= 0 -> true + ; Coeff >= 0 + ), + all_nonnegative(As, Is). + +pivot_column(Tableau, PCol) :- + tableau_objective(Tableau, row(_, Left, _)), + tableau_indicators(Tableau, Indicators), + first_negative(Left, Indicators, 0, Index0, Val, RestL, RestI), + Index1 is Index0 + 1, + pivot_column(RestL, RestI, Val, Index1, Index0, PCol). + +first_negative([L|Ls], [I|Is], Index0, N, Val, RestL, RestI) :- + Index1 is Index0 + 1, + ( I =:= 0 -> first_negative(Ls, Is, Index1, N, Val, RestL, RestI) + ; ( L < 0 -> N = Index0, Val = L, RestL = Ls, RestI = Is + ; first_negative(Ls, Is, Index1, N, Val, RestL, RestI) + ) + ). + + +pivot_column([], _, _, _, N, N). +pivot_column([L|Ls], [I|Is], Coeff0, Index0, N0, N) :- + ( I =:= 0 -> Coeff1 = Coeff0, N1 = N0 + ; ( L < Coeff0 -> Coeff1 = L, N1 = Index0 + ; Coeff1 = Coeff0, N1 = N0 + ) + ), + Index1 is Index0 + 1, + pivot_column(Ls, Is, Coeff1, Index1, N1, N). + + +pivot_row(Tableau, PCol, PRow) :- + tableau_rows(Tableau, Rows), + pivot_row(Rows, PCol, false, _, 0, 0, PRow). + +pivot_row([], _, Bounded, _, _, Row, Row) :- Bounded. +pivot_row([Row|Rows], PCol, Bounded0, Min0, Index0, PRow0, PRow) :- + Row = row(_Var, Left, B), + nth0(PCol, Left, Ae), + ( Ae > 0 -> + Bounded1 = true, + Bound is B rdiv Ae, + ( Bounded0 -> + ( Bound < Min0 -> Min1 = Bound, PRow1 = Index0 + ; Min1 = Min0, PRow1 = PRow0 + ) + ; Min1 = Bound, PRow1 = Index0 + ) + ; Bounded1 = Bounded0, Min1 = Min0, PRow1 = PRow0 + ), + Index1 is Index0 + 1, + pivot_row(Rows, PCol, Bounded1, Min1, Index1, PRow1, PRow). + + +row_divide(row(Var, Left0, Right0), Div, row(Var, Left, Right)) :- + all_divide(Left0, Div, Left), + Right is Right0 rdiv Div. + + +all_divide([], _, []). +all_divide([R|Rs], Div, [DR|DRs]) :- + DR is R rdiv Div, + all_divide(Rs, Div, DRs). + +gauss_elimination([], _, _, []). +gauss_elimination([Row0|Rows0], PivotRow, Col, [Row|Rows]) :- + PivotRow = row(PVar, _, _), + Row0 = row(Var, _, _), + ( PVar == Var -> Row = PivotRow + ; row_eliminate(Row0, PivotRow, Col, Row) + ), + gauss_elimination(Rows0, PivotRow, Col, Rows). + +row_eliminate(row(Var, Ls0, R0), row(_, PLs, PR), Col, row(Var, Ls, R)) :- + nth0(Col, Ls0, Coeff), + ( Coeff =:= 0 -> Ls = Ls0, R = R0 + ; MCoeff is -Coeff, + all_times_plus([PR|PLs], MCoeff, [R0|Ls0], [R|Ls]) + ). + +all_times_plus([], _, _, []). +all_times_plus([A|As], T, [B|Bs], [AT|ATs]) :- + AT is A * T + B, + all_times_plus(As, T, Bs, ATs). + +all_times([], _, []). +all_times([A|As], T, [AT|ATs]) :- + AT is A * T, + all_times(As, T, ATs). + +simplex(Tableau0, Tableau) :- + ( tableau_optimal(Tableau0) -> Tableau0 = Tableau + ; pivot_column(Tableau0, PCol), + pivot_row(Tableau0, PCol, PRow), + Tableau0 = tableau(Obj0,Variables,Inds,Matrix0), + nth0(PRow, Matrix0, Row0), + Row0 = row(Leaving, Left0, _Right0), + nth0(PCol, Left0, PivotElement), + row_divide(Row0, PivotElement, Row1), + gauss_elimination([Obj0|Matrix0], Row1, PCol, [Obj|Matrix1]), + nth0(PCol, Variables, Entering), + swap_basic(Matrix1, Leaving, Entering, Matrix), + simplex(tableau(Obj,Variables,Inds,Matrix), Tableau) + ). + +swap_basic([Row0|Rows0], Old, New, Matrix) :- + Row0 = row(Var, Left, Right), + ( Var == Old -> Matrix = [row(New, Left, Right)|Rows0] + ; Matrix = [Row0|Rest], + swap_basic(Rows0, Old, New, Rest) + ). + +constraints_variables([]) --> []. +constraints_variables([c(_,Left,_)|Cs]) --> + variables(Left), + constraints_variables(Cs). + +variables([]) --> []. +variables([_Coeff*Var|Rest]) --> [Var], variables(Rest). + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + A dual algorithm ("algorithm alpha-beta" in Papadimitriou and + Steiglitz) is used for transportation and assignment problems. The + arising max-flow problem is solved with Edmonds-Karp, itself a dual + algorithm. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + An attributed variable is introduced for each node. Attributes: + node: Original name of the node. + edges: arc_to(To,F,Capacity) (F has an attribute "flow") or + arc_from(From,F,Capacity) + parent: used in breadth-first search +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +:- attribute + node/1, + edges/1, + flow/1, + parent/1. + +arcs([], Assoc, Assoc). +arcs([arc(From0,To0,C)|As], Assoc0, Assoc) :- + ( get_assoc(From0, Assoc0, From) -> Assoc1 = Assoc0 + ; put_assoc(From0, Assoc0, From, Assoc1), + put_atts(From, node(From0)) + ), + ( get_atts(From, edges(Es)) -> true + ; Es = [] + ), + put_atts(F, flow(0)), + put_atts(From, edges([arc_to(To,F,C)|Es])), + ( get_assoc(To0, Assoc1, To) -> Assoc2 = Assoc1 + ; put_assoc(To0, Assoc1, To, Assoc2), + put_atts(To, node(To0)) + ), + ( get_atts(To, edges(Es1)) -> true + ; Es1 = [] + ), + put_atts(To, edges([arc_from(From,F,C)|Es1])), + arcs(As, Assoc2, Assoc). + + +edmonds_karp(Arcs0, Arcs) :- + empty_assoc(E), + arcs(Arcs0, E, Assoc), + get_assoc(s, Assoc, S), + get_assoc(t, Assoc, T), + maximum_flow(S, T), + % fetch attvars before deleting visited edges + term_attributed_variables(S, AttVars), + phrase(flow_to_arcs(S), Ls), + arcs_assoc(Ls, Arcs), + maplist(del_attrs, AttVars). + +del_attrs(V) :- + put_atts(V, [-node(_),-edges(_),-flow(_),-parent(_)]). + +flow_to_arcs(V) --> + ( { get_atts(V, edges(Es)) } -> + { put_atts(V, -edges(_)), + get_atts(V, node(Name)) }, + flow_to_arcs_(Es, Name) + ; [] + ). + +flow_to_arcs_([], _) --> []. +flow_to_arcs_([E|Es], Name) --> + edge_to_arc(E, Name), + flow_to_arcs_(Es, Name). + +edge_to_arc(arc_from(_,_,_), _) --> []. +edge_to_arc(arc_to(To,F,C), Name) --> + { get_atts(To, node(NTo)), + get_atts(F, flow(Flow)) }, + [arc(Name,NTo,Flow,C)], + flow_to_arcs(To). + +arcs_assoc(Arcs, Hash) :- + empty_assoc(E), + arcs_assoc(Arcs, E, Hash). + +arcs_assoc([], Hs, Hs). +arcs_assoc([arc(From,To,F,C)|Rest], Hs0, Hs) :- + ( get_assoc(From, Hs0, As) -> Hs1 = Hs0 + ; put_assoc(From, Hs0, [], Hs1), + empty_assoc(As) + ), + put_assoc(To, As, arc(From,To,F,C), As1), + put_assoc(From, Hs1, As1, Hs2), + arcs_assoc(Rest, Hs2, Hs). + + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Strategy: Breadth-first search until we find a free right vertex in + the value graph, then find an augmenting path in reverse. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +maximum_flow(S, T) :- + ( augmenting_path([[S]], Levels, T) -> + phrase(augmenting_path(S, T), Path), + Path = [augment(_,First,_)|Rest], + path_minimum(Rest, First, Min), + % format("augmenting path: ~w\n", [Min]), + maplist(augment(Min), Path), + maplist(maplist(clear_parent), Levels), + maximum_flow(S, T) + ; true + ). + +clear_parent(V) :- put_atts(V, -parent(_)). + +augmenting_path(Levels0, Levels, T) :- + Levels0 = [Vs|_], + Levels1 = [Tos|Levels0], + phrase(reachables(Vs), Tos), + Tos = [_|_], + ( member(To, Tos), To == T -> Levels = Levels1 + ; augmenting_path(Levels1, Levels, T) + ). + +reachables([]) --> []. +reachables([V|Vs]) --> + { get_atts(V, edges(Es)) }, + reachables_(Es, V), + reachables(Vs). + +reachables_([], _) --> []. +reachables_([E|Es], V) --> + reachable(E, V), + reachables_(Es, V). + +reachable(arc_from(V,F,_), P) --> + ( { \+ get_atts(V, parent(_)), + get_atts(F, flow(Flow)), + Flow > 0 } -> + { put_atts(V, parent(P-augment(F,Flow,-))) }, + [V] + ; [] + ). +reachable(arc_to(V,F,C), P) --> + ( { \+ get_atts(V, parent(_)), + get_atts(F, flow(Flow)), + ( C == inf ; Flow < C )} -> + { ( C == inf -> Diff = inf + ; Diff is C - Flow + ), + put_atts(V, parent(P-augment(F,Diff,+))) }, + [V] + ; [] + ). + + +path_minimum([], Min, Min). +path_minimum([augment(_,A,_)|As], Min0, Min) :- + ( A == inf -> Min1 = Min0 + ; Min1 is min(Min0,A) + ), + path_minimum(As, Min1, Min). + +augment(Min, augment(F,_,Sign)) :- + get_atts(F, flow(Flow0)), + flow_(Sign, Flow0, Min, Flow), + put_atts(F, flow(Flow)). + +flow_(+, F0, A, F) :- F is F0 + A. +flow_(-, F0, A, F) :- F is F0 - A. + +augmenting_path(S, V) --> + ( { V == S } -> [] + ; { get_atts(V, parent(V1-Augment)) }, + [Augment], + augmenting_path(S, V1) + ). + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +naive_init(Supplies, _, Costs, Alphas, Betas) :- + same_length(Supplies, Alphas), + maplist(=(0), Alphas), + transpose(Costs, TCs), + maplist(min_list, TCs, Betas). + +min_list([L|Ls], Min) :- + foldl(min_, Ls, L, Min). + +min_(E, M0, M) :- M is min(E,M0). + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + TODO: use attributed variables throughout +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + +%% 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_. + +transportation(Supplies, Demands, Costs, Transport) :- + length(Supplies, LAs), + length(Demands, LBs), + naive_init(Supplies, Demands, Costs, Alphas, Betas), + network_head(Supplies, 1, SArcs, []), + network_tail(Demands, 1, DArcs, []), + numlist(1, LAs, Sources), + numlist(1, LBs, Sinks0), + maplist(make_sink, Sinks0, Sinks), + append(SArcs, DArcs, Torso), + alpha_beta(Torso, Sources, Sinks, Costs, Alphas, Betas, Flow), + flow_transport(Supplies, 1, Demands, Flow, Transport). + +flow_transport([], _, _, _, []). +flow_transport([_|Rest], N, Demands, Flow, [Line|Lines]) :- + transport_line(Demands, N, 1, Flow, Line), + N1 is N + 1, + flow_transport(Rest, N1, Demands, Flow, Lines). + +transport_line([], _, _, _, []). +transport_line([_|Rest], I, J, Flow, [L|Ls]) :- + ( get_assoc(I, Flow, As), get_assoc(p(J), As, arc(I,p(J),F,_)) -> L = F + ; L = 0 + ), + J1 is J + 1, + transport_line(Rest, I, J1, Flow, Ls). + + +make_sink(N, p(N)). + +network_head([], _) --> []. +network_head([S|Ss], N) --> + [arc(s,N,S)], + { N1 is N + 1 }, + network_head(Ss, N1). + +network_tail([], _) --> []. +network_tail([D|Ds], N) --> + [arc(p(N),t,D)], + { N1 is N + 1 }, + network_tail(Ds, N1). + +network_connections([], _, _, _) --> []. +network_connections([A|As], Betas, [Cs|Css], N) --> + network_connections(Betas, Cs, A, N, 1), + { N1 is N + 1 }, + network_connections(As, Betas, Css, N1). + +network_connections([], _, _, _, _) --> []. +network_connections([B|Bs], [C|Cs], A, N, PN) --> + ( { C =:= A + B } -> [arc(N,p(PN),inf)] + ; [] + ), + { PN1 is PN + 1 }, + network_connections(Bs, Cs, A, N, PN1). + +alpha_beta(Torso, Sources, Sinks, Costs, Alphas, Betas, Flow) :- + network_connections(Alphas, Betas, Costs, 1, Cons, []), + append(Torso, Cons, Arcs), + edmonds_karp(Arcs, MaxFlow), + mark_hashes(MaxFlow, MArcs, MRevArcs), + all_markable(MArcs, MRevArcs, Markable), + mark_unmark(Sources, Markable, MarkSources, UnmarkSources), + ( MarkSources == [] -> Flow = MaxFlow + ; mark_unmark(Sinks, Markable, MarkSinks0, UnmarkSinks0), + maplist(un_p, MarkSinks0, MarkSinks), + maplist(un_p, UnmarkSinks0, UnmarkSinks), + MarkSources = [FirstSource|_], + UnmarkSinks = [FirstSink|_], + theta(FirstSource, FirstSink, Costs, Alphas, Betas, TInit), + theta(MarkSources, UnmarkSinks, Costs, Alphas, Betas, TInit, Theta), + duals_add(MarkSources, Alphas, Theta, Alphas1), + duals_add(UnmarkSinks, Betas, Theta, Betas1), + Theta1 is -Theta, + duals_add(UnmarkSources, Alphas1, Theta1, Alphas2), + duals_add(MarkSinks, Betas1, Theta1, Betas2), + alpha_beta(Torso, Sources, Sinks, Costs, Alphas2, Betas2, Flow) + ). + +mark_hashes(MaxFlow, Arcs, RevArcs) :- + assoc_to_list(MaxFlow, FlowList), + maplist(un_arc, FlowList, FlowList1), + flatten(FlowList1, FlowList2), + empty_assoc(E), + mark_arcs(FlowList2, E, Arcs), + mark_revarcs(FlowList2, E, RevArcs). + +un_arc(_-Ls0, Ls) :- + assoc_to_list(Ls0, Ls1), + maplist(un_arc_, Ls1, Ls). + +un_arc_(_-Ls, Ls). + +mark_arcs([], Arcs, Arcs). +mark_arcs([arc(From,To,F,C)|Rest], Arcs0, Arcs) :- + ( get_assoc(From, Arcs0, As) -> true + ; As = [] + ), + ( C == inf -> As1 = [To|As] + ; F < C -> As1 = [To|As] + ; As1 = As + ), + put_assoc(From, Arcs0, As1, Arcs1), + mark_arcs(Rest, Arcs1, Arcs). + +mark_revarcs([], Arcs, Arcs). +mark_revarcs([arc(From,To,F,_)|Rest], Arcs0, Arcs) :- + ( get_assoc(To, Arcs0, Fs) -> true + ; Fs = [] + ), + ( F > 0 -> Fs1 = [From|Fs] + ; Fs1 = Fs + ), + put_assoc(To, Arcs0, Fs1, Arcs1), + mark_revarcs(Rest, Arcs1, Arcs). + + +un_p(p(N), N). + +duals_add([], Alphas, _, Alphas). +duals_add([S|Ss], Alphas0, Theta, Alphas) :- + add_to_nth(1, S, Alphas0, Theta, Alphas1), + duals_add(Ss, Alphas1, Theta, Alphas). + +add_to_nth(N, N, [A0|As], Theta, [A|As]) :- !, + A is A0 + Theta. +add_to_nth(N0, N, [A|As0], Theta, [A|As]) :- + N1 is N0 + 1, + add_to_nth(N1, N, As0, Theta, As). + + +theta(Source, Sink, Costs, Alphas, Betas, Theta) :- + nth1(Source, Costs, Row), + nth1(Sink, Row, C), + nth1(Source, Alphas, A), + nth1(Sink, Betas, B), + Theta is (C - A - B) rdiv 2. + +theta([], _, _, _, _, Theta, Theta). +theta([Source|Sources], Sinks, Costs, Alphas, Betas, Theta0, Theta) :- + theta_(Sinks, Source, Costs, Alphas, Betas, Theta0, Theta1), + theta(Sources, Sinks, Costs, Alphas, Betas, Theta1, Theta). + +theta_([], _, _, _, _, Theta, Theta). +theta_([Sink|Sinks], Source, Costs, Alphas, Betas, Theta0, Theta) :- + theta(Source, Sink, Costs, Alphas, Betas, Theta1), + Theta2 is min(Theta0, Theta1), + theta_(Sinks, Source, Costs, Alphas, Betas, Theta2, Theta). + + +mark_unmark(Nodes, Hash, Mark, Unmark) :- + mark_unmark(Nodes, Hash, Mark, [], Unmark, []). + +mark_unmark([], _, Mark, Mark, Unmark, Unmark). +mark_unmark([Node|Nodes], Markable, Mark0, Mark, Unmark0, Unmark) :- + ( memberchk(Node, Markable) -> + Mark0 = [Node|Mark1], + Unmark0 = Unmark1 + ; Mark0 = Mark1, + Unmark0 = [Node|Unmark1] + ), + mark_unmark(Nodes, Markable, Mark1, Mark, Unmark1, Unmark). + +all_markable(Flow, RevArcs, Markable) :- + phrase(markable(s, [], _, Flow, RevArcs), Markable). + +all_markable([], Visited, Visited, _, _) --> []. +all_markable([To|Tos], Visited0, Visited, Arcs, RevArcs) --> + ( { memberchk(To, Visited0) } -> { Visited0 = Visited1 } + ; markable(To, [To|Visited0], Visited1, Arcs, RevArcs) + ), + all_markable(Tos, Visited1, Visited, Arcs, RevArcs). + +markable(Current, Visited0, Visited, Arcs, RevArcs) --> + { ( Current = p(_) -> + ( get_assoc(Current, RevArcs, Fs) -> true + ; Fs = [] + ) + ; ( get_assoc(Current, Arcs, Fs) -> true + ; Fs = [] + ) + ) }, + [Current], + all_markable(Fs, [Current|Visited0], Visited, Arcs, RevArcs). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + solve(File) -- read input from File. + + Format (NS = number of sources, ND = number of demands): + + NS + ND + S1 S2 S3 ... + D1 D2 D3 ... + C11 C12 C13 ... + C21 C22 C23 ... + ... ... ... ... +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +input(Ss, Ds, Costs) --> + integer(NS), + integer(ND), + n_integers(NS, Ss), + n_integers(ND, Ds), + n_kvectors(NS, ND, Costs). + +n_kvectors(0, _, []) --> !. +n_kvectors(N, K, [V|Vs]) --> + n_integers(K, V), + { N1 is N - 1 }, + n_kvectors(N1, K, Vs). + +n_integers(0, []) --> !. +n_integers(N, [I|Is]) --> integer(I), { N1 is N - 1 }, n_integers(N1, Is). + + +number([D|Ds]) --> digit(D), number(Ds). +number([D]) --> digit(D). + +digit(D) --> [D], { char_type(D, decimal_digit) }. + +integer(N) --> number(Ds), !, ws, { number_chars(N, Ds) }. + +ws --> [W], { char_type(W, whitespace) }, !, ws. +ws --> []. + +solve(File) :- + time((phrase_from_file(input(Supplies, Demands, Costs), File), + transportation(Supplies, Demands, Costs, Matrix), + maplist(print_row, Matrix))), + halt. + +print_row(R) :- maplist(print_row_, R), nl. + +print_row_(N) :- format("~w ", [N]). + +nth1(N, Es, E) :- + N1 is N - 1, + nth0(N1, Es, E). + +%?- transportation([1,1], [1,1], [[1,1],[1,1]], Ms). + +%?- transportation([12,7,14], [3,15,9,6], [[20,50,10,60],[70,40,60,30],[40,80,70,40]], Ms). + +% ?- simplex: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). +%@ Ms = [[0,3,9,0],[0,7,0,0],[3,5,0,6]], Vs = []. + + +%?- call_residue_vars(simplex:solve('instance_80_80.txt'), Vs). + +%?- call_residue_vars(simplex:solve('instance_3_4.txt'), Vs). + +%?- call_residue_vars(simplex:solve('instance_100_100.txt'), Vs). + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% 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. + + +% Assignment problem - for now, reduce to transportation problem +assignment(Costs, Assignment) :- + length(Costs, LC), + length(Supply, LC), + all_one(Supply), + transportation(Supply, Supply, Costs, Assignment). + -- 2.54.0