mercurylang/mercury
1
0
Fork 0
mirror of https://github.com/Mercury-Language/mercury.git synced 2026-04-15 01:13:30 +00:00
Code Issues Projects Releases Wiki Activity
Files
cdf10e673cf66188e3facae0f458aea71f9ee438
mercury/compiler/lp_rational.m
Zoltan Somogyi 97b9084d4e Rename some functions to avoid ambiguity. 
compiler/lp_rational.m:
compiler/polyhedron.m:
    Rename a function in each module.

compiler/term_constr_fixpoint.m:
compiler/term_constr_pass2.m:
    Conform to the renames.
2025-08-10 20:31:46 +02:00

2486 lines
84 KiB
Mathematica
Raw Blame History

%---------------------------------------------------------------------------%
% vim: ft=mercury ts=4 sw=4 et
%---------------------------------------------------------------------------%
% Copyright (C) 1997-2002, 2005-2007, 2009-2012 The University of Melbourne.
% Copyright (C) 2015-2021, 2023-2025 The Mercury team.
% This file may only be copied under the terms of the GNU General
% Public License - see the file COPYING in the Mercury distribution.
%---------------------------------------------------------------------------%
%
% File: lp_rational.m.
% Main authors: conway, juliensf, vjteag.
%
% This module contains code for creating and manipulating systems of rational
% linear arithmetic constraints. It provides the following operations:
%
% * optimization (using the simplex method)
%
% * projection (using Fourier elimination).
%
% * an entailment test (using the above linear optimizer).
%
%---------------------------------------------------------------------------%
:- module libs.lp_rational.
:- interface.
:- import_module libs.rat.
:- import_module io.
:- import_module list.
:- import_module map.
:- import_module maybe.
:- import_module pair.
:- import_module set.
:- import_module term.
:- import_module varset.
%---------------------------------------------------------------------------%
%
% Linear constraints over Q^n.
%
:- type lp_constant == rat.
:- type lp_coefficient == rat.
:- type lp_var == var.
:- type lp_varset == varset.
:- type lp_term == pair(lp_var, lp_coefficient).
% Create a term with a coefficient of 1.
% For use with ho functions.
%
:- func lp_term(lp_var) = lp_term.
:- type lp_operator
---> lp_lt_eq
; lp_eq
; lp_gt_eq.
:- inst lp_op_lt_eq_or_eq for lp_operator/0
---> lp_lt_eq
; lp_eq.
% A primitive linear arithmetic constraint.
%
:- type lp_constraint.
% A conjunction of primitive constraints.
%
:- type lp_constraint_conj == list(lp_constraint).
% Create a constraint from the given components.
%
:- func construct_constraint(list(lp_term), lp_operator, lp_constant) =
lp_constraint.
% Create a constraint from the given components.
% Throws an exception if the resulting constraint is trivially false.
%
:- func construct_non_false_constraint(list(lp_term), lp_operator, lp_constant)
= lp_constraint.
% Deconstruct the given constraint.
%
:- pred deconstruct_constraint(lp_constraint::in,
list(lp_term)::out, lp_operator::out, lp_constant::out) is det.
% As above but throws an exception if the constraint is false.
%
:- pred deconstruct_non_false_constraint(lp_constraint::in,
list(lp_term)::out, lp_operator::out(lp_op_lt_eq_or_eq), lp_constant::out)
is det.
% Succeeds iff the given constraint contains a single variable and
% that variable is constrained to be a nonnegative value.
%
:- pred nonneg_constr(lp_constraint::in) is semidet.
% Create a constraint that constrains the argument
% have a non-negative value.
%
:- func make_nonneg_constr(lp_var) = lp_constraint.
% Create a constraint that equates two variables.
%
:- func make_vars_eq_constraint(lp_var, lp_var) = lp_constraint.
% Create constraints of the form:
%
% Var = Constant or Var >= Constant
%
% These functions are useful with higher-order code.
%
:- func make_var_const_eq_constraint(lp_var, rat) = lp_constraint.
:- func make_var_const_gte_constraint(lp_var, rat) = lp_constraint.
% Create a constraint that is trivially true.
%
:- func true_constraint = lp_constraint.
% Create a constraint that is trivially false.
%
:- func false_constraint = lp_constraint.
% Succeeds if the constraint is trivially true.
%
:- pred is_true(lp_constraint::in) is semidet.
% Succeeds if the constraint is trivially false.
%
:- pred is_false(lp_constraint::in) is semidet.
% Takes a list of constraints and looks for equality constraints
% that may be implicit in any inequalities.
%
% NOTE: this is only a syntactic check so it may miss
% some equalities that are implicit in the system.
%
:- pred restore_equalities(lp_constraint_conj::in, lp_constraint_conj::out)
is det.
% Succeed iff the given system of constraints is inconsistent.
%
:- pred inconsistent(lp_varset::in, lp_constraint_conj::in) is semidet.
% Remove those constraints from the system whose redundancy can be
% trivially detected.
%
% NOTE: the resulting system of constraints may not be minimal.
%
:- func simplify_constraints(lp_constraint_conj) = lp_constraint_conj.
% substitute_vars(RenameMap, Constraints0) = Constraints:
% substitute_corresponding_vars(VarsA, VarsB, Constraints0) = Constraints:
%
% Perform variable substitution on the given system of constraints
% based upon either the explicitly-given RenameMap, or the mapping
% that is implicit between the corresponding elements of VarsA and VarsB.
% If length(VarsA) \= length(VarsB), then throw an exception.
%
% XXX As of 2025 aug 10, substitute_corresponding_vars is unused,
% except by polyhedron.substitute_corresponding_vars, which is
% itself unused.
%
:- func substitute_vars(map(lp_var, lp_var), lp_constraint_conj) =
lp_constraint_conj.
:- func substitute_corresponding_vars(list(lp_var), list(lp_var),
lp_constraint_conj) = lp_constraint_conj.
% Make the values of all the variables in the set zero.
%
:- func set_vars_to_zero(set(lp_var), lp_constraint_conj) = lp_constraint_conj.
%---------------------------------------------------------------------------%
%
% Bounding boxes and other approximations.
%
% Approximate the solution space of a set of constraints using
% a bounding box. If the system is inconsistent then the resulting
% system will also be inconsistent.
%
:- func bounding_box(lp_varset, lp_constraint_conj) = lp_constraint_conj.
% Create non-negativity constraints for all of the variables in the
% given list of constraints, except for the variables specified
% in the first argument.
%
:- func nonneg_box(list(lp_var), lp_constraint_conj) = lp_constraint_conj.
%---------------------------------------------------------------------------%
%
% Linear solver.
%
:- type objective == list(lp_term).
:- type direction
---> max
; min.
:- type lp_result
---> lp_res_unbounded
; lp_res_inconsistent
; lp_res_satisfiable(rat, map(lp_var, rat)).
% lp_res_satisfiable(ObjVal, MapFromObjVarsToVals)
% Maximize (or minimize - depending on `direction') `objective'
% subject to the given constraints. The variables in the objective
% and the constraints *must* be from the given `lp_varset'.
% This is passed to the solver so that it can allocate fresh variables
% as required.
%
% The result is `unbounded' if the objective is not bounded by
% the constraints, `inconsistent' if the given constraints are
% inconsistent, or `satisfiable/2' otherwise.
%
:- func solve(lp_constraint_conj, direction, objective, lp_varset) = lp_result.
%---------------------------------------------------------------------------%
%
% Projection.
%
:- type projection_result
---> pr_res_ok(lp_constraint_conj) % projection succeeded.
; pr_res_inconsistent % matrix was inconsistent.
; pr_res_aborted. % ran out of time/space and backed out.
% project(Constraints0, Vars, VarSet, Result):
%
% Takes a list of constraints, `Constraints0', and eliminates the
% variables in the list `Vars' using Fourier elimination.
%
% Returns `ok(Constraints)' if `Constraints' is the projection
% of `Constraints0' over `Vars'. Returns `inconsistent' if `Constraints0'
% is inconsistent. Returns `aborted' if the intermediate matrices grow
% too large while performing Fourier elimination.
%
% NOTE: this does not always detect that a constraint set is inconsistent,
% so if callers of this procedure require the resulting system of
% constraints to be consistent, they will need to do a consistency check
% on the result themselves.
%
:- pred project_constraints(lp_varset::in, list(lp_var)::in,
lp_constraint_conj::in, projection_result::out) is det.
% project_constraints_maybe_size_limit(VarSet, MaybeMaxMatrixSize, Vars,
% Matrix, Result):
%
% Same as above, but if MaybeMaxMatrixSize = yes(MaxMatrixSize), then
% if the size of the matrix ever exceeds `MaxMatrixSize', we back out
% of the computation.
%
:- pred project_constraints_maybe_size_limit(lp_varset::in, maybe(int)::in,
list(lp_var)::in, lp_constraint_conj::in, projection_result::out) is det.
%---------------------------------------------------------------------------%
%
% Entailment.
%
:- type entailment_result
---> entailed
; not_entailed
; inconsistent.
% entailed(VarSet, Cs, C):
%
% Determines if the constraint `C' is implied by the set of
% constraints `Cs'. Uses the simplex method to find the point `P'
% satisfying `Cs' which maximizes (if `C' contains '=<') or
% minimizes (if `C' contains '>=') a function parallel to `C'.
% Returns `entailed' if `P' satisfies `C', `not_entailed' if it does not
% and `inconsistent' if `Cs' is not a consistent system of constraints.
%
% This assumes that all variables are non-negative.
%
:- func entailed(lp_varset, lp_constraint_conj, lp_constraint)
= entailment_result.
% entailed(VarSet, Cs, C):
%
% As above but fails if `C' is not implied by `Cs' and
% throws an exception if `Cs' is inconsistent.
%
:- pred entailed(lp_varset::in, lp_constraint_conj::in, lp_constraint::in)
is semidet.
% Check if a constraint is entailed by all the others in the set.
% If it is, then remove it from the set.
%
% NOTE: this can be very slow - also due to the order in which
% the constraints are processed, it may not produce a minimal set.
%
% Fails if the system of constraints is inconsistent.
%
:- pred remove_some_entailed_constraints(lp_varset::in,
lp_constraint_conj::in, lp_constraint_conj::out) is semidet.
%---------------------------------------------------------------------------%
%
% Stuff for intermodule optimization.
%
% A function that converts an lp_var into a string.
% XXX This is *not* a good name for this type.
%
:- type output_var == (func(lp_var) = string).
% Write out the constraints in a form we can read in using the
% term parser from the standard library.
%
:- pred output_constraints(io.text_output_stream::in, output_var::in,
lp_constraint_conj::in, io::di, io::uo) is det.
%---------------------------------------------------------------------------%
%
% Debugging predicates.
%
% Return the set of variables that are present in a list of constraints.
%
% XXX This shouldn't be exported but it is currently needed by the
% workaround for the problem with head variables in term_constr_fixpoint.m.
%
:- func get_vars_from_constraints(lp_constraint_conj) = set(lp_var).
% Print out the constraints using the names in the varset. If the variable
% has no name it will be given the name Temp<n>, where <n> is the
% variable number.
%
:- pred write_constraints(io.text_output_stream::in, lp_varset::in,
lp_constraint_conj::in, io::di, io::uo) is det.
%---------------------------------------------------------------------------%
%---------------------------------------------------------------------------%
:- implementation.
:- import_module parse_tree.
:- import_module parse_tree.parse_tree_output.
% undesirable dependency, for write_out_list
:- import_module assoc_list.
:- import_module bool.
:- import_module int.
:- import_module require.
:- import_module solutions.
:- import_module string.
%---------------------------------------------------------------------------%
%
% Constraints.
%
% The following properties should hold for each constraint:
% - there is one instance of each variable in the term list.
% - the terms are sorted in increasing order by variable.
% - the terms should be normalized so that the leading term
% has a coefficient of +/-1 (unless all terms have a coefficient
% of zero - in which case the term list is empty).
% - variables with coefficient zero are *not* included in the list
% of terms.
:- type lp_constraint
---> lte(list(lp_term), lp_constant) % sumof(Terms) =< Constant
; eq(list(lp_term), lp_constant) % sumof(Terms) = Constant
; gte(list(lp_term), lp_constant). % sumof(Terms) >= Constant
%---------------------------------------------------------------------------%
%
% Procedures for constructing/deconstructing constraints.
%
lp_term(Var) = Var - one.
construct_constraint(Terms0, Op, Const0) = Constraint :-
(
Terms0 = [],
(
Op = lp_lt_eq,
Constraint = lte([], Const0)
;
Op = lp_eq,
Constraint = eq([], Const0)
;
Op = lp_gt_eq,
Constraint = lte([], -Const0)
)
;
Terms0 = [_ | _],
(
Op = lp_lt_eq,
Terms1 = sum_like_terms(Terms0),
normalize_terms_and_const(yes, Terms1, Const0, Terms, Const),
Constraint = lte(Terms, Const)
;
Op = lp_eq,
Terms1 = sum_like_terms(Terms0),
normalize_terms_and_const(no, Terms1, Const0, Terms, Const),
Constraint = eq(Terms, Const)
;
Op = lp_gt_eq,
Terms1 = sum_like_terms(Terms0),
normalize_terms_and_const(yes, Terms1, Const0, Terms, Const),
Constraint = lte(negate_lp_terms(Terms), -Const)
)
).
% This is for internal use only - it builds a constraint out of the parts
% but does *not* attempt to perform any standardization. It is intended for
% use in operations such as normalization.
%
:- func unchecked_construct_constraint(list(lp_term), lp_operator,
lp_constant) = lp_constraint.
unchecked_construct_constraint(Terms, lp_lt_eq, Constant) =
lte(Terms, Constant).
unchecked_construct_constraint(Terms, lp_eq, Constant) =
eq(Terms, Constant).
unchecked_construct_constraint(Terms, lp_gt_eq, Constant) =
gte(Terms, Constant).
:- func sum_like_terms(list(lp_term)) = list(lp_term).
sum_like_terms(Terms) = map.to_assoc_list(lp_terms_to_map(Terms)).
% Convert an association list of lp_vars and coefficients to a map
% of the same. If there are duplicate keys in the list make sure that
% eventual value in the map is the sum of the two coefficients.
% Also if a coefficient is (or ends up being) zero, make sure that
% the variable doesn't end up in the resulting map.
%
:- func lp_terms_to_map(assoc_list(lp_var, lp_coefficient)) =
map(lp_var, lp_coefficient).
lp_terms_to_map(Terms) = Map :-
list.foldl(lp_terms_to_map_2, Terms, map.init, Map).
:- pred lp_terms_to_map_2(pair(lp_var, lp_coefficient)::in,
map(lp_var, lp_coefficient)::in, map(lp_var, lp_coefficient)::out) is det.
lp_terms_to_map_2(Var - Coeff0, !Map) :-
( if map.search(!.Map, Var, MapCoeff) then
Coeff = MapCoeff + Coeff0,
( if Coeff = zero then
map.delete(Var, !Map)
else
map.set(Var, Coeff, !Map)
)
else
( if Coeff0 = zero then
true
else
map.set(Var, Coeff0, !Map)
)
).
construct_non_false_constraint(Terms, Op, Constant) = Constraint :-
Constraint = construct_constraint(Terms, Op, Constant),
( if is_false(Constraint) then
unexpected($pred, "false constraint")
else
true
).
deconstruct_constraint(lte(Terms, Constant), Terms, lp_lt_eq, Constant).
deconstruct_constraint(eq(Terms, Constant), Terms, lp_eq, Constant).
deconstruct_constraint(gte(Terms, Constant), Terms, lp_gt_eq, Constant).
deconstruct_non_false_constraint(Constraint, Terms, Operator, Constant) :-
( if is_false(Constraint) then
unexpected($pred, "false_constraint")
else
true
),
(
Constraint = lte(Terms, Constant),
Operator = lp_lt_eq
;
Constraint = eq(Terms, Constant),
Operator = lp_eq
;
Constraint = gte(_, _),
unexpected($pred, "gte encountered")
).
:- func lp_terms(lp_constraint) = list(lp_term).
lp_terms(lte(Terms, _)) = Terms.
lp_terms(eq(Terms, _)) = Terms.
lp_terms(gte(Terms, _)) = Terms.
:- func constant(lp_constraint) = lp_constant.
constant(lte(_, Constant)) = Constant.
constant(eq(_, Constant)) = Constant.
constant(gte(_, Constant)) = Constant.
:- func operator(lp_constraint) = lp_operator.
operator(lte(_, _)) = lp_lt_eq.
operator(eq(_, _)) = lp_eq.
operator(gte(_,_)) = unexpected($pred, "gte").
:- func negate_operator(lp_operator) = lp_operator.
negate_operator(lp_lt_eq) = lp_gt_eq.
negate_operator(lp_eq) = lp_eq.
negate_operator(lp_gt_eq) = lp_lt_eq.
nonneg_constr(lte([_ - (-rat.one)], rat.zero)).
nonneg_constr(gte(_, _)) :-
unexpected($pred, "gte").
make_nonneg_constr(Var) =
construct_constraint([Var - (-rat.one)], lp_lt_eq, rat.zero).
make_vars_eq_constraint(Var1, Var2) =
construct_constraint([Var1 - rat.one, Var2 - (-rat.one)], lp_eq, rat.zero).
make_var_const_eq_constraint(Var, Constant) =
construct_constraint([Var - rat.one], lp_eq, Constant).
make_var_const_gte_constraint(Var, Constant) =
construct_constraint([Var - rat.one], lp_gt_eq, Constant).
true_constraint = eq([], rat.zero).
false_constraint = eq([], rat.one).
is_true(gte([], Const)) :- Const =< rat.zero.
is_true(lte([], Const)) :- Const >= rat.zero.
is_true(eq([], Const)) :- Const = rat.zero.
is_false(gte([], Const)) :- Const > rat.zero.
is_false(lte([], Const)) :- Const < rat.zero.
is_false(eq([], Const)) :- Const \= rat.zero.
%---------------------------------------------------------------------------%
restore_equalities([], []).
restore_equalities([E0 | Es0], [E | Es]) :-
( if check_for_equalities(E0, Es0, [], E1, Es1) then
E = E1,
Es2 = Es1
else
Es2 = Es0,
E = E0
),
restore_equalities(Es2, Es).
:- pred check_for_equalities(lp_constraint::in, lp_constraint_conj::in,
lp_constraint_conj::in,
lp_constraint::out, lp_constraint_conj::out) is semidet.
check_for_equalities(Eqn0, [Eqn | Eqns], SoFar, NewEqn, NewEqnSet) :-
( if opposing_inequalities(Eqn0 @ lte(Coeffs, Constant), Eqn) then
NewEqn = standardize_constraint(eq(Coeffs, Constant)),
NewEqnSet = SoFar ++ Eqns
else
check_for_equalities(Eqn0, Eqns, [Eqn | SoFar], NewEqn, NewEqnSet)
).
% Checks if a pair of constraints are inequalities of the form:
%
% -ax1 - ax2 - ... - axN =< -C
% ax1 + ax2 + ... + axN =< C
%
% These can be converted into the equality:
%
% ax1 + ... + axN = C
%
% NOTE: we don't check for gte constraints because these should
% have been transformed away when we converted to standard form.
%
:- pred opposing_inequalities(lp_constraint::in, lp_constraint::in) is semidet.
opposing_inequalities(lte(TermsA, Const), lte(TermsB, -Const)) :-
TermsB = list.map((func(V - X) = V - (-X)), TermsA).
%---------------------------------------------------------------------------%
% Put a constraint into standard form. Every constraint has its terms list
% in increasing order of variable name and then multiplied so that
% the absolute value of the leading coefficient is one.
% op_ge is converted to op_le by multiplying through by negative one.
% op_eq constraints should have an initial coefficient of (positive) 1.
%
:- func standardize_constraint(lp_constraint) = lp_constraint.
standardize_constraint(gte(Terms0, Const0)) = Constraint :-
normalize_terms_and_const(yes, Terms0, Const0, Terms, Const),
Constraint = lte(negate_lp_terms(Terms), -Const).
standardize_constraint(eq(Terms0, Const0)) = eq(Terms, Const) :-
normalize_terms_and_const(no, Terms0, Const0, Terms, Const).
standardize_constraint(lte(Terms0, Const0)) = lte(Terms, Const) :-
normalize_terms_and_const(yes, Terms0, Const0, Terms, Const).
% Sort the list of terms in ascending order by variable and then
% multiply through so that the first term has a coefficient of
% one or negative one. If the first argument is `yes', then we multiply
% through by the reciprocal of the absolute value of the coefficient,
% otherwise we multiply through by the reciprocal of the value.
%
:- pred normalize_terms_and_const(bool::in, list(lp_term)::in, lp_constant::in,
list(lp_term)::out, lp_constant::out) is det.
normalize_terms_and_const(AbsVal, !.Terms, !.Const, !:Terms, !:Const) :-
CompareTerms =
( func(VarA - _, VarB - _) = Result :-
compare(Result, VarA, VarB)
),
!:Terms = list.sort(CompareTerms, !.Terms),
(
!.Terms = [_ - Coefficient0 | _],
(
AbsVal = yes,
Coefficient = rat.abs(Coefficient0)
;
AbsVal = no,
Coefficient = Coefficient0
),
( if Coefficient = rat.zero then
unexpected($pred, "zero coefficient")
else
true
),
DivideBy = (func(Var - Coeff) = Var - (Coeff / Coefficient)),
!:Terms = list.map(DivideBy, !.Terms),
!:Const = !.Const / Coefficient
;
!.Terms = []
).
% Succeeds iff the constraint is implied by the assumption that
% all variables are non-negative *and* the constraint is not one
% used to force non-negativity of the variables.
%
:- pred obvious_constraint(lp_constraint::in) is semidet.
obvious_constraint(lte(Terms, Constant)) :-
Constant >= rat.zero,
list.length(Terms) >= 2,
all [Term] list.member(Term, Terms) => snd(Term) < zero.
obvious_constraint(gte(Terms, Constant)) :-
Constant =< rat.zero,
list.length(Terms) >= 2,
all [Term] (
list.member(Term, Terms)
=>
snd(Term) > zero
).
inconsistent(Vars, Constraints @ [Constraint | _]) :-
(
is_false(Constraint)
;
( Constraint = lte([Term | _], _)
; Constraint = eq([Term | _], _)
; Constraint = gte([Term | _], _)
),
DummyObjective = [Term],
lp_rational.solve(Constraints, max, DummyObjective, Vars) =
lp_res_inconsistent
).
simplify_constraints(Constraints) = remove_weaker(remove_trivial(Constraints)).
:- func remove_trivial(lp_constraint_conj) = lp_constraint_conj.
remove_trivial([]) = [].
remove_trivial([Constraint | Constraints]) = Result :-
( if is_false(Constraint) then
Result = [false_constraint]
else
Result0 = remove_trivial(Constraints),
( if
Result0 = [C],
is_false(C)
then
Result = Result0
else
% Remove the constraint if it is trivially true or the result
% of all the variables being non-negative.
( if
( is_true(Constraint)
; obvious_constraint(Constraint)
)
then
Result = Result0
else
Result = [Constraint | Result0]
)
)
).
:- func remove_weaker(lp_constraint_conj) = lp_constraint_conj.
remove_weaker([]) = [].
remove_weaker([C | Cs0]) = Result :-
list.foldl2(remove_weaker_2(C), Cs0, [], Cs, yes, Keep),
Result0 = remove_weaker(Cs),
(
Keep = yes,
Result = [C | Result0]
;
Keep = no,
Result = Result0
).
:- pred remove_weaker_2(lp_constraint::in, lp_constraint::in,
lp_constraint_conj::in, lp_constraint_conj::out,
bool::in, bool::out) is det.
remove_weaker_2(A, B, !Acc, !Keep) :-
( if is_stronger(A, B) then
true
else if is_stronger(B, A) then
list.cons(B, !Acc),
!:Keep = no
else
list.cons(B, !Acc)
).
:- pred is_stronger(lp_constraint::in, lp_constraint::in) is semidet.
is_stronger(eq(Terms, Const), gte(Terms, Const)).
is_stronger(eq(Terms, Const), lte(Terms, Const)).
is_stronger(eq(Terms, Const), gte(negate_lp_terms(Terms), -Const)).
is_stronger(eq(Terms, Const), lte(negate_lp_terms(Terms), -Const)).
is_stronger(lte([Var - (-one)], ConstA), lte([Var - (-one)], ConstB)) :-
ConstA =< zero, ConstA =< ConstB.
is_stronger(eq(Terms, ConstA), lte(negate_lp_terms(Terms), ConstB)) :-
ConstA >= (-one) * ConstB.
is_stronger(lte(Terms, ConstA), lte(Terms, ConstB)) :-
ConstB =< zero, ConstA =< ConstB.
substitute_vars(RenameMap, Constraints0) = Constraints :-
Constraints = list.map(substitute_vars_in_constraint(RenameMap),
Constraints0).
substitute_corresponding_vars(Old, New, Constraints0) = Constraints :-
map.from_corresponding_lists(Old, New, RenameMap),
Constraints = list.map(substitute_vars_in_constraint(RenameMap),
Constraints0).
:- func substitute_vars_in_constraint(map(lp_var, lp_var), lp_constraint)
= lp_constraint.
substitute_vars_in_constraint(RenameMap, lte(Terms0, Const)) = Result :-
Terms = list.map(substitute_term(RenameMap), Terms0),
Result = lte(sum_like_terms(Terms), Const).
substitute_vars_in_constraint(RenameMap, eq(Terms0, Const)) = Result :-
Terms = list.map(substitute_term(RenameMap), Terms0),
Result = eq(sum_like_terms(Terms), Const).
substitute_vars_in_constraint(_, gte(_, _)) =
unexpected($pred, "gte").
:- func substitute_term(map(lp_var, lp_var), lp_term) = lp_term.
substitute_term(RenameMap, Term0) = Term :-
Term0 = Var0 - Coeff,
map.lookup(RenameMap, Var0, Var),
Term = Var - Coeff.
set_vars_to_zero(Vars, Constraints) =
list.map(set_vars_to_zero_2(Vars), Constraints).
:- func set_vars_to_zero_2(set(lp_var), lp_constraint) = lp_constraint.
set_vars_to_zero_2(Vars, lte(Terms0, Const)) = lte(Terms, Const) :-
Terms = set_terms_to_zero(Vars, Terms0).
set_vars_to_zero_2(Vars, eq(Terms0, Const)) = eq(Terms, Const) :-
Terms = set_terms_to_zero(Vars, Terms0).
set_vars_to_zero_2(Vars, gte(Terms0, Const)) = gte(Terms, Const) :-
Terms = set_terms_to_zero(Vars, Terms0).
:- func set_terms_to_zero(set(lp_var), list(lp_term)) = list(lp_term).
set_terms_to_zero(Vars, Terms0) = Terms :-
IsNonZero =
( pred(Term::in) is semidet :-
Term = Var - _Coeff,
not set.member(Var, Vars)
),
Terms = list.filter(IsNonZero, Terms0).
%---------------------------------------------------------------------------%
%
% Bounding boxes and other weaker approximations of the convex union.
%
bounding_box(VarSet, Constraints) = BoundingBox :-
Vars = set.to_sorted_list(get_vars_from_constraints(Constraints)),
CallProject =
( func(Var, Constrs0) = Constrs :-
project_constraints(VarSet, [Var], Constrs0, Result),
(
Result = pr_res_inconsistent,
Constrs = [false_constraint]
;
% If we needed to abort this computation we will just
% approximate the whole lot by `true'.
Result = pr_res_aborted,
Constrs = []
;
Result = pr_res_ok(Constrs)
)
),
BoundingBox = list.foldl(CallProject, Vars, Constraints).
nonneg_box(VarsToIgnore, Constraints) = NonNegConstraints :-
Vars0 = get_vars_from_constraints(Constraints),
MakeConstr =
( pred(Var::in, !.C::in, !:C::out) is det :-
( if list.member(Var, VarsToIgnore) then
true
else
list.cons(make_nonneg_constr(Var), !C)
)
),
set.fold(MakeConstr, Vars0, [], NonNegConstraints).
%---------------------------------------------------------------------------%
%---------------------------------------------------------------------------%
%
% Linear solver.
%
% XXX Most of this came from lp.m. We should try to remove a lot of
% nondeterminism here.
:- type lpr_info
---> lpr_info(
lpr_varset :: lp_varset,
lpr_slack_vars :: list(lp_var),
lpr_art_vars :: list(lp_var)
).
solve(Constraints, Direction, Objective, VarSet) = Result :-
Info0 = lpr_info_init(VarSet),
solve_2(Constraints, Direction, Objective, Result, Info0, _).
% solve_2(Eqns, Dir, Obj, Res, LPRInfo0, LPRInfo) takes
% a list of inequalities `Eqns', a direction for optimization `Dir',
% an objective function `Obj' and an lpr_info structure `LPRInfo0'.
% See inline comments for details on the algorithm.
%
:- pred solve_2(lp_constraint_conj::in, direction::in, objective::in,
lp_result::out, lpr_info::in, lpr_info::out) is det.
solve_2(!.Constraints, Direction, !.Objective, Result, !LPRInfo) :-
% Simplify the inequalities and convert them to standard form by
% introducing slack/artificial variables.
Obj = !.Objective,
lp_standardize_constraints(!Constraints, !LPRInfo),
% If we are maximizing the objective function then we need to negate
% all the coefficients in the objective.
(
Direction = max,
ObjTerms = negate_constraint(eq(!.Objective, zero)),
!:Objective = lp_terms(ObjTerms)
;
Direction = min
),
Rows = list.length(!.Constraints),
Vars = collect_vars(!.Constraints, Obj),
VarList = set.to_sorted_list(Vars),
Columns = list.length(VarList),
VarNums = number_vars(VarList, 0),
ArtVars = !.LPRInfo ^ lpr_art_vars,
Tableau0 = init_tableau(Rows, Columns, VarNums),
insert_constraints(!.Constraints, 1, Columns, VarNums, Tableau0, Tableau),
(
ArtVars = [_ | _],
% There are one or more artificial variables, so we use
% the two-phase method for solving the system.
Result0 = two_phase(Obj, !.Objective, ArtVars, VarNums, Tableau)
;
ArtVars = [],
Result0 = one_phase(Obj, !.Objective, VarNums, Tableau)
),
(
Direction = max,
Result = Result0
;
Direction = min,
(
( Result0 = lp_res_unbounded
; Result0 = lp_res_inconsistent
),
Result = Result0
;
Result0 = lp_res_satisfiable(OptVal, OptCoffs),
Result = lp_res_satisfiable(-OptVal, OptCoffs)
)
).
%---------------------------------------------------------------------------%
:- func one_phase(list(lp_term), list(lp_term), map(lp_var, int), tableau)
= lp_result.
one_phase(Obj0, Obj, VarNums, !.Tableau) = Result :-
insert_terms(Obj, 0, VarNums, !Tableau),
get_vars_from_terms(Obj0, set.init, ObjVars0),
ObjVars = set.to_sorted_list(ObjVars0),
optimize(ObjVars, Result, !.Tableau, _).
%---------------------------------------------------------------------------%
:- func two_phase(list(lp_term), list(lp_term), list(lp_var),
map(lp_var, int), tableau) = lp_result.
two_phase(Obj0, Obj, ArtVars, VarNums, !.Tableau) = Result :-
% Phase 1: minimize the sum of the artificial variables.
ArtObj = list.map(lp_term, ArtVars),
insert_terms(ArtObj, 0, VarNums, !Tableau),
ensure_zero_obj_coeffs(ArtVars, !Tableau),
optimize(ArtVars, Result0, !Tableau),
(
Result0 = lp_res_unbounded,
Result = lp_res_unbounded
;
Result0 = lp_res_inconsistent,
Result = lp_res_inconsistent
;
Result0 = lp_res_satisfiable(Val, _ArtRes),
( if Val = zero then
fix_basis_and_rem_cols(ArtVars, !.Tableau, Tableau1),
% Phase 2:
% Insert the real objective, zero the objective coefficients
% of the basis variables and optimize the objective.
insert_terms(Obj, 0, VarNums, Tableau1, Tableau2),
BasisVars = get_basis_vars(Tableau2),
ensure_zero_obj_coeffs(BasisVars, Tableau2, Tableau3),
get_vars_from_terms(Obj0, set.init, ObjVars0),
ObjVars = set.to_sorted_list(ObjVars0),
optimize(ObjVars, Result, Tableau3, _)
else
Result = lp_res_inconsistent
)
).
%---------------------------------------------------------------------------%
:- pred lp_standardize_constraints(
lp_constraint_conj::in, lp_constraint_conj::out,
lpr_info::in, lpr_info::out) is det.
lp_standardize_constraints(!Constraints, !LPRInfo) :-
list.map_foldl(lp_standardize_constraint, !Constraints, !LPRInfo).
% standardize_constraint performs the following operations on a
% constraint:
%
% - ensures the constant is >= 0 (multiplying by -1 if necessary)
% - introduces slack and artificial variables
%
:- pred lp_standardize_constraint(lp_constraint::in, lp_constraint::out,
lpr_info::in, lpr_info::out) is det.
lp_standardize_constraint(Constr0 @ lte(Coeffs, Const), Constr, !LPRInfo) :-
( if Const < zero then
Constr1 = negate_constraint(Constr0),
lp_standardize_constraint(Constr1, Constr, !LPRInfo)
else
new_slack_var(Var, !LPRInfo),
Constr = lte([Var - one | Coeffs], Const)
).
lp_standardize_constraint(Eqn0 @ eq(Coeffs, Const), Eqn, !LPRInfo) :-
( if Const < zero then
Eqn1 = negate_constraint(Eqn0),
lp_standardize_constraint(Eqn1, Eqn, !LPRInfo)
else
new_art_var(Var, !LPRInfo),
Eqn = lte([Var - one | Coeffs], Const)
).
lp_standardize_constraint(Eqn0 @ gte(Coeffs, Const), Eqn, !LPRInfo) :-
( if Const < zero then
Eqn1 = negate_constraint(Eqn0),
lp_standardize_constraint(Eqn1, Eqn, !LPRInfo)
else
new_slack_var(SVar, !LPRInfo),
new_art_var(AVar, !LPRInfo),
Eqn = gte([AVar - one, SVar - (-one) | Coeffs], Const)
).
:- func negate_constraint(lp_constraint) = lp_constraint.
negate_constraint(lte(Terms, Const)) = gte(negate_lp_terms(Terms), -Const).
negate_constraint(eq(Terms, Const)) = eq(negate_lp_terms(Terms), -Const).
negate_constraint(gte(Terms, Const)) = lte(negate_lp_terms(Terms), -Const).
:- func negate_lp_terms(list(lp_term)) = list(lp_term).
negate_lp_terms(Terms) = assoc_list.map_values_only((func(X) = (-X)), Terms).
%---------------------------------------------------------------------------%
:- func collect_vars(lp_constraint_conj, objective) = set(lp_var).
collect_vars(Eqns, Obj) = Vars :-
GetVar =
( pred(Var::out) is nondet :-
(
list.member(Eqn, Eqns),
Coeffs = lp_terms(Eqn),
list.member(Pair, Coeffs)
;
list.member(Pair, Obj)
),
Var = fst(Pair)
),
solutions.solutions(GetVar, VarList),
Vars = set.list_to_set(VarList).
:- type var_num_map == map(lp_var, int).
:- func number_vars(list(lp_var), int) = var_num_map.
number_vars(Vars, N) = VarNum :-
number_vars_2(Vars, N, map.init, VarNum).
:- pred number_vars_2(list(lp_var)::in, int::in,
var_num_map::in, var_num_map::out) is det.
number_vars_2([], _, !VarNums).
number_vars_2([Var | Vars], N, !VarNums) :-
map.det_insert(Var, N, !VarNums),
number_vars_2(Vars, N + 1, !VarNums).
:- pred insert_constraints(lp_constraint_conj::in, int::in, int::in,
var_num_map::in, tableau::in, tableau::out) is det.
insert_constraints([], _, _, _, !Tableau).
insert_constraints([C | Cs], Row, ConstCol, VarNums, !Tableau) :-
insert_terms(lp_terms(C), Row, VarNums, !Tableau),
set_cell(Row, ConstCol, constant(C), !Tableau),
insert_constraints(Cs, Row + 1, ConstCol, VarNums, !Tableau).
:- pred insert_terms(list(lp_term)::in, int::in, var_num_map::in,
tableau::in, tableau::out) is det.
insert_terms([], _, _, !Tableau).
insert_terms([Var - Const | Coeffs], Row, VarNums, !Tableau) :-
map.lookup(VarNums, Var, Col),
set_cell(Row, Col, Const, !Tableau),
insert_terms(Coeffs, Row, VarNums, !Tableau).
%---------------------------------------------------------------------------%
:- pred optimize(list(lp_var)::in, lp_result::out, tableau::in, tableau::out)
is det.
optimize(ObjVars, Result, !Tableau) :-
simplex(Result0, !Tableau),
(
Result0 = no,
Result = lp_res_unbounded
;
Result0 = yes,
ObjVal = !.Tableau ^ elem(0, !.Tableau ^ cols),
ObjMap = extract_objective(ObjVars, !.Tableau),
Result = lp_res_satisfiable(ObjVal, ObjMap)
).
:- func extract_objective(list(lp_var), tableau) = map(lp_var, rat).
extract_objective(ObjVars, Tableau) = Objective :-
Objective = list.foldl(extract_obj_var(Tableau), ObjVars, map.init).
:- func extract_obj_var(tableau, lp_var, map(lp_var, rat))
= map(lp_var, rat).
extract_obj_var(Tableau, Var, Map0) = Map :-
extract_obj_var2(Tableau, Var, Val),
map.set(Var, Val, Map0, Map).
:- pred extract_obj_var2(tableau::in, lp_var::in, rat::out) is det.
extract_obj_var2(Tableau, Var, Val) :-
Col = var_col(Tableau, Var),
GetCell =
( pred(Val0::out) is nondet :-
all_rows(Tableau, Row),
one = Tableau ^ elem(Row, Col),
Val0 = Tableau ^ elem(Row, Tableau ^ cols)
),
solutions.solutions(GetCell, Solns),
( if Solns = [Val1] then Val = Val1 else Val = zero ).
:- pred simplex(bool::out, tableau::in, tableau::out) is det.
simplex(Result, !Tableau) :-
AllColumns = all_cols(!.Tableau),
MinAgg =
( pred(Col::in, !.Min::in, !:Min::out) is det :-
(
!.Min = no,
MinVal = !.Tableau ^ elem(0, Col),
( if MinVal < zero then
!:Min = yes(Col - MinVal)
else
!:Min = no
)
;
!.Min = yes(_ - MinVal0),
CellVal = !.Tableau ^ elem(0, Col),
( if CellVal < MinVal0 then
!:Min = yes(Col - CellVal)
else
true
)
)
),
solutions.aggregate(AllColumns, MinAgg, no, MinResult),
(
MinResult = no,
Result = yes
;
MinResult = yes(Q - _Val),
AllRows = all_rows(!.Tableau),
MaxAgg =
( pred(Row::in, !.Max::in, !:Max::out) is det :-
(
!.Max = no,
MaxVal = !.Tableau ^ elem(Row, Q),
( if MaxVal > zero then
Col = !.Tableau ^ cols,
MVal = !.Tableau ^ elem(Row, Col),
( if MaxVal = zero then
unexpected($pred, "zero divisor")
else
true
),
CVal = MVal / MaxVal,
!:Max = yes(Row - CVal)
else
!:Max = no
)
;
!.Max = yes(_ - MaxVal0),
CellVal = !.Tableau ^ elem(Row, Q),
RHSC = rhs_col(!.Tableau),
MVal = !.Tableau ^ elem(Row, RHSC),
( if CellVal =< zero then
% CellVal = 0 => multiple optimal sol'ns.
true
else
( if CellVal = zero then
unexpected($pred, "zero divisor")
else
true
),
MaxVal1 = MVal / CellVal,
( if MaxVal1 =< MaxVal0 then
!:Max = yes(Row - MaxVal1)
else
true
)
)
)
),
solutions.aggregate(AllRows, MaxAgg, no, MaxResult),
(
MaxResult = no,
Result = no
;
MaxResult = yes(P - _),
pivot(P, Q, !Tableau),
disable_warning [suspicious_recursion] (
simplex(Result, !Tableau)
)
)
).
%---------------------------------------------------------------------------%
:- pred ensure_zero_obj_coeffs(list(lp_var)::in,
tableau::in, tableau::out) is det.
ensure_zero_obj_coeffs([], !Tableau).
ensure_zero_obj_coeffs([Var | Vars], !Tableau) :-
Col = var_col(!.Tableau, Var),
Val = !.Tableau ^ elem(0, Col),
( if Val = zero then
ensure_zero_obj_coeffs(Vars, !Tableau)
else
FindOne =
( pred(P::out) is nondet :-
all_rows(!.Tableau, R),
ValF0 = !.Tableau ^ elem(R, Col),
ValF0 \= zero,
P = R - ValF0
),
solutions.solutions(FindOne, Ones),
(
Ones = [Row - Fac0 | _],
( if Fac0 = zero then
unexpected($pred, "zero divisor")
else
true
),
Fac = -Val / Fac0,
row_op(Fac, Row, 0, !Tableau),
ensure_zero_obj_coeffs(Vars, !Tableau)
;
Ones = [],
unexpected($pred, "problem with artificial variable")
)
).
:- pred fix_basis_and_rem_cols(list(lp_var)::in,
tableau::in, tableau::out) is det.
fix_basis_and_rem_cols([], !Tableau).
fix_basis_and_rem_cols([Var | Vars], !Tableau) :-
Col = var_col(!.Tableau, Var),
BasisAgg =
( pred(R::in, Ones0::in, Ones::out) is det :-
Val = !.Tableau ^ elem(R, Col),
Ones = ( if Val = zero then Ones0 else [Val - R | Ones0] )
),
solutions.aggregate(all_rows(!.Tableau), BasisAgg, [], Res),
( if Res = [one - Row] then
PivGoal =
( pred(Col1::out) is nondet :-
all_cols(!.Tableau, Col1),
Col \= Col1,
Zz = !.Tableau ^ elem(Row, Col1),
Zz \= zero
),
solutions.solutions(PivGoal, PivSolns),
(
PivSolns = [],
remove_col(Col, !Tableau),
remove_row(Row, !Tableau)
;
PivSolns = [Col2 | _],
pivot(Row, Col2, !Tableau),
remove_col(Col, !Tableau)
)
else
true
),
remove_col(Col, !Tableau),
fix_basis_and_rem_cols(Vars, !Tableau).
%---------------------------------------------------------------------------%
:- type cell
---> cell(int, int).
:- pred pivot(int::in, int::in, tableau::in, tableau::out) is det.
pivot(P, Q, !Tableau) :-
Apq = !.Tableau ^ elem(P, Q),
MostCells =
( pred(Cell::out) is nondet :-
all_rows0(!.Tableau, J),
J \= P,
all_cols0(!.Tableau, K),
K \= Q,
Cell = cell(J, K)
),
ScaleCell =
( pred(Cell::in, T0::in, T::out) is det :-
Cell = cell(J, K),
Ajk = T0 ^ elem(J, K),
Ajq = T0 ^ elem(J, Q),
Apk = T0 ^ elem(P, K),
( if Apq = zero then
unexpected($pred, "ScaleCell: zero divisor")
else
true
),
T = T0 ^ elem(J, K) := Ajk - Apk * Ajq / Apq
),
solutions.aggregate(MostCells, ScaleCell, !Tableau),
QColumn =
( pred(Cell::out) is nondet :-
all_rows0(!.Tableau, J),
Cell = cell(J, Q)
),
Zero =
( pred(Cell::in, T0::in, T::out) is det :-
Cell = cell(J, K),
T = T0 ^ elem(J, K) := zero
),
solutions.aggregate(QColumn, Zero, !Tableau),
PRow = all_cols0(!.Tableau),
ScaleRow =
( pred(K::in, T0::in, T::out) is det :-
Apk = T0 ^ elem(P, K),
( if Apq = zero then
unexpected($pred, "ScaleRow: zero divisor")
else
true
),
T = T0 ^ elem(P, K) := Apk / Apq
),
solutions.aggregate(PRow, ScaleRow, !Tableau),
set_cell(P, Q, one, !Tableau).
:- pred row_op(rat::in, int::in, int::in, tableau::in,
tableau::out) is det.
row_op(Scale, From, To, !Tableau) :-
AllCols = all_cols0(!.Tableau),
AddRow =
( pred(Col::in, T0::in, T::out) is det :-
X = T0 ^ elem(From, Col),
Y = T0 ^ elem(To, Col),
Z = Y + (Scale * X),
T = T0 ^ elem(To, Col) := Z
),
solutions.aggregate(AllCols, AddRow, !Tableau).
%---------------------------------------------------------------------------%
% XXX We should try using arrays or version_arrays for the simplex tableau.
% (We should try this in lp.m as well).
:- type tableau
---> tableau(
rows :: int,
cols :: int,
var_nums :: map(lp_var, int),
shunned_rows :: list(int),
shunned_cols :: list(int),
cells :: map(pair(int), rat)
).
:- func init_tableau(int, int, map(lp_var, int)) = tableau.
init_tableau(Rows, Cols, VarNums) = Tableau :-
Tableau = tableau(Rows, Cols, VarNums, [], [], map.init).
:- func tableau ^ elem(int, int) = rat.
Tableau ^ elem(Row, Col) = get_cell(Tableau, Row, Col).
:- func tableau ^ elem(int, int) := rat = tableau.
Tableau0 ^ elem(Row, Col) := Cell = Tableau :-
set_cell(Row, Col, Cell, Tableau0, Tableau).
:- func get_cell(tableau, int, int) = rat.
get_cell(Tableau, Row, Col) = Cell :-
( if
( list.member(Row, Tableau ^ shunned_rows)
; list.member(Col, Tableau ^ shunned_cols)
)
then
unexpected($pred, "attempt to address shunned row/col")
else
true
),
( if map.search(Tableau ^ cells, Row - Col, Cell0) then
Cell = Cell0
else
Cell = zero
).
:- pred set_cell(int::in, int::in, rat::in, tableau::in,
tableau::out) is det.
set_cell(J, K, R, Tableau0, Tableau) :-
Tableau0 = tableau(Rows, Cols, VarNums, SR, SC, Cells0),
( if
( list.member(J, SR)
; list.member(K, SC)
)
then
unexpected($pred, "Attempt to write shunned row/col")
else
true
),
( if R = zero then
Cells = map.delete(Cells0, J - K)
else
Cells = map.set(Cells0, J - K, R)
),
Tableau = tableau(Rows, Cols, VarNums, SR, SC, Cells).
% Returns the number of the RHS column in the tableau.
%
:- func rhs_col(tableau) = int.
rhs_col(Tableau) = Tableau ^ cols.
:- pred all_rows0(tableau::in, int::out) is nondet.
all_rows0(Tableau, Row) :-
between(0, Tableau ^ rows, Row),
not list.member(Row, Tableau ^ shunned_rows).
:- pred all_rows(tableau::in, int::out) is nondet.
all_rows(Tableau, Row) :-
between(1, Tableau ^ rows, Row),
not list.member(Row, Tableau ^ shunned_rows).
:- pred all_cols0(tableau::in, int::out) is nondet.
all_cols0(Tableau, Col) :-
between(0, Tableau ^ cols, Col),
not list.member(Col, Tableau ^ shunned_cols).
:- pred all_cols(tableau::in, int::out) is nondet.
all_cols(Tableau, Col) :-
Cols1 = Tableau ^ cols - 1,
between(0, Cols1, Col),
not list.member(Col, Tableau ^ shunned_cols).
:- func var_col(tableau, lp_var) = int.
var_col(Tableau, Var) = (Tableau ^ var_nums) ^ det_elem(Var).
:- pred remove_row(int::in, tableau::in, tableau::out) is det.
remove_row(Row, !Tableau) :-
SR = !.Tableau ^ shunned_rows,
!Tableau ^ shunned_rows := [Row | SR].
:- pred remove_col(int::in, tableau::in, tableau::out) is det.
remove_col(C, Tableau0, Tableau) :-
Tableau0 = tableau(Rows, Cols, VarNums, SR, SC, Cells),
Tableau = tableau(Rows, Cols, VarNums, SR, [C | SC], Cells).
:- func get_basis_vars(tableau) = list(lp_var).
get_basis_vars(Tableau) = Vars :-
BasisCol =
( pred(C::out) is nondet :-
all_cols(Tableau, C),
NonZeroGoal =
( pred(P::out) is nondet :-
all_rows(Tableau, R),
Z = Tableau ^ elem(R, C),
Z \= zero,
P = R - Z
),
solutions.solutions(NonZeroGoal, Solns),
Solns = [_ - one]
),
solutions.solutions(BasisCol, Cols),
BasisVars =
( pred(V::out) is nondet :-
list.member(Col, Cols),
map.member(Tableau ^ var_nums, V, Col)
),
solutions.solutions(BasisVars, Vars).
%---------------------------------------------------------------------------%
:- func lpr_info_init(lp_varset) = lpr_info.
lpr_info_init(VarSet) = lpr_info(VarSet, [], []).
:- pred new_slack_var(lp_var::out, lpr_info::in, lpr_info::out) is det.
new_slack_var(Var, !LPRInfo) :-
varset.new_var(Var, !.LPRInfo ^ lpr_varset, VarSet),
!LPRInfo ^ lpr_varset := VarSet,
Vars = !.LPRInfo ^ lpr_slack_vars,
!LPRInfo ^ lpr_slack_vars := [Var | Vars].
:- pred new_art_var(lp_var::out, lpr_info::in, lpr_info::out) is det.
new_art_var(Var, !LPRInfo) :-
varset.new_var(Var, !.LPRInfo ^ lpr_varset, VarSet),
!LPRInfo ^ lpr_varset := VarSet,
Vars = !.LPRInfo ^ lpr_art_vars,
!LPRInfo ^ lpr_art_vars := [Var | Vars].
%---------------------------------------------------------------------------%
:- pred between(int::in, int::in, int::out) is nondet.
between(Min, Max, I) :-
Min =< Max,
(
I = Min
;
between(Min + 1, Max, I)
).
%---------------------------------------------------------------------------%
%---------------------------------------------------------------------------%
%
% Projection.
%
%
% The following code more or less follows the algorithm described in:
% Joxan Jaffar, Michael Maher, Peter Stuckey and Roland Yap.
% Projecting CLP(R) Constraints. New Generation Computing 11(3): 449-469.
%
% * Linear equations (Gaussian elimination)
% - substitutions need to be performed on the inequalities as well.
% * Linear inequalities (Fourier elimination)
%
% We next convert any remaining equations into opposing inequalities and
% then use Fourier elimination to try and eliminate any remaining target
% variables. The main problem here is ensuring that we don't get
% swamped by redundant constraints.
%
% The implementation here uses the extensions to FM elimination described by
% Cernikov as well as some other redundancy checks. Note that in general
% arbitrarily mixing redundancy elimination techniques with the Cernikov
% methods is unsound (See the above article for an example).
%
% In addition to Cernikov's methods and quasi-syntactic redundancy checks
% we also use a heuristic developed by Duffin to choose the order in
% which we eliminate variables (See below).
%
%---------------------------------------------------------------------------%
:- type vector
---> vector(
% The vector's label is for redundancy checking
% during Fourier elimination - see below.
label :: set(int),
% A map from each variable in the vector to its coefficient.
terms :: map(lp_var, lp_coefficient),
const :: lp_constant
).
:- type matrix == list(vector).
project_constraints(VarSet, Vars, Constraints, Result) :-
project_constraints_maybe_size_limit(VarSet, no, Vars,
Constraints, Result).
% For the first clause, the `Constraints' may actually be
% an inconsistent system - we don't bother checking that here though.
% We instead delay that until we need to perform an entailment check.
project_constraints_maybe_size_limit(_, _, [],
Constraints, pr_res_ok(Constraints)).
project_constraints_maybe_size_limit(VarSet, MaybeThreshold,
!.Vars @ [_ | _], Constraints0, Result) :-
eliminate_equations(!Vars, Constraints0, EqlResult),
(
EqlResult = pr_res_inconsistent,
Result = pr_res_inconsistent
;
% Elimination of equations should not cause an abort, since we always
% make the matrix smaller.
EqlResult = pr_res_aborted,
unexpected($pred, "abort from eliminate_equations")
;
EqlResult = pr_res_ok(Constraints1),
% Skip the call to fourier_elimination/6 if there are no variables to
% project - this avoids the transformation to vector form.
(
!.Vars = [_ | _],
Matrix0 = constraints_to_matrix(Constraints1),
fourier_elimination(!.Vars, VarSet, MaybeThreshold, 0,
Matrix0, FourierResult),
(
FourierResult = yes(Matrix),
Constraints = matrix_to_constraints(Matrix),
Result = pr_res_ok(Constraints)
;
FourierResult = no,
Result = pr_res_aborted
)
;
% NOTE: the matrix `Constraints1' may actually be inconsistent here
% - we don't bother checking at this point because that would mean
% traversing the matrix, so we wait until the next operation that
% needs to traverse it anyway or until the next entailment check.
!.Vars = [],
Result = pr_res_ok(Constraints1)
)
).
%---------------------------------------------------------------------------%
%
% Convert each constraint into `=<' form and give each an initial label.
%
:- func constraints_to_matrix(lp_constraint_conj) = matrix.
constraints_to_matrix(Constraints) = Matrix :-
list.foldl2(fm_standardize, Constraints, 0, _, [], Matrix).
:- pred fm_standardize(lp_constraint::in, int::in, int::out, matrix::in,
matrix::out) is det.
fm_standardize(lte(Terms0, Constant), !Labels, !Matrix) :-
Terms = lp_terms_to_map(Terms0),
make_label(Label, !Labels),
list.cons(vector(Label, Terms, Constant), !Matrix).
fm_standardize(eq(Terms, Constant), !Labels, !Matrix) :-
make_label(Label1, !Labels),
make_label(Label2, !Labels),
Vector1 = vector(Label1, lp_terms_to_map(Terms), Constant),
Vector2 = vector(Label2, lp_terms_to_map(negate_lp_terms(Terms)),
-Constant),
list.append([Vector1, Vector2], !Matrix).
fm_standardize(gte(Terms0, Constant), !Labels, !Matrix) :-
make_label(Label, !Labels),
Terms = lp_terms_to_map(negate_lp_terms(Terms0)),
list.cons(vector(Label, Terms, -Constant), !Matrix).
:- pred make_label(set(int)::out, int::in, int::out) is det.
make_label(Label, Labels, Labels + 1) :-
Label = set.make_singleton_set(Labels).
:- func matrix_to_constraints(matrix) = lp_constraint_conj.
matrix_to_constraints(Matrix) = list.map(vector_to_constraint, Matrix).
:- func vector_to_constraint(vector) = lp_constraint.
vector_to_constraint(vector(_, Terms0, Constant0)) = Constraint :-
Terms1 = map.to_assoc_list(Terms0),
normalize_terms_and_const(yes, Terms1, Constant0, Terms, Constant),
Constraint = lte(Terms, Constant).
%---------------------------------------------------------------------------%
%
% Predicates for eliminating equations from the constraints.
% (Gaussian elimination)
%
% Split the constraints into a set of inequalities and a set of equalities.
% For every variable in the set of target variables (i.e. those we are
% eliminating), check if there is at least one equality that contains
% that variable. If so, then substitute the value of that variable
% into the other constraints. Return the set of target variables
% that do not occur in any equality.
%
:- pred eliminate_equations(list(lp_var)::in, list(lp_var)::out,
lp_constraint_conj::in, projection_result::out) is det.
eliminate_equations(!Vars, Constraints0, Result) :-
Constraints = simplify_constraints(Constraints0),
list.filter((pred(eq(_, _)::in) is semidet), Constraints,
Equalities0, Inequalities0),
( if
eliminate_equations_2(!Vars, Equalities0, Equalities,
Inequalities0, Inequalities)
then
Result = pr_res_ok(Equalities ++ Inequalities)
else
Result = pr_res_inconsistent
).
:- pred eliminate_equations_2(list(lp_var)::in, list(lp_var)::out,
lp_constraint_conj::in, lp_constraint_conj::out, lp_constraint_conj::in,
lp_constraint_conj::out) is semidet.
eliminate_equations_2([], [], !Equations, !Inequations).
eliminate_equations_2([Var | !.Vars], !:Vars, !Equations, !Inequations) :-
eliminate_equations_2(!Vars, !Equations, !Inequations),
( if find_target_equality(Var, Target, !Equations) then
substitute_variable(Target, Var, !Equations, !Inequations, Succeeded),
(
Succeeded = no,
list.cons(Var, !Vars),
list.cons(Target, !Equations)
;
Succeeded = yes
)
else
list.cons(Var, !Vars)
).
% Find an equation that constrains a variable we are trying to eliminate.
%
:- pred find_target_equality(lp_var::in, lp_constraint::out,
lp_constraint_conj::in, lp_constraint_conj::out) is semidet.
find_target_equality(Var, Target, Constraints0, Constraints) :-
Result = find_target_equality(Var, Constraints0),
Result = yes(Target - Constraints).
:- func find_target_equality(lp_var, lp_constraint_conj) =
maybe(pair(lp_constraint, lp_constraint_conj)).
find_target_equality(Var, Eqns) = find_target_equality_2(Var, Eqns, []).
:- func find_target_equality_2(lp_var, lp_constraint_conj, lp_constraint_conj)
= maybe(pair(lp_constraint, lp_constraint_conj)).
find_target_equality_2(_, [], _) = no.
find_target_equality_2(Var, [Eqn | Eqns], Acc) = MaybeTargetEqn :-
( if operator(Eqn) = lp_eq then
true
else
unexpected($pred, "inequality encountered")
),
Coeffs = lp_terms(Eqn),
( if list.member(Var - _, Coeffs) then
MaybeTargetEqn = yes(Eqn - (Acc ++ Eqns))
else
MaybeTargetEqn = find_target_equality_2(Var, Eqns, [Eqn | Acc])
).
% Given a target equation of the form a1x1 + .. + aNxN = C and
% a target variable, say `x1', notionally rewrite the equation as:
%
% x1 = C - ... aN/a1 xN
%
% and then substitute that value for x1 in the supplied sets
% of equations and inequations.
%
:- pred substitute_variable(lp_constraint::in, lp_var::in,
lp_constraint_conj::in, lp_constraint_conj::out, lp_constraint_conj::in,
lp_constraint_conj::out,
bool::out) is semidet.
substitute_variable(Target0, Var, !Equations, !Inequations, Flag) :-
normalize_constraint(Var, Target0, Target),
deconstruct_constraint(Target, TargetCoeffs, Op, TargetConst),
expect(unify(Op, lp_eq), $pred, "inequality encountered"),
fix_coeff_and_const(Var, TargetCoeffs, TargetConst, Coeffs, Const),
substitute_into_constraints(Var, Coeffs, Const, !Equations, EqlFlag),
substitute_into_constraints(Var, Coeffs, Const, !Inequations, IneqlFlag),
Flag = bool.or(EqlFlag, IneqlFlag).
% Multiply the terms and constant except for the term containing
% the specified variable in preparation for making a substitution
% for that variable. Notionally this converts a constraint of the form:
% t + z + w = C ... C is a constant
%
% into:
%
% t = C - z - w
%
:- pred fix_coeff_and_const(lp_var::in, list(lp_term)::in, lp_constant::in,
list(lp_term)::out, lp_constant::out) is det.
fix_coeff_and_const(_, [], Const, [], -Const).
fix_coeff_and_const(Var, [Var1 - Coeff1 | Coeffs], Const0, FixedCoeffs,
Const) :-
fix_coeff_and_const(Var, Coeffs, Const0, FCoeffs0, Const),
( if Var = Var1 then
FixedCoeffs = FCoeffs0
else
FixedCoeffs = [Var1 - (-Coeff1) | FCoeffs0]
).
% The `Flag' argument is `yes' if one or more substitutions were made,
% `no' otherwise. substitute_into_constraints/7 fails if a false constraint
% is generated as a result of a substitution. This means that the original
% matrix was inconsistent.
%
:- pred substitute_into_constraints(lp_var::in, list(lp_term)::in,
lp_constant::in, lp_constraint_conj::in, lp_constraint_conj::out,
bool::out) is semidet.
substitute_into_constraints(_, _, _, [], [], no).
substitute_into_constraints(Var, Coeffs, Const, [Constr0 | Constrs0], Result,
Flag) :-
substitute_into_constraint(Var, Coeffs, Const, Constr0, Constr, Flag0),
not is_false(Constr),
substitute_into_constraints(Var, Coeffs, Const, Constrs0, Constrs, Flag1),
Result = ( if is_true(Constr) then Constrs else [Constr | Constrs] ),
Flag = bool.or(Flag0, Flag1).
:- pred substitute_into_constraint(lp_var::in, list(lp_term)::in,
lp_constant::in, lp_constraint::in, lp_constraint::out, bool::out) is det.
substitute_into_constraint(Var, SubCoeffs, SubConst, !Constraint, Flag) :-
normalize_constraint(Var, !Constraint),
deconstruct_constraint(!.Constraint, TargetCoeffs, Op, TargetConst),
( if list.member(Var - one, TargetCoeffs) then
FinalCoeffs0 = lp_terms_to_map(TargetCoeffs ++ SubCoeffs),
% Delete the target variable from both constraints.
FinalCoeffs1 = map.delete(FinalCoeffs0, Var),
FinalCoeffs = map.to_assoc_list(FinalCoeffs1),
FinalConst = TargetConst + SubConst,
!:Constraint = construct_constraint(FinalCoeffs, Op, FinalConst),
Flag = yes
else
Flag = no
).
%---------------------------------------------------------------------------%
%
% Fourier elimination.
%
% Will return `no' if it aborts otherwise `yes(Matrix)', where
% `Matrix' is the result of the projection.
%
:- pred fourier_elimination(list(lp_var)::in, lp_varset::in, maybe(int)::in,
int::in, matrix::in, maybe(matrix)::out) is det.
fourier_elimination([], _, _, _, Matrix, yes(Matrix)).
fourier_elimination(Vars @ [Var0 | Vars0], VarSet, MaybeThreshold, !.Step,
Matrix0, Result) :-
% Use Duffin's heuristic to try and find a "nice" variable to eliminate.
%
% NOTE: the heuristic will fail if none of the variables being projected
% actually occur in the constraints. In that case, we just pick
% the first one - it doesn't really matter since the projection
% will be trivial.
( if duffin_heuristic(Vars, Matrix0, TargetVar0, OtherVars0) then
Var = TargetVar0,
OtherVars = OtherVars0
else
Var = Var0,
OtherVars = Vars0
),
separate_vectors(Matrix0, Var, PosMatrix, NegMatrix, ZeroMatrix,
SizeZeroMatrix),
% `Step' counts active Fourier eliminations only. An elimination is active
% if at least one constraint contains a term that has a non-zero
% coefficient for the variable being eliminated.
( if
PosMatrix = [_ | _],
NegMatrix = [_ | _]
then
!:Step = !.Step + 1,
( if
list.foldl2(eliminate_var(!.Step, MaybeThreshold, NegMatrix),
PosMatrix, ZeroMatrix, ResultMatrix, SizeZeroMatrix, _)
then
NewMatrix = yes(ResultMatrix)
else
NewMatrix = no
)
else
NewMatrix = yes(ZeroMatrix)
),
(
NewMatrix = yes(Matrix),
fourier_elimination(OtherVars, VarSet, MaybeThreshold, !.Step,
Matrix, Result)
;
NewMatrix = no,
Result = no
).
% separate_vectors(Matrix, Var, Positive, Negative, Zero, Num).
% `Positive' is a matrix containing those constraints of `Matrix' for
% which the coefficient of `Var' is positive. `Negative' similarly
% for those which the coefficient of `Var' is negative and `Zero'
% those for which the coefficient of `Var' is zero. `Num' is the
% number of constraints in `Zero'.
%
:- pred separate_vectors(matrix::in, lp_var::in, matrix::out, matrix::out,
matrix::out, int::out) is det.
separate_vectors(Matrix, Var, Pos, Neg, Zero, NumZeros) :-
list.foldl4(classify_vector(Var), Matrix, [], Pos, [], Neg, [], Zero,
0, NumZeros).
:- pred classify_vector(lp_var::in, vector::in, matrix::in,
matrix::out, matrix::in, matrix::out, matrix::in, matrix::out,
int::in, int::out) is det.
classify_vector(Var, Vector0, !Pos, !Neg, !Zero, !Num) :-
( if map.search(Vector0 ^ terms, Var, Coefficient) then
Vector0 = vector(Label, Terms0, Const0),
normalize_vector(Var, Terms0, Terms, Const0, Const),
Vector1 = vector(Label, Terms, Const),
( if Coefficient > zero then
list.cons(Vector1, !Pos)
else
list.cons(Vector1, !Neg)
)
else
list.cons(Vector0, !Zero),
!:Num = !.Num + 1
).
:- pred eliminate_var(int::in, maybe(int)::in, matrix::in,
vector::in, matrix::in, matrix::out, int::in, int::out) is semidet.
eliminate_var(Step, MaybeThreshold, NegMatrix, PosVector, !Zeros,
!ZerosSize) :-
list.foldl2(combine_vectors(Step, MaybeThreshold, PosVector),
NegMatrix, !Zeros, !ZerosSize).
:- pred combine_vectors(int::in, maybe(int)::in, vector::in,
vector::in, matrix::in, matrix::out, int::in, int::out) is semidet.
combine_vectors(Step, MaybeThreshold, vector(LabelPos, TermsPos, ConstPos),
vector(LabelNeg, TermsNeg, ConstNeg), !Zeros, !Num) :-
LabelNew = set.union(LabelPos, LabelNeg),
( if
% If the cardinality of the label set is greater than `Step + 2'
% then the constraint we are trying to add is redundant.
set.count(LabelNew) < Step + 2
then
add_vectors(TermsPos, ConstPos, TermsNeg, ConstNeg, Coeffs, Const),
New = vector(LabelNew, Coeffs, Const),
( if
(
% Do not bother adding the new constraint
% if it is just `true'.
map.is_empty(Coeffs),
Const >= zero
;
list.member(Vec, !.Zeros),
quasi_syntactic_redundant(New, Vec)
)
then
% If the new constraint is `true' or is quasi-syntactic redundant
% with something already there.
true
else
% Remove anything in the matrix that is quasi-syntactic redundant
% w.r.t the new constraint.
filter_and_count(
( pred(Vec2::in) is semidet :-
not quasi_syntactic_redundant(Vec2, New)
),
!.Zeros, [], !:Zeros, 0, !:Num),
( if
list.member(Vec, !.Zeros),
label_subsumed(New, Vec)
then
% Do not add the new constraint because it is label subsumed
% by something already in the matrix.
true
else
filter_and_count(
( pred(Vec2::in) is semidet :-
not label_subsumed(Vec2, New)
),
!.Zeros, [], !:Zeros, 0, !:Num),
list.cons(New, !Zeros),
!:Num = !.Num + 1
)
)
else
true
),
% Check that the size of the matrix does not exceed the threshold
% for aborting the projection.
not (
MaybeThreshold = yes(Threshold),
!.Num > Threshold
).
%---------------------------------------------------------------------------%
:- pred filter_and_count(pred(vector)::in(pred(in) is semidet),
matrix::in, matrix::in, matrix::out, int::in, int::out) is det.
filter_and_count(_, [], !Acc, !Count).
filter_and_count(P, [X | Xs], !Acc, !Count) :-
( if P(X) then
list.cons(X, !Acc),
!:Count = !.Count + 1
else
true
),
filter_and_count(P, Xs, !Acc, !Count).
%---------------------------------------------------------------------------%
%
% Detection of quasi-syntactic redundancy.
%
% Succeeds if the first vector is quasi-syntactic redundant wrt to the
% second. That is c = c' + (0 < e), for e > 0.
%
:- pred quasi_syntactic_redundant(vector::in, vector::in) is semidet.
quasi_syntactic_redundant(VecA, VecB) :-
VecB ^ const < VecA ^ const,
all [Var] (
map.member(VecA ^ terms, Var, Coeff)
<=>
map.member(VecB ^ terms, Var, Coeff)
).
%---------------------------------------------------------------------------%
%
% Label subsumption.
%
% label_subsumed(A, B):
%
% Succeeds iff constraint A is label subsumed by constraint B.
%
:- pred label_subsumed(vector::in, vector::in) is semidet.
label_subsumed(VectorA, VectorB) :-
set.subset(VectorB ^ label, VectorA ^ label).
%---------------------------------------------------------------------------%
%
% Duffin's heuristic.
%
%
% This attempts to find an order in which to eliminate variables such that
% the minimal number of redundant constraints are generated at each
% Fourier step. For each variable, x_h, to be eliminated, we
% calculate E(x_h) which is defined as follows:
%
% E(x_h) = p(x_h)q(x_h) + r(x_h) ... if p(x_h) + q(x_h) > 0
% E(x_h) = 0 ... if p(x_h) + q(x_h) = 0
%
% p, q, r are the number of positive, negative and zero coefficients
% of the variable x_h respectively in the system of constraints under
% consideration. E(x_h) is called the expansion number of x_h.
%
% We eliminate the variable that has minimal expansion number.
%
% For further details see:
% R.J. Duffin. On Fourier's Analysis of Linear Inequality Systems.
% Mathematical Programming Study 1, 71 - 95 (1974).
%
%---------------------------------------------------------------------------%
% We only count the occurrences of positive and negative coefficients.
% We can work out the zero occurrences by subtracting the two
% previous totals from the total number of constraints.
%
:- type coeff_info
---> coeff_info(
pos :: int,
neg :: int
).
:- type cc_map == map(lp_var, coeff_info).
% Calculates the variable with the minimal expansion number and
% returns that variable. (Removes those variables that have an
% expansion number of zero, because there are no constraints on them
% anyway). Fails if it can't find such a variable, ie. none of the
% variables being eliminated actually occurs in the constraints.
%
:- pred duffin_heuristic(list(lp_var)::in, matrix::in, lp_var::out,
list(lp_var)::out) is semidet.
duffin_heuristic([Var], _, Var, []).
duffin_heuristic(Vars0 @ [_, _ | _], Matrix, TargetVar, Vars) :-
VarsAndNums0 = generate_expansion_nums(Vars0, Matrix),
VarsAndNums1 = list.filter(relevant, VarsAndNums0),
VarsAndNums1 \= [],
TargetVar = find_max(VarsAndNums1),
Vars = collect_remaining_vars(VarsAndNums1, TargetVar).
:- func collect_remaining_vars(assoc_list(lp_var, int), lp_var) = list(lp_var).
collect_remaining_vars([], _) = [].
collect_remaining_vars([Var - _ | Rest], TargetVar) = Result :-
( if Var = TargetVar then
Result = collect_remaining_vars(Rest, TargetVar)
else
Result = [Var | collect_remaining_vars(Rest, TargetVar)]
).
:- func find_max(list(pair(lp_var, int))) = lp_var.
find_max([]) = unexpected($pred, "empty list").
find_max([Var0 - ExpnNum0 | Vars]) = fst(find_max_2(Vars, Var0 - ExpnNum0)).
:- func find_max_2(assoc_list(lp_var, int), pair(lp_var, int)) =
pair(lp_var, int).
find_max_2([], Best) = Best.
find_max_2([Var1 - ExpnNum1 | Vars], Var0 - ExpnNum0) =
( if ExpnNum1 < ExpnNum0 then
find_max_2(Vars, Var1 - ExpnNum1)
else
find_max_2(Vars, Var0 - ExpnNum0)
).
:- pred relevant(pair(lp_var, int)::in) is semidet.
relevant(Var) :-
Var \= _ - 0.
% Given a list of variables and a system of linear inequalities
% generate the expansion number for each of the variables in the list.
%
:- func generate_expansion_nums(list(lp_var), matrix)
= assoc_list(lp_var, int).
generate_expansion_nums(Vars0, Matrix) = ExpansionNums :-
Vars = list.sort_and_remove_dups(Vars0),
CoeffMap0 = init_cc_map(Vars),
CoeffMap = list.foldl(count_coeffs_in_vector, Matrix, CoeffMap0),
CoeffList = map.to_assoc_list(CoeffMap),
ConstrNum = list.length(Matrix),
ExpansionNums = list.map(make_expansion_num(ConstrNum), CoeffList).
:- func make_expansion_num(int, pair(lp_var, coeff_info)) = pair(lp_var, int).
make_expansion_num(ConstrNum, Var - coeff_info(Pos, Neg)) = Var - ExpnNum :-
PosAndNeg = Pos + Neg,
( if PosAndNeg = 0 then
ExpnNum = 0
else
ExpnNum = (Pos * Neg) + (ConstrNum - PosAndNeg)
).
:- func count_coeffs_in_vector(vector, cc_map) = cc_map.
count_coeffs_in_vector(Vector, Map0) = Map :-
CoeffList = map.to_assoc_list(Vector ^ terms),
list.foldl(count_coeff, CoeffList, Map0, Map).
:- pred count_coeff(lp_term::in, cc_map::in, cc_map::out) is det.
count_coeff(Var - Coeff, !Map) :-
( if map.search(!.Map, Var, coeff_info(Pos0, Neg0)) then
( if Coeff > zero then
Pos = Pos0 + 1,
Neg = Neg0
else if Coeff < zero then
Pos = Pos0,
Neg = Neg0 + 1
else
unexpected($pred, "zero coefficient")
),
map.det_update(Var, coeff_info(Pos, Neg), !Map)
else
true
% If the variable in the term was not in the map then it is not
% one of the ones that is being eliminated.
).
:- func init_cc_map(list(lp_var)) = cc_map.
init_cc_map(Vars) = list.foldl(InitMap, Vars, map.init) :-
InitMap = (func(Var, Map) =
map.det_insert(Map, Var, coeff_info(0, 0))
).
%---------------------------------------------------------------------------%
%
% Predicates for normalizing vectors and constraints.
%
% normalize_vector(Var, Terms0, Terms, Const0, Const):
%
% Multiply the given vector by a scalar appropriate to make the
% coefficient of the given variable in the vector one. Throws an exception
% if `Var' has a zero coefficient.
%
:- pred normalize_vector(lp_var::in,
map(lp_var, lp_coefficient)::in, map(lp_var, lp_coefficient)::out,
lp_constant::in, lp_constant::out) is det.
normalize_vector(Var, !Terms, !Constant) :-
( if map.search(!.Terms, Var, Coefficient) then
( if Coefficient = zero then
unexpected($pred, "zero coefficient in vector")
else
true
),
DivVal = rat.abs(Coefficient),
!:Terms = map.map_values_only((func(C) = C / DivVal), !.Terms),
!:Constant = !.Constant / DivVal
else
% In this case the coefficient of the variable was zero
% (implicit in the fact that it is not in the map).
true
).
% Multiply the given constraint by a scalar appropriate to make the
% coefficient of the given variable in the constraint one. If the variable
% does not occur in the constraint then the constraint is unchanged.
% If the constraint is an inequality the sign may be changed.
% Throws an exception if the variable is found in the constraint
% and it has a coefficient of zero.
%
:- pred normalize_constraint(lp_var::in, lp_constraint::in, lp_constraint::out)
is det.
normalize_constraint(Var, Constraint0, Constraint) :-
deconstruct_constraint(Constraint0, Terms0, Op0, Constant0),
( if assoc_list.search(Terms0, Var, Coefficient) then
( if Coefficient = zero then
unexpected($pred, "zero coefficient constraint")
else
true
),
Terms = list.map((func(V - C) = V - (C / Coefficient)), Terms0),
Constant = Constant0 / Coefficient,
Op = ( if Coefficient < zero then negate_operator(Op0) else Op0 )
else
% In this case the coefficient of the variable was zero
% (implicit in the fact that it is not in the list).
Terms = Terms0,
Op = Op0,
Constant = Constant0
),
Constraint = unchecked_construct_constraint(Terms, Op, Constant).
:- pred add_vectors(map(lp_var, lp_coefficient)::in, lp_constant::in,
map(lp_var, lp_coefficient)::in, lp_constant::in,
map(lp_var, lp_coefficient)::out, lp_constant::out) is det.
add_vectors(TermsA, ConstA, TermsB, ConstB, Terms, ConstA + ConstB) :-
IsMapKey =
( pred(Var::out) is nondet :-
map.member(TermsA, Var, _)
),
AddVal =
( pred(Var::in, Coeffs0::in, Coeffs::out) is det :-
map.lookup(TermsA, Var, NumA),
( if map.search(Coeffs0, Var, Num1) then
( if NumA + Num1 = zero then
Coeffs = map.delete(Coeffs0, Var)
else
Coeffs = map.det_update(Coeffs0, Var, NumA + Num1)
)
else
Coeffs = map.det_insert(Coeffs0, Var, NumA)
)
),
solutions.aggregate(IsMapKey, AddVal, TermsB, Terms).
%---------------------------------------------------------------------------%
%---------------------------------------------------------------------------%
%
% Entailment test.
%
entailed(VarSet, Constraints, lte(Objective, Constant)) = Result :-
SolverResult = lp_rational.solve(Constraints, max, Objective, VarSet),
(
SolverResult = lp_res_satisfiable(MaxVal, _),
Result = ( if MaxVal =< Constant then entailed else not_entailed )
;
SolverResult = lp_res_unbounded,
Result = not_entailed
;
SolverResult = lp_res_inconsistent,
Result = inconsistent
).
entailed(VarSet, Constraints, eq(Objective, Constant)) = Result :-
Result0 = entailed(VarSet, Constraints, lte(Objective, Constant)),
(
Result0 = entailed,
Result = entailed(VarSet, Constraints, gte(Objective, Constant))
;
( Result0 = not_entailed
; Result0 = inconsistent
),
Result0 = Result
).
entailed(VarSet, Constraints, gte(Objective, Constant)) = Result :-
SolverResult = lp_rational.solve(Constraints, min, Objective, VarSet),
(
SolverResult = lp_res_satisfiable(MinVal, _),
Result = ( if MinVal >= Constant then entailed else not_entailed )
;
SolverResult = lp_res_unbounded,
Result = not_entailed
;
SolverResult = lp_res_inconsistent,
Result = inconsistent
).
entailed(VarSet, Constraints, Constraint) :-
Result = entailed(VarSet, Constraints, Constraint),
(
Result = entailed
;
Result = inconsistent,
unexpected($pred, "inconsistent constraint set")
;
Result = not_entailed,
fail
).
%---------------------------------------------------------------------------%
%
% Redundancy checking using the linear solver.
%
% Check if each constraint in the set is entailed by all the others.
% XXX It would be preferable not to use this as it can be very slow.
%
remove_some_entailed_constraints(VarSet, Constraints0, Constraints) :-
remove_some_entailed_constraints_2(VarSet, Constraints0, [], Constraints).
:- pred remove_some_entailed_constraints_2(lp_varset::in,
lp_constraint_conj::in, lp_constraint_conj::in,
lp_constraint_conj::out) is semidet.
remove_some_entailed_constraints_2(_, [], !Constraints).
remove_some_entailed_constraints_2(_, [ E ], !Constraints) :-
list.cons(E, !Constraints).
remove_some_entailed_constraints_2(VarSet, [E, X | Es], !Constraints) :-
( if obvious_constraint(E) then
true
else
RestOfMatrix = [X | Es] ++ !.Constraints,
Result = entailed(VarSet, RestOfMatrix, E),
(
Result = entailed
;
Result = not_entailed,
list.cons(E, !Constraints)
;
Result = inconsistent,
fail
)
),
remove_some_entailed_constraints_2(VarSet, [X | Es], !Constraints).
%---------------------------------------------------------------------------%
%
% Intermodule optimization stuff.
%
% The following predicates write out constraints in a form that is useful
% for (transitive) intermodule optimization.
% XXX This should not be needed; (transitive) intermodule optimization
% should output these constraints only as parts of termination pragmas,
% and that should be done by parse_tree_out_pragma.m.
output_constraints(Stream, OutputVar, Constraints, !IO) :-
io.write_char(Stream, '[', !IO),
write_out_list(output_constraint(OutputVar), ", ", Constraints,
Stream, !IO),
io.write_char(Stream, ']', !IO).
:- pred output_constraint(output_var::in, lp_constraint::in,
io.text_output_stream::in, io::di, io::uo) is det.
output_constraint(OutputVar, lte(Terms, Constant), Stream, !IO) :-
io.write_string(Stream, "le(", !IO),
output_constraint_2(OutputVar, Terms, Constant, Stream, !IO).
output_constraint(OutputVar, eq(Terms, Constant), Stream, !IO) :-
io.write_string(Stream, "eq(", !IO),
output_constraint_2(OutputVar, Terms, Constant, Stream, !IO).
output_constraint(_, gte(_,_), _, _, _) :-
unexpected($pred, "gte").
:- pred output_constraint_2(output_var::in, list(lp_term)::in, lp_constant::in,
io.text_output_stream::in, io::di, io::uo) is det.
output_constraint_2(OutputVar, Terms, Constant, Stream, !IO) :-
output_terms(OutputVar, Terms, Stream, !IO),
io.format(Stream, ", %s)", [s(rat.to_rat_string(Constant))], !IO).
:- pred output_terms(output_var::in, list(lp_term)::in,
io.text_output_stream::in, io::di, io::uo) is det.
output_terms(OutputVar, Terms, Stream, !IO) :-
io.write_char(Stream, '[', !IO),
write_out_list(output_term(OutputVar), ", ", Terms, Stream, !IO),
io.write_char(Stream, ']', !IO).
:- pred output_term(output_var::in, lp_term::in,
io.text_output_stream::in, io::di, io::uo) is det.
output_term(OutputVar, Var - Coefficient, Stream, !IO) :-
io.format(Stream, "term(%s, %s)",
[s(OutputVar(Var)), s(rat.to_rat_string(Coefficient))], !IO).
%---------------------------------------------------------------------------%
get_vars_from_constraints(Constraints) = Vars :-
list.foldl(get_vars_from_constraint, Constraints, set.init, Vars).
:- pred get_vars_from_constraint(lp_constraint::in, set(lp_var)::in,
set(lp_var)::out) is det.
get_vars_from_constraint(Constraint, !SetVar) :-
get_vars_from_terms(lp_terms(Constraint), !SetVar).
:- pred get_vars_from_terms(list(lp_term)::in, set(lp_var)::in,
set(lp_var)::out) is det.
get_vars_from_terms([], !SetVar).
get_vars_from_terms([Var - _ | Coeffs], !SetVar) :-
set.insert(Var, !SetVar),
get_vars_from_terms(Coeffs, !SetVar).
%---------------------------------------------------------------------------%
%
% Debugging predicates for writing out constraints.
%
write_constraints(Stream, VarSet, Constraints, !IO) :-
list.foldl(write_constraint(Stream, VarSet), Constraints, !IO).
:- pred write_constraint(io.text_output_stream::in, lp_varset::in,
lp_constraint::in, io::di, io::uo) is det.
write_constraint(Stream, VarSet, Constr, !IO) :-
deconstruct_constraint(Constr, Coeffs, Operator, Constant),
io.write_char(Stream, '\t', !IO),
list.foldl(write_constr_term(Stream, VarSet), Coeffs, !IO),
io.format(Stream, "%s %s\n",
[s(operator_to_string(Operator)), s(rat.to_arith_string(Constant))],
!IO).
:- pred write_constr_term(io.text_output_stream::in, lp_varset::in,
lp_term::in, io::di, io::uo) is det.
write_constr_term(Stream, VarSet, Var - Coeff, !IO) :-
VarName = varset.lookup_name(VarSet, Var),
io.format(Stream, "%s%s ",
[s(rat.to_arith_string(Coeff)), s(VarName)], !IO).
:- func operator_to_string(lp_operator) = string.
operator_to_string(lp_lt_eq) = "=<".
operator_to_string(lp_eq ) = "=".
operator_to_string(lp_gt_eq) = ">=".
:- pred write_vars(io.text_output_stream::in, varset::in, list(lp_var)::in,
io::di, io::uo) is det.
:- pragma consider_used(pred(write_vars/5)).
write_vars(Stream, VarSet, Vars, !IO) :-
VarStrs = list.map(var_to_string(VarSet), Vars),
VarsStr = string.join_list(", ", VarStrs),
io.format(Stream, "[%s]", [s(VarsStr)], !IO).
:- func var_to_string(lp_varset, lp_var) = string.
var_to_string(VarSet, Var) = Name :-
varset.lookup_name_default_prefix(VarSet, Var, "Unnamed", Name).
% Write out the matrix used during fourier elimination.
% If `Labels' is `yes', then write out the label for each vector as well.
%
:- pred write_matrix(io.text_output_stream::in, lp_varset::in,
bool::in, matrix::in, io::di, io::uo) is det.
:- pragma consider_used(pred(write_matrix/6)).
write_matrix(Stream, VarSet, WriteLabels, Matrix, !IO) :-
list.foldl(write_vector(Stream, VarSet, WriteLabels), Matrix, !IO).
:- pred write_vector(io.text_output_stream::in, lp_varset::in, bool::in,
vector::in, io::di, io::uo) is det.
write_vector(Stream, VarSet, _WriteLabels, Vector, !IO) :-
Vector = vector(_Label, Terms0, Constant),
Terms = map.to_assoc_list(Terms0),
list.foldl(write_constr_term(Stream, VarSet), Terms, !IO),
io.format(Stream, " (=<) %s\n", [s(rat.to_arith_string(Constant))], !IO).
% Write out a term - outputs the empty string if the term
% has a coefficient of zero.
%
:- pred write_term(io.text_output_stream::in, lp_varset::in, lp_term::in,
io::di, io::uo) is det.
:- pragma consider_used(pred(write_term/5)).
write_term(Stream, VarSet, Var - Coefficient, !IO) :-
( if Coefficient > zero then
Sign = "+"
else
Sign = "-"
),
Numerator = abs(numer(Coefficient)),
( if denom(Coefficient) = 1 then
MaybeDenumerator = ""
else
MaybeDenumerator =
string.format("/%s", [s(int_to_string(denom(Coefficient)))])
),
varset.lookup_name(VarSet, Var, VarName),
io.format(Stream, "%s (%d%s) %s",
[s(Sign), i(Numerator), s(MaybeDenumerator), s(VarName)], !IO).
%---------------------------------------------------------------------------%
:- end_module libs.lp_rational.
%---------------------------------------------------------------------------%
Reference in New Issue View Git Blame Copy Permalink