for review: changes to linear programming module
Vanessa Joy TEAGUE
vjteag at students.cs.mu.OZ.AU
Fri Feb 27 21:33:42 AEDT 1998
Hi Tom,
Could you review this please?
I haven't compiled the mercury compiler with these
changes because the current termination analysis uses
floats and this version of lp uses rats.
Thanks,
Vanessa.
2c2
< % Copyright (C) 1994-1995, 1997 The University of Melbourne.
---
> % Copyright (C) 1997-1998 The University of Melbourne.
8a9
> % rewrite by vjteag, Feb 1998
10,26c11
< % This module implements a linear constraint solver that finds an
< % optimal solution to a set of linear [in]equalities with respect
< % to some objective function. It does this using the simplex method.
< %
< % The form of an [in]equation is
< % a1.x1 + a2.x2 + ... + an.xn {=<,=,>=} b
< % where alll the numbers are floats.
< % By defaut, there is an additional constraint on each of the `xn's:
< % xn >= 0
< % If you want xn to take on any value, you can include it in the list
< % of URS (UnRestricted in Sign) variables.
< %
< % The objective function is simply a weighted sum of the variables.
< %
< % The `x's are represented by `term__var's. The varset from which
< % they are allocated is passed to the solver because it needs to
< % introduce new variables as part of the solving algorithm.
---
> % This module implements:
27a13,54
> % A linear constraint solver that finds an optimal solution to a set of
> % linear [in]equalities with respect to some objective function. It does
> % this using the simplex method.
> %
> % The form of an [in]equation is
> % a1.x1 + a2.x2 + ... + an.xn {=<,=,>=} b
> % where all the numbers are rats, a variable xn may occur multiple
> % times in an equation (or the objective function) - the solver simplifies
> % all the inequations.
> % By defaut, there is an additional constraint on each of the `xn's:
> % xn >= 0
> % If you want xn to take on any value, you can include it in the list
> % of URS (UnRestricted in Sign) variables.
> %
> % The objective function is simply a weighted sum of the variables.
> %
> % The `x's are represented by `term__var's. The varset from which
> % they are allocated is passed to the solver because it needs to
> % introduce new variables as part of the solving algorithm.
> %
> %
> % Fourier-Motzkin elimination, which takes a polyhedron represented by
> % a list of inequations (represented in the same way as in the
> % constraint solver) and projects it onto a subset of the variables
> % contained in the inequations.
> %
> %
> % A convex hull algorithm, which takes a list of polyhedrons (represented
> % by lists of inequations), and computes the smallest convex set which is
> % a superset of all of them.
> % For details of the algorithm, see Benoy and King.
> %
> %
> % A test for entailment of one constraint (a single inequality) by a
> % conjunction of other inequalities (a polyhedron). This uses the constraint
> % solver to find the point on the polyhedron closest to the boundary of the
> % constraint, then checks whether it is on the correct side of the boundary.
> % It assumes that all variables are non-negative.
> %
> % A widening algorithm, which takes two polyhedra and retains only the
> % inequalities in the first polyhedron which are entailed by the second.
> %
35c62
< :- import_module float, io, list, map, std_util, term, varset.
---
> :- import_module rat, list, map, std_util, term, varset.
37c64
< :- type coeff == pair(var, float).
---
> :- type coeff == pair(var, rat).
40c67
< ---> eqn(list(coeff), operator, float).
---
> ---> eqn(list(coeff), operator, rat).
53,54c80,82
< ---> unsatisfiable
< ; satisfiable(float, map(var, float))
---
> ---> unbounded
> ; inconsistent
> ; satisfiable(rat, map(var, rat))
56a85,87
> :- type vector ---> pair(map(var, rat), rat).
>
> :- type matrix == list(vector).
60c91
< % Result, IO0, IO)
---
> % Result)
66,70c97,100
< % lp_solve binds Result either to `unsatisfiable' if the there
< % was no optimum value of the objective function (ie the
< % constraints were inconsistent, or the objective function
< % is unbounded by the constraints), or `satisfiable(ObjVal,
< % MapFromObjVarsToVals)'.
---
> % lp_solve binds Result to `unbounded' if the the objective
> % function is not bounded by the constraints, to 'inconsistent'
> % if the given constraints are inconsistent, or to
> % `satisfiable(ObjVal, MapFromObjVarsToVals)'.
73,74c103,104
< lp__result, io__state, io__state).
< :- mode lp_solve(in, in, in, in, in, out, di, uo) is det.
---
> lp__result).
> :- mode lp_solve(in, in, in, in, in, out) is det.
75a106,126
> :- pred project(equations, list(var), equations).
> :- mode project(in, in, out) is det.
>
> :- pred convex_hull(list(equations), equations, varset, varset).
> :- mode convex_hull(in, out, in, out) is det.
>
> :- pred widen(equations, equations, varset, equations).
> :- mode widen(in, in, in, out) is det.
>
> :- pred entailed(equation, equations, varset).
> :- mode entailed(in, in, in) is semidet.
>
> % Exported for testing purposes.
> :- pred eqns_to_vectors(equations, list(vector)).
> :- mode eqns_to_vectors(in, out) is det.
>
> :- pred is_false(equation).
> :- mode is_false(in) is semidet.
>
> :- pred false_eqn(equation).
> :- mode false_eqn(out) is det.
81c132
< :- import_module bool, int, require, set.
---
> :- import_module bool, int, io, require, set, string, array.
86c137,141
< map(var, pair(var)), % variables with urs.
---
> map(var, pair(var)), % map from variables with URS to
> % the corresponding pair of variables
> % that represent that variable in
> % the standard form (x = x' - x'',
> % x', x'' >= 0).
91,93c146,158
< lp_solve(Eqns, Dir, Obj, Varset0, URSVars, Result, IO0, IO) :-
< lp_info_init(Varset0, URSVars, Info0),
< lp_solve2(Eqns, Dir, Obj, Result, IO0, IO, Info0, _).
---
> is_false(eqn([], (>=), Const)) :-
> Const > zero.
> is_false(eqn([], (=<), Const)) :-
> Const < zero.
> is_false(eqn([], (=), Const)) :-
> not (Const = zero).
>
> false_eqn(eqn([], (=<), rat(-1,1))).
>
>
> lp_solve(Eqns, Dir, Obj, Varset, URSVars, Result) :-
> lp_info_init(Varset, URSVars, Info0),
> lp_solve2(Eqns, Dir, Obj, Result, Info0, _).
95,97c160,165
< :- pred lp_solve2(equations, direction, objective, lp__result,
< io__state, io__state, lp_info, lp_info).
< :- mode lp_solve2(in, in, in, out, di, uo, in, out) is det.
---
> %
> % lp_solve2(Eqns, Dir, Obj, Res, LPInfo0, LPInfo) takes a list
> % of inequations `Eqns', a direction for optimization `Dir', an objective
> % function `Obj', an I/O state `IO0' and an lp_info structure `LPInfo0'.
> % See inline comments for details on the algorithm.
> %
99c167,175
< lp_solve2(Eqns0, Dir, Obj0, Result, IO00, IO, Info0, Info) :-
---
> :- pred lp_solve2(equations, direction, objective, lp__result,
> lp_info, lp_info).
> :- mode lp_solve2(in, in, in, out, in, out) is det.
> lp_solve2(Eqns0, Dir, Obj0, Result, Info0, Info) :-
> % simplify the inequations and convert them
> % to standard form by introducing slack/excess/
> % artificial variables. We also expand URS variables
> % by replacing them with the difference of two
> % fresh variables.
100a177,179
>
> % If we're maximizing the objective function then we need
> % to negate all the coefficients in the objective.w
103c182
< negate_equation(eqn(Obj0, (=), 0.0), eqn(Obj1, _, _))
---
> negate_equation(eqn(Obj0, (=), zero), eqn(Obj1, _, _))
108a188
>
117,119c197,198
< get_art_vars(ArtVars, Info1, Info2),
< get_urs_vars(URSVars, Info2, Info),
< init_tableau(Rows, Cols, VarNums, URSVars, Tableau0),
---
> get_art_vars(ArtVars, Info1, Info),
> init_tableau(Rows, Cols, VarNums, URS, Tableau0),
122,123d200
< io__write_string("first tableau:\n", IO00, IO00a),
< show_tableau(Tableau1, IO00a, IO0),
128,130c205,206
< two_phase(Obj0, Obj, Dir, ArtVars, VarNums, Tableau1, Result,
< IO0, IO)
< ;
---
> two_phase(Obj0, Obj, ArtVars, VarNums, Tableau1, Result0)
> ;
132,148c208,209
< one_phase(Obj0, Obj, Dir, VarNums, Tableau1, Result, IO0, IO)
< ).
<
< %------------------------------------------------------------------------------%
<
< :- pred one_phase(list(coeff), list(coeff), direction, map(var, int), tableau,
< lp__result, io__state, io__state).
< :- mode one_phase(in, in, in, in, in, out, di, uo) is det.
<
< one_phase(Obj0, Obj, Dir, VarNums, Tableau0, Result, IO0, IO) :-
< insert_coeffs(Obj, 0, VarNums, Tableau0, Tableau1),
< GetObjVar = lambda([V::out] is nondet, (
< list__member(X, Obj0),
< X = V - _Cof
< )),
< solutions(GetObjVar, ObjVars),
< optimize(ObjVars, Tableau1, _, Result0, IO0, IO),
---
> one_phase(Obj0, Obj, VarNums, Tableau1, Result0)
> ),
155c216,219
< Result0 = unsatisfiable,
---
> Result0 = unbounded,
> Result = Result0
> ;
> Result0 = inconsistent,
166,168c230,245
< :- pred two_phase(list(coeff), list(coeff), direction, list(var), map(var, int),
< tableau, lp__result, io__state, io__state).
< :- mode two_phase(in, in, in, in, in, in, out, di, uo) is det.
---
> :- pred one_phase(list(coeff), list(coeff), map(var, int), tableau, lp__result).
> :- mode one_phase(in, in, in, in, out) is det.
> one_phase(Obj0, Obj, VarNums, Tableau0, Result) :-
> insert_coeffs(Obj, 0, VarNums, Tableau0, Tableau1),
> GetObjVar = lambda([V::out] is nondet, (
> list__member(X, Obj0),
> X = V - _Cof
> )),
> solutions(GetObjVar, ObjVars),
> optimize(ObjVars, Tableau1, _, Result).
> %------------------------------------------------------------------------------%
>
> :- pred two_phase(list(coeff), list(coeff), list(var), map(var, int),
> tableau, lp__result).
> :- mode two_phase(in, in, in, in, in, out) is det.
> two_phase(Obj0, Obj, ArtVars, VarNums, Tableau0, Result) :-
170d246
< two_phase(Obj0, Obj, Dir, ArtVars, VarNums, Tableau0, Result, IO0, IO) :-
176,183c252,259
< io__write_string("art tableau:\n", IO0, IO0a),
< show_tableau(Tableau1b, IO0a, IO0b),
< optimize(ArtVars, Tableau1b, Tableau1c,
< Res0, IO0b, IO1),
< (
< Res0 = unsatisfiable,
< Result = unsatisfiable,
< IO = IO1
---
> optimize(ArtVars, Tableau1b, Tableau1c, Res0),
> (
> Res0 = unbounded,
> Result = unbounded
>
> ;
> Res0 = inconsistent,
> Result = inconsistent
186,188c262,263
< ( Val \= 0.0 ->
< Result = unsatisfiable,
< IO = IO1
---
> ( Val \= zero ->
> Result = inconsistent
205,219c280
< optimize(ObjVars, Tableau4, _, Res1, IO1, IO),
< (
< Dir = max,
< Result = Res1
< ;
< Dir = min,
< (
< Res1 = unsatisfiable,
< Result = Res1
< ;
< Res1 = satisfiable(NOptVal, OptCoffs),
< OptVal is -NOptVal,
< Result = satisfiable(OptVal, OptCoffs)
< )
< )
---
> optimize(ObjVars, Tableau4, _, Result)
229c290
< construct_art_objective([V|Vs], [V - (1.0)|Rest]) :-
---
> construct_art_objective([V|Vs], [V - (one)|Rest]) :-
239a301,307
> % standardize_equation peforms the following operations on an
> % equation:
> % - ensures the constant is >= 0 (multiplying by -1 if
> % necessary)
> % - introduces slack, excess and artificial variables
> % - replace the URS variables with their corresponding
> % difference pair
245c313
< ( { Const0 < 0.0 } ->
---
> ( { Const0 < zero } ->
250c318
< { Coeffs = [Var - 1.0|Coeffs0] },
---
> { Coeffs = [Var - one|Coeffs0] },
258c326
< ( { Const0 < 0.0 } ->
---
> ( { Const0 < zero } ->
263c331
< { Coeffs = [Var - 1.0|Coeffs0] },
---
> { Coeffs = [Var - one|Coeffs0] },
271c339
< ( { Const0 < 0.0 } ->
---
> ( { Const0 < zero } ->
277c345
< { Coeffs = [SVar - (-1.0), AVar - (1.0)|Coeffs0] },
---
> { Coeffs = [SVar - (-one), AVar - (one)|Coeffs0] },
320c388
< :- pred add_var(map(var, float), var, float, map(var, float)).
---
> :- pred add_var(map(var, rat), var, rat, map(var, rat)).
327c395
< Acc1 = 0.0
---
> Acc1 = zero
411,417c479,482
< :- pred optimize(list(var), tableau, tableau, lp__result,
< io__state, io__state).
< :- mode optimize(in, in, out, out, di, uo) is det.
<
< optimize(ObjVars, A0, A, Result) -->
< io__write_string("initial tableau:\n"),
< show_tableau(A0),
---
> :- pred optimize(list(var), tableau, tableau, lp__result).
> :- mode optimize(in, in, out, out) is det.
>
> optimize(ObjVars, A0, A, Result) :-
419,420d483
< io__write_string("final tableau:\n"),
< show_tableau(A),
422,423c485,486
< { Res0 = no },
< { Result = unsatisfiable }
---
> Res0 = no ,
> Result = unbounded
425,429c488,492
< { Res0 = yes },
< { rhs_col(A, M) },
< { index(A, 0, M, ObjVal) },
< { extract_objective(ObjVars, A, ObjMap) },
< { Result = satisfiable(ObjVal, ObjMap) }
---
> Res0 = yes,
> rhs_col(A, M),
> index(A, 0, M, ObjVal),
> extract_objective(ObjVars, A, ObjMap),
> Result = satisfiable(ObjVal, ObjMap)
432c495
< :- pred extract_objective(list(var), tableau, map(var, float)).
---
> :- pred extract_objective(list(var), tableau, map(var, rat)).
439c502
< :- pred extract_obj_var(tableau, var, map(var, float), map(var, float)).
---
> :- pred extract_obj_var(tableau, var, map(var, rat), map(var, rat)).
453c516
< :- pred extract_obj_var2(tableau, var, float).
---
> :- pred extract_obj_var2(tableau, var, rat).
460c523
< index(Tab, Row, Col, 1.0),
---
> index(Tab, Row, Col, one),
468c531
< Val = 0.0
---
> Val = zero
471,472c534,535
< :- pred simplex(tableau, tableau, bool, io__state, io__state).
< :- mode simplex(in, out, out, di, uo) is det.
---
> :- pred simplex(tableau, tableau, bool).
> :- mode simplex(in, out, out) is det.
474c537
< simplex(A0, A, Result, IO0, IO) :-
---
> simplex(A0, A, Result) :-
480c543
< ( MinVal < 0.0 ->
---
> ( MinVal < zero ->
499d561
< IO = IO0,
508c570
< ( MaxVal > 0.0 ->
---
> ( MaxVal > zero ->
510a573,574
> require(nz(MaxVal), "simplex: zero
> divisor MaxVal"),
521,524c585,586
< MaxVal1 is MVal/CellVal,
< ( CellVal > 0.0, MaxVal1 =< MaxVal0 ->
< Max = yes(Row - MaxVal1)
< ;
---
> ( CellVal =< zero ->
> % CellVal = 0 indicates multiple optimal sol'ns
525a588,596
> ;
> require(nz(CellVal), "simplex: zero
> divisor CellVal"),
> MaxVal1 is MVal/CellVal,
> ( MaxVal1 =< MaxVal0 ->
> Max = yes(Row - MaxVal1)
> ;
> Max = Max0
> )
533d603
< IO = IO0,
537,540d606
< index(A0, P, Q, Apq),
< string__format("pivot on (%d, %d) = %f\n",
< [i(P), i(Q), f(Apq)], Str),
< io__write_string(Str, IO0, IO1),
542,543c608
< show_tableau(A1, IO1, IO2),
< simplex(A1, A, Result, IO2, IO)
---
> simplex(A1, A, Result)
556c621
< ( Val = 0.0 ->
---
> ( Val = zero ->
562c627
< ValF0 \= 0.0,
---
> ValF0 \= zero,
567a633,634
> require(nz(Fac0), "ensure_zero_obj_coeffs:
> zero divisor Fac0"),
585c652
< ( Val = 0.0 ->
---
> ( Val = zero ->
593c660
< Res = [1.0 - Row]
---
> Res = [one - Row]
599c666
< Zz \= 0.0
---
> Zz \= zero
637a705
> require(nz(Apq), "ScaleCell: zero divisor"),
648c716
< set_index(T0, J, K, 0.0, T)
---
> set_index(T0, J, K, zero, T)
653a722
> require(nz(Apq), "ScaleRow: zero divisor"),
658c727
< set_index(A3, P, Q, 1.0, A).
---
> set_index(A3, P, Q, one, A).
660c729
< :- pred row_op(float, int, int, tableau, tableau).
---
> :- pred row_op(rat, int, int, tableau, tableau).
683c752
< map(pair(int), float)
---
> map(pair(int), rat)
693c762
< :- pred index(tableau, int, int, float).
---
> :- pred index(tableau, int, int, rat).
712c781
< R = 0.0
---
> R = zero
715c784
< :- pred set_index(tableau, int, int, float, tableau).
---
> :- pred set_index(tableau, int, int, rat, tableau).
729c798
< ( R = 0.0 ->
---
> ( R = zero ->
810c879
< Z \= 0.0,
---
> Z \= zero,
814c883
< Solns = [_ - 1.0]
---
> Solns = [_ - one]
826c895
< :- import_module string.
---
> % For debugging ....
831,834d899
< /*
< show_tableau(_N, _M, _Tableau) --> [].
< */
<
853,854c918,920
< { string__format("%2.2f\t", [f(Val)], Str) },
< io__write_string(Str).
---
> { string__format("%2.2f\t", [f(float(Val))], Str) },
> io__write_string(Str),
> io__write(Val).
950a1017,1496
>
>
> % Fourier-Motzkin elimination.
>
> % project(Equations0, Variables, Equations) takes a list of inequalities
> % and eliminates the variables in the input list, i.e. it projects onto
> % those variables not in the list.
> project(Eqns, [], Eqns).
> project(Eqns0, Vars, Eqns) :-
> Vars = [_|_],
> fm_standardize_equations(Eqns0, Eqns1),
> eqns_to_vectors(Eqns1, Vecs),
> eliminate_vars(Vars, Vecs, Matrix),
> vectors_to_eqns(Matrix, Eqns).
>
>
> % fm_standardize_equations ensures all inequations are of the form
> % a1x1 + ... + anxn =< C. Inequalities containing >= are multiplied
> % through by -1, and equations are converted into two inequations
> % of the form a1x1 + ... + anxn =< C, -a1x1 - ... - anxn =< -C.
> :- pred fm_standardize_equations(equations, equations).
> :- mode fm_standardize_equations(in, out) is det.
>
> fm_standardize_equations(Eqns0, Eqns1) :-
> fm_standardize_equations_acc(Eqns0, [], Eqns1).
>
> :- pred fm_standardize_equations_acc(equations, equations, equations).
> :- mode fm_standardize_equations_acc(in, in, out) is det.
>
> fm_standardize_equations_acc([], Acc, Acc).
> fm_standardize_equations_acc([Eqn0|Eqns], Acc0, Acc) :-
> (
> Eqn0 = eqn(_, (=<), _),
> Acc1 = [Eqn0|Acc0]
> ;
> Eqn0 = eqn(Coeffs, (=), Const),
> Acc1 = [eqn(Coeffs, (=<), Const), Neg_eqn | Acc0],
> negate_equation(eqn(Coeffs, (>=), Const), Neg_eqn)
> ;
> Eqn0 = eqn(Coeffs, (>=), Const),
> Acc1 = [Neg_eqn | Acc0],
> negate_equation(eqn(Coeffs, (>=), Const), Neg_eqn)
> ),
> fm_standardize_equations_acc(Eqns, Acc1, Acc).
>
>
> eqns_to_vectors(Eqns, Vecs) :-
> Eqn_to_vec = lambda([eqn(List, Op, Num)::in,
> pair(Map, Num)::out] is det, (
> ( Op = (=<) ->
> map__init(Map0),
> coeff_list_to_map(List, Map0, Map)
> ;
> error("Equality or >= passed to eqns_to_vectors\n")
> )
> )),
> list__map(Eqn_to_vec, Eqns, Vecs).
>
>
> :- pred vectors_to_eqns(list(vector), equations).
> :- mode vectors_to_eqns(in, out) is det.
> vectors_to_eqns(Vectors, Equations) :-
> Vec_to_eqn = lambda([pair(Map,Rat)::in, Eqn::out] is det, (
> map__to_assoc_list(Map, List),
> Eqn = eqn(List, (=<), Rat)
> )),
> list__map(Vec_to_eqn, Vectors, Equations).
>
>
> % Ensures that there are no zeros in the resulting map.
> :- pred coeff_list_to_map(list(coeff), map(var, rat), map(var, rat)).
> :- mode coeff_list_to_map(in, in, out) is det.
> coeff_list_to_map([], Map, Map).
> coeff_list_to_map([Var-Num | Coeffs], Map0, Map) :-
> ( Num = zero ->
> Map1 = Map0
> ;
> ( map__search(Map0, Var, Num1) ->
> ( Num1 + Num = zero ->
> map__delete(Map0, Var, Map1)
> ;
> map__det_update(Map0, Var, Num+Num1, Map1)
> )
> ;
> map__det_insert(Map0, Var, Num, Map1)
> )
> ),
> coeff_list_to_map(Coeffs, Map1, Map).
>
>
> % normalize_vector(Vec0, Var, Vec1) divides all coefficients in
> % Vec0 by the absolute value of the coefficient of Var.
> % It is an error if the map contains a zero coefficient.
> :- pred normalize_vector(vector, var, vector).
> :- mode normalize_vector(in, in, out) is det.
> normalize_vector(pair(Map0, Num0), Var0, pair(Map, Num)) :-
> ( map__search(Map0, Var0, Coeff) ->
> Is_map_key = lambda([Var::out] is nondet, (
> map__member(Map0, Var, _)
> )),
> Divide_map_val = lambda([Var::in, Map1::in, Map2::out] is det, (
> map__lookup(Map1, Var, Coeff0),
> ( Coeff0 = zero ->
> error("Divide by zero in normalize_vector")
> ;
> require(nz(Coeff),
> "normalize_vector: zero divisor"),
> Coeff1 = Coeff0 / abs(Coeff),
> map__det_update(Map1, Var, Coeff1, Map2)
> )
> )),
> aggregate(Is_map_key, Divide_map_val, Map0, Map),
> require(nz(Coeff), "normalize_vector: zero divisor"),
> Num = Num0 / abs(Coeff)
> ;
> Map = Map0,
> Num = Num0
> ).
>
> % separate_vectors(List1, Variable, Posititves, Negatives, Zeroes):
> % Breaks a list of vectors up into three lists of vectors according to
> % whether the coefficient of Variable is positive, negative or zero.
> % Applies normalize_vector to each vector.
> :- pred separate_vectors(matrix, var, matrix, matrix, matrix).
> :- mode separate_vectors(in, in, out, out, out) is det.
> separate_vectors(Matrix, Var, Pos, Neg, Zero) :-
> separate_vectors_acc(Matrix, Var, [], [], [], Pos, Neg, Zero).
>
> :- pred separate_vectors_acc(matrix, var, matrix, matrix,
> matrix, matrix, matrix, matrix).
> :- mode separate_vectors_acc(in, in, in, in, in, out, out, out) is det.
>
> separate_vectors_acc([], _, Pos, Neg, Zeros, Pos, Neg, Zeros).
> separate_vectors_acc([Vec0|Vectors], Var, Pos0, Neg0, Zeros0, Pos, Neg, Zeros):-
> ( get_var_coeff(Vec0, Var, Coeff) ->
> normalize_vector(Vec0, Var, Vec1),
> ( Coeff > zero ->
> Pos1 = [Vec1 | Pos0],
> Neg1 = Neg0,
> Zeros1 = Zeros0
> ;
> Pos1 = Pos0,
> Neg1 = [Vec1 | Neg0],
> Zeros1 = Zeros0
> )
> ;
> Pos1 = Pos0,
> Neg1 = Neg0,
> Zeros1 = [Vec0 | Zeros0]
> ),
> separate_vectors_acc(Vectors, Var, Pos1, Neg1, Zeros1, Pos, Neg, Zeros).
>
>
> % Fails if the variable searched for is not in the vector
> % (i.e. if it has coefficient zero).
> :- pred get_var_coeff(vector, var, rat).
> :- mode get_var_coeff(in, in, out) is semidet.
> get_var_coeff(pair(Varmap, _), Var, Num) :-
> map__search(Varmap, Var, Num).
> :- pred add_vectors(vector, vector, vector).
> :- mode add_vectors(in, in, out) is det.
>
>
> add_vectors(pair(Map0, Float0), pair(Map1, Float1), pair(Map2, Float2)) :-
> Float2 = Float0 + Float1,
> Is_map_key = lambda([Var::out] is nondet, (
> map__member(Map0, Var, _)
> )),
> Add_val = lambda([Var::in, Map_i::in, Map_o::out] is det, (
> map__lookup(Map0, Var, Num0),
> ( map__search(Map_i, Var, Num1) ->
> ( Num0 + Num1 = zero ->
> map__delete(Map_i, Var, Map_o)
> ;
> map__det_update(Map_i, Var, Num0+Num1, Map_o)
> )
> ;
> map__det_insert(Map_i, Var, Num0, Map_o)
> )
> )),
> aggregate(Is_map_key, Add_val, Map1, Map2).
>
> % Eliminates a variable V from three sets of linear inequalities,
> % assuming that the inequations in the first list have a positive
> % coefficient for V, and that the second and third lists contain negative
> % and zero coefficients respectively.
> :- pred eliminate_var(var, matrix, matrix, matrix, matrix).
> :- mode eliminate_var(in, in, in, in, out) is det.
> eliminate_var(_, [], _, Matrix, Matrix).
> eliminate_var(Var, [Vec_P|Pos], Neg, Zeros, Result) :-
> Add_vec = lambda([Vec0::in, Vec1::out] is semidet, (
> add_vectors(Vec_P, Vec0, Vec1),
> Vec1 = pair(Map1, Const),
> not ( map__is_empty(Map1), Const >= zero )
> )),
> list__filter_map(Add_vec, Neg, Zeros1),
> list__append(Zeros1, Zeros, New_zeros),
> eliminate_var(Var, Pos, Neg, New_zeros, Result).
>
>
> :- pred eliminate_vars(list(var), matrix, matrix).
> :- mode eliminate_vars(in, in, out) is det.
> eliminate_vars([], Matrix, Matrix).
> eliminate_vars([Var|Vars], Matrix, Result) :-
> separate_vectors(Matrix, Var, Pos, Neg, Zeros),
> eliminate_var(Var, Pos, Neg, Zeros, New_zeros1),
> list__sort_and_remove_dups(New_zeros1, New_zeros2),
> remove_some_redundancies(New_zeros2, [], New_zeros),
> eliminate_vars(Vars, New_zeros, Result).
>
>
> % Removes all constraints of the form a1x1 + ... + anxn <= A
> % where there exists another constraint a1x1 + ... + anxn <= B
> % and B <= A.
> :- pred remove_some_redundancies(matrix, matrix, matrix).
> :- mode remove_some_redundancies(in, in, out) is det.
> remove_some_redundancies([], Constraints, Constraints).
> remove_some_redundancies([Constr0|Equations0], Constrs, Equations) :-
> find_strongest_remove_weaker(Constr0, Constr, Equations0, Equations1),
> remove_some_redundancies(Equations1, [Constr|Constrs], Equations).
>
>
> % Removes all equations in the given list which are
> % weaker than the given constraint Vec0. Returns the strongest vector
> % which is comparable to Vec0, and the list with the redundant
> % constraints removed.
> :- pred find_strongest_remove_weaker(vector, vector, matrix, matrix).
> :- mode find_strongest_remove_weaker(in, out, in, out) is det.
> find_strongest_remove_weaker(Vec0, Vec, Matrix0, Matrix) :-
> Comp_constraints = lambda([Vec1::in, Str_vec0::in, Str_vec1::out,
> List1::in, List2::out] is det, (
> ( compare_constraints(Str_vec0, Vec1, Str_vec) ->
> ( Str_vec = Str_vec0 ->
> List2 = List1,
> Str_vec1 = Str_vec0
> ;
> List2 = List1,
> Str_vec1 = Vec1
> )
> ;
> List2 = [Vec1|List1],
> Str_vec1 = Str_vec0
> )
> )),
> list__foldl2(Comp_constraints, Matrix0, Vec0, Vec, [], Matrix).
>
>
> % Fails if neither constraint implies the other.
> % Otherwise, returns the stronger constraint.
> :- pred compare_constraints(vector, vector, vector).
> :- mode compare_constraints(in, in, out) is semidet.
> compare_constraints(pair(Map0,Float0), pair(Map1,Float1), pair(Map3,Float)) :-
> map__keys(Map0, Keylist),
> Del_same_keys = lambda([Var0::in, Map4::in, Map5::out] is semidet, (
> map__lookup(Map0, Var0, Num0),
> map__remove(Map4, Var0, Num4, Map5),
> Num4 = Num0
> )),
> list__foldl(Del_same_keys, Keylist, Map1, Map2),
> map__is_empty(Map2),
> Map3 = Map0,
> ( Float0 < Float1 ->
> Float = Float0
> ;
> Float = Float1
> ).
>
>
>
>
>
>
>
> %------------------------------------------------------------------------------%
> %------------------------------------------------------------------------------%
>
> % Convex Hull
>
> :- type var_info
> ---> var_info(
> list(map(var,var)), % Maps from original variables
> % to new (temporary) ones.
> % A variables which occurs in
> % more than one polygon is
> % mapped to a separate variable
> % for each one.
> % This list contains one map
> % for each polygon.
>
> list(var), % List of sigma variables.
>
> varset
> ).
>
> :- type eqn_info
> ---> eqn_info(
> map(var, var),
> varset
> ).
>
> % convex_hull takes a list of polyhedra (represented as lists of constraints)
> % and computes the smallest convex set which is a superset of all of them.
> % For details of the algorithm, see Benoy and King.
> convex_hull(Polys0, Poly, Vars1, Vars) :-
> list__map(fm_standardize_equations, Polys0, Polys),
> transform_polys(Polys, Eqns1, var_info([], [], Vars1),
> var_info(Map_list, Sigmas, Vars)),
> add_sigma_eqns(Eqns1, Sigmas, Eqns2),
> add_last_eqns(Eqns2, Map_list, Eqns3),
> Append_values = lambda([Map::in, Varlist0::in, Varlist::out] is det, (
> map__values(Map, List),
> append(List, Varlist0, Varlist)
> )),
> list__foldl(Append_values, Map_list, [], Values),
> append(Sigmas, Values, Variables),
> project(Eqns3, Variables, Poly).
>
>
> :- pred transform_polys(list(equations), equations, var_info, var_info).
> :- mode transform_polys(in, out, in, out) is det.
> transform_polys(Polys, Eqns, Var_info0, Var_info) :-
> list__foldl2(transform_poly, Polys, [], Eqns, Var_info0, Var_info).
>
>
> % transform_poly: takes a polygon (with original variables), a list
> % of the polygons transformed so far and variable information.
> % It substitutes new (temporary) variables for all of the variables
> % in the polygon while constructing a map from old variables to new
> % ones. The new polygon is then added to the list of transformed ones.
> % It then updates the variable info by adding the new map and a new
> % sigma variable.
> :- pred transform_poly(equations, equations, equations, var_info, var_info).
> :- mode transform_poly(in, in, out, in, out) is det.
> transform_poly(Poly, Polys0, Polys, var_info(Maps0, Sigmas0, Vars0), Var_info):-
> map__init(Newmap),
> varset__new_var(Vars0, Sigma, Vars1),
> Trans_eqn = lambda([Eqn0::in, Eqn::out, Eq_inf0::in,
> Eq_inf1::out] is det, (
> transform_eqn(Eqn0, Eqn, Sigma, Eq_inf0, Eq_inf1)
> )),
> list__map_foldl(Trans_eqn, Poly, New_eqns, eqn_info(Newmap, Vars1),
> eqn_info(Map, Vars)),
> append(New_eqns, Polys0, Polys),
> Var_info = var_info([Map|Maps0], [Sigma|Sigmas0], Vars).
>
>
> :- pred transform_eqn(equation, equation, var, eqn_info, eqn_info).
> :- mode transform_eqn(in, out, in, in, out) is det.
> transform_eqn(eqn(Coeffs0, Op, Num0), eqn(Coeffs, Op, Num), Sigma,
> eqn_info(Map0, Vars0), eqn_info(Map, Vars)) :-
> list__map_foldl(change_var, Coeffs0, Coeffs1,
> eqn_info(Map0, Vars0), eqn_info(Map, Vars)),
> Coeffs = [Sigma - -Num0|Coeffs1],
> Num = zero.
>
> % change_var: takes a Var-Num pair with an old variable and returns one with
> % a new variable which the old variable maps to. Updates the map of old to
> % new variables if necessary.
> :- pred change_var(coeff, coeff, eqn_info, eqn_info).
> :- mode change_var(in, out, in, out) is det.
> change_var(Var0-Num, Var-Num, eqn_info(Map0, Varset0), eqn_info(Map, Varset)) :-
> % Have we already mapped this original variable to a new one?
> ( map__search(Map0, Var0, Var2) ->
> Var = Var2,
> Map = Map0,
> Varset = Varset0
> ;
> varset__new_var(Varset0, Nvar, Varset),
> map__det_insert(Map0, Var0, Nvar, Map),
> Var = Nvar
> ).
>
>
>
> % add_sigma_eqns: Takes the list of equations so far and appends
> % the non-negativity equation for each sigma variable,
> % and an equation stating that the sum of the sigmas is one.
> :- pred add_sigma_eqns(equations, list(var), equations).
> :- mode add_sigma_eqns(in, in, out) is det.
> add_sigma_eqns(Eqns0, Sigmas, Eqns) :-
> Non_neg = lambda([Var::in, Eqn::out] is det, (
> Eqn = eqn([Var - -one], (=<), zero)
> )),
> list__map(Non_neg, Sigmas, Eqns1),
> append(Eqns1, Eqns0, Eqns2),
> Var_to_coeff = lambda([Var::in, Coeff::out] is det, (
> Coeff = Var-one
> )),
> list__map(Var_to_coeff, Sigmas, Coefflist),
> fm_standardize_equations([eqn(Coefflist, (=), one)], Eqns3),
> append(Eqns3, Eqns2, Eqns).
>
> % add_last_eqns: Adds the equations stating that each original
> % variable is the sum of the temporary variables to which it has
> % been mapped.
> :- pred add_last_eqns(equations, list(map(var,var)), equations).
> :- mode add_last_eqns(in, in, out) is det.
> add_last_eqns(Eqns0, Old_to_new_var_maps, Eqns) :-
> get_map_keys(Old_to_new_var_maps, [], Keys1),
> VarComp = lambda([Var1::in, Var2::in, Res::out] is det, (
> compare(Res, Var1, Var2)
> )),
> list__sort_and_remove_dups(VarComp, Keys1, Keys),
> Original_var_to_eqn = lambda([Orig_var::in, Eqn::out] is semidet, (
> original_var_to_eqn(Orig_var, Old_to_new_var_maps, Eqn)
> )),
> list__filter_map(Original_var_to_eqn, Keys, Eqns1),
> fm_standardize_equations(Eqns1, Eqns2),
> append(Eqns2, Eqns0, Eqns).
>
>
> :- pred original_var_to_eqn(var, list(map(var,var)), equation).
> :- mode original_var_to_eqn(in, in, out) is semidet.
> original_var_to_eqn(Original_var, Maplist, Eqn) :-
> Corresp_vars = lambda([Map::in, Coeffs0::in, Coeffs::out] is semidet, (
> map__search(Map, Original_var, Newvar),
> Coeffs = [Newvar - -one|Coeffs0]
> )),
> list__foldl(Corresp_vars, Maplist, [], Coefflist),
> Eqn = eqn([Original_var-one|Coefflist], (=), zero).
>
>
> :- pred get_map_keys(list(map(var,var)), list(var), list(var)).
> :- mode get_map_keys(in, in, out) is det.
> get_map_keys([], Keys, Keys).
> get_map_keys([Map|Maps], Keys0, Keys) :-
> map__keys(Map, Keys1),
> append(Keys1, Keys0, NewKeys),
> get_map_keys(Maps, NewKeys, Keys).
>
> %-----------------------------------------------------------------------------%
> % Widening
>
> % widen(Cs1, Cs2, Cs, Vars): Relaxes the constraints Cs1 by selecting
> % those constraints from Cs1 which are implied by Cs2.
> widen(Poly1, Poly2, Vars, Wide_poly) :-
> Entailed_by_Poly2 = lambda([Constraint::in] is semidet, (
> entailed(Constraint, Poly2, Vars)
> )),
> list__filter(Entailed_by_Poly2, Poly1, Wide_poly).
>
>
> % entailed(C, Cs, Vars): Determines whether the constraint C is
> % implied by the set of constraints Cs.
> % Uses the simplex method to find the point P satisfying Cs
> % which maximises (or minimises if C contains >= ) a function parallel
> % to C. Then tests whether P satisfies C.
> % This assumes that all the variables are non-negative.
> entailed(eqn(Coeff_list, (=<), Const), Poly, Vars) :-
> Objective = Coeff_list,
> lp_solve(Poly, max, Objective, Vars, [], Result),
> (
> Result = satisfiable(Max_val, _),
> Max_val =< Const
> ;
> Result = inconsistent,
> error("inconsistent polygon passed to entailed")
> ).
>
>
> entailed(eqn(Coeff_list, (>=), Const), Poly, Vars) :-
> Objective = Coeff_list,
> lp_solve(Poly, min, Objective, Vars, [], Result),
> (
> Result = satisfiable(Min_val, _),
> Min_val >= Const
> ;
> Result = inconsistent,
> error("inconsistent polygon passed to entails")
> ).
>
>
> entailed(eqn(Coeff_list, (=), Const), Poly, Vars) :-
> entailed(eqn(Coeff_list, (=<), Const), Poly, Vars),
> entailed(eqn(Coeff_list, (>=), Const), Poly, Vars).
>
> :- pred nz(rat::in) is semidet.
>
> nz(X) :- X \= zero.
>
More information about the developers
mailing list