Files
mercury/compiler/lp_rational.m
Zoltan Somogyi 672f77c4ec Add a new compiler option. --inform-ite-instead-of-switch.
Estimated hours taken: 20
Branches: main

Add a new compiler option. --inform-ite-instead-of-switch. If this is enabled,
the compiler will generate informational messages about if-then-elses that
it thinks should be converted to switches for the sake of program reliability.

Act on the output generated by this option.

compiler/simplify.m:
	Implement the new option.

	Fix an old bug that could cause us to generate warnings about code
	that was OK in one duplicated copy but not in another (where a switch
	arm's code is duplicated due to the case being selected for more than
	one cons_id).

compiler/options.m:
	Add the new option.

	Add a way to test for the bug fix in simplify.

doc/user_guide.texi:
	Document the new option.

NEWS:
	Mention the new option.

library/*.m:
mdbcomp/*.m:
browser/*.m:
compiler/*.m:
deep_profiler/*.m:
	Convert if-then-elses to switches at most of the sites suggested by the
	new option. At the remaining sites, switching to switches would have
	nontrivial downsides. This typically happens with the switched-on type
	has many functors, and we treat one or two specially (e.g. cons/2 in
	the cons_id type).

	Perform misc cleanups in the vicinity of the if-then-else to switch
	conversions.

	In a few cases, improve the error messages generated.

compiler/accumulator.m:
compiler/hlds_goal.m:
	(Rename and) move insts for particular kinds of goal from
	accumulator.m to hlds_goal.m, to allow them to be used in other
	modules. Using these insts allowed us to eliminate some if-then-elses
	entirely.

compiler/exprn_aux.m:
	Instead of fixing some if-then-elses, delete the predicates containing
	them, since they aren't used, and (as pointed out by the new option)
	would need considerable other fixing if they were ever needed again.

compiler/lp_rational.m:
	Add prefixes to the names of the function symbols on some types,
	since without those prefixes, it was hard to figure out what type
	the switch corresponding to an old if-then-else was switching on.

tests/invalid/reserve_tag.err_exp:
	Expect a new, improved error message.
2007-11-23 07:36:01 +00:00

2409 lines
82 KiB
Mathematica

%-----------------------------------------------------------------------------%
% vim: ft=mercury ts=4 sw=4 et
%-----------------------------------------------------------------------------%
% Copyright (C) 1997-2002, 2005-2007 The University of Melbourne.
% 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 constant == rat.
:- type coefficient == rat.
:- type lp_var == var.
:- type lp_vars == list(lp_var).
:- type lp_varset == varset.
:- type lp_term == pair(lp_var, coefficient).
:- type lp_terms == list(lp_term).
% Create a term with a coefficient of 1.
% For use with ho functions.
%
:- func lp_rational.lp_term(lp_var) = lp_term.
:- type operator ---> (=<) ; (=) ; (>=).
% A primitive linear arithmetic constraint.
%
:- type constraint.
% A conjunction of primitive constraints.
%
:- type constraints == list(constraint).
% Create a constraint from the given components.
%
:- func lp_rational.constraint(lp_terms, operator, constant) = constraint.
% Create a constraint from the given components.
% Throws an exception if the resulting constraint is trivially false.
%
:- func lp_rational.non_false_constraint(lp_terms, operator, constant)
= constraint.
% Deconstruct the given constraint.
%
:- pred lp_rational.constraint(constraint::in, lp_terms::out, operator::out,
constant::out) is det.
% As above but throws an exception if the constraint is false.
%
:- pred lp_rational.non_false_constraint(constraint::in, lp_terms::out,
operator::out, constant::out) is det.
% Succeeds iff the given constraint contains a single variable and
% that variable is constrained to be a nonnegative value.
%
:- pred lp_rational.nonneg_constr(constraint::in) is semidet.
% Create a constraint that constrains the argument
% have a non-negative value.
%
:- func lp_rational.make_nonneg_constr(lp_var) = constraint.
% Create a constraint that equates two variables.
%
:- func lp_rational.make_vars_eq_constraint(lp_var, lp_var) = constraint.
% Create constraints of the form:
%
% Var = Constant or Var >= Constant
%
% These functions are useful with higher-order code.
%
:- func lp_rational.make_var_const_eq_constraint(lp_var, rat) = constraint.
:- func lp_rational.make_var_const_gte_constraint(lp_var, rat) = constraint.
% Create a constraint that is trivially false.
%
:- func lp_rational.false_constraint = constraint.
% Create a constraint that is trivially true.
%
:- func lp_rational.true_constraint = constraint.
% Succeeds if the constraint is trivially false.
%
:- pred lp_rational.is_false(constraint::in) is semidet.
% Succeeds if the constraint is trivially true.
%
:- pred lp_rational.is_true(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 lp_rational.restore_equalities(constraints::in, constraints::out)
is det.
% This checks if a constraint is entailed by all the others
% in the set. If it is then it is removed 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, constraints::in,
constraints::out) is semidet.
% Succeeds iff the given system of constraints is inconsistent.
%
:- pred lp_rational.inconsistent(lp_varset::in, constraints::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 lp_rational.simplify_constraints(constraints) = constraints.
% substitute_vars(VarsA, VarsB, Constraints0) = Constraints.
% Perform variable substitution on the given system of constriants
% based upon the mapping that is implicit between the corresponding
% elements of the variable lists `VarsA' and `VarsB'.
%
% If length(VarsA) \= length(VarsB) then an exception is thrown.
%
:- func lp_rational.substitute_vars(lp_vars, lp_vars, constraints)
= constraints.
:- func lp_rational.substitute_vars(map(lp_var, lp_var), constraints)
= constraints.
% Make the values of all the variables in the set zero.
%
:- func lp_rational.set_vars_to_zero(set(lp_var), constraints) = constraints.
%-----------------------------------------------------------------------------%
%
% 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 lp_rational.bounding_box(lp_varset, constraints) = constraints.
% Create non-negativity constraints for all of the variables
% in list of constraints except for the variables specified
% in the argument.
%
:- func lp_rational.nonneg_box(lp_vars, constraints) = constraints.
%-----------------------------------------------------------------------------%
%
% Linear solver
%
:- type objective == lp_terms.
:- 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 lp_rational.solve(constraints, direction, objective, lp_varset)
= lp_result.
%-----------------------------------------------------------------------------%
%
% Projection
%
:- type projection_result
---> pr_res_ok(constraints) % 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 callers to this procedure may need
% to do a consistency check on the result if they require
% the resulting system of constraints to be consistent.
%
:- func lp_rational.project(lp_vars, lp_varset, constraints)
= projection_result.
:- pred lp_rational.project(lp_vars::in, lp_varset::in, constraints::in,
projection_result::out) is det.
% project(Vars, Varset, maybe(MaxMatrixSize), Matrix, Result).
% Same as above but if the size of the matrix ever exceeds
% `MaxMatrixSize' we back out of the computation.
%
:- pred lp_rational.project(lp_vars::in, lp_varset::in, maybe(int)::in,
constraints::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 lp_rational.entailed(lp_varset, constraints, 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 lp_rational.entailed(lp_varset::in, constraints::in,
constraint::in) is semidet.
%-----------------------------------------------------------------------------%
%
% Stuff for intermodule optimization
%
% A function that converts an lp_var into a string.
%
:- 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 lp_rational.output_constraints(output_var::in, constraints::in,
io::di, io::uo) is det.
%-----------------------------------------------------------------------------%
%
% Debugging predicates
%
% 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 lp_rational.write_constraints(constraints::in, lp_varset::in, io::di,
io::uo) is det.
% Return the set of variables that are present in a list of constraints.
%
% XXX This shouldn't be exported but it's currently needed by the
% workaround for the problem with head variables in term_constr_fixpoint.m
%
:- func get_vars_from_constraints(constraints) = set(lp_var).
%-----------------------------------------------------------------------------%
%-----------------------------------------------------------------------------%
:- implementation.
:- import_module libs.compiler_util.
:- import_module assoc_list.
:- import_module bool.
:- import_module int.
:- import_module solutions.
:- import_module string.
:- import_module svmap.
:- import_module svset.
%-----------------------------------------------------------------------------%
%
% 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 constraint
---> lte(lp_terms, constant) % sumof(Terms) =< Constant
; eq(lp_terms, constant) % sumof(Terms) = Constant
; gte(lp_terms, constant). % sumof(Terms) >= Constant
%-----------------------------------------------------------------------------%
%
% Procedures for constructing/deconstructing constraints.
%
lp_term(Var) = Var - one.
constraint([], (=<), Const) = lte([], Const).
constraint([], (=), Const) = eq([], Const).
constraint([], (>=), Const) = lte([], -Const).
constraint(Terms0 @ [_ | _], (=<), Const0) = Constraint :-
Terms1 = sum_like_terms(Terms0),
normalize_terms_and_const(yes, Terms1, Const0, Terms, Const),
Constraint = lte(Terms, Const).
constraint(Terms0 @ [_ | _], (=), Const0) = Constraint :-
Terms1 = sum_like_terms(Terms0),
normalize_terms_and_const(no, Terms1, Const0, Terms, Const),
Constraint = eq(Terms, Const).
constraint(Terms0 @ [_ | _], (>=), Const0) = Constraint :-
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_constraint(lp_terms, operator, constant) = constraint.
unchecked_constraint(Terms, (=<), Constant) = lte(Terms, Constant).
unchecked_constraint(Terms, (=), Constant) = eq(Terms, Constant).
unchecked_constraint(Terms, (>=), Constant) = gte(Terms, Constant).
:- func sum_like_terms(lp_terms) = lp_terms.
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, coefficient)) =
map(lp_var, 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, coefficient)::in,
map(lp_var, coefficient)::in, map(lp_var, coefficient)::out) is det.
lp_terms_to_map_2(Var - Coeff0, !Map) :-
( MapCoeff = !.Map ^ elem(Var) ->
Coeff = MapCoeff + Coeff0,
( if Coeff = zero
then svmap.delete(Var, !Map)
else svmap.set(Var, Coeff, !Map)
)
;
( if Coeff0 \= zero
then svmap.set(Var, Coeff0, !Map)
else true
)
).
non_false_constraint(Terms, Op, Constant) = Constraint :-
Constraint = constraint(Terms, Op, Constant),
( if is_false(Constraint)
then unexpected(this_file,
"non_false_constraints/3: false constraint.")
else true
).
constraint(lte(Terms, Constant), Terms, (=<), Constant).
constraint(eq(Terms, Constant), Terms, (=), Constant).
constraint(gte(Terms, Constant), Terms, (>=), Constant).
non_false_constraint(Constraint, Terms, Operator, Constant) :-
( if is_false(Constraint)
then unexpected(this_file,
"non_false_constraint/4: false_constraint.")
else true
),
(
Constraint = lte(Terms, Constant),
Operator = (=<)
;
Constraint = eq(Terms, Constant),
Operator = (=)
;
Constraint = gte(_, _),
unexpected(this_file, "non_false_constraint/4: gte encountered.")
).
:- func lp_terms(constraint) = lp_terms.
lp_terms(lte(Terms, _)) = Terms.
lp_terms(eq(Terms, _)) = Terms.
lp_terms(gte(Terms, _)) = Terms.
:- func constant(constraint) = constant.
constant(lte(_, Constant)) = Constant.
constant(eq(_, Constant)) = Constant.
constant(gte(_, Constant)) = Constant.
:- func operator(constraint) = operator.
operator(lte(_, _)) = (=<).
operator(eq(_, _)) = (=).
operator(gte(_,_)) = unexpected(this_file, "operator/1: gte.").
:- func negate_operator(operator) = operator.
negate_operator((=<)) = (>=).
negate_operator((=)) = (=).
negate_operator((>=)) = (=<).
nonneg_constr(lte([_ - (-rat.one)], rat.zero)).
nonneg_constr(gte(_, _)) :-
unexpected(this_file, "nonneg_constr/1: gte.").
make_nonneg_constr(Var) = constraint([Var - (-rat.one)], (=<), rat.zero).
make_vars_eq_constraint(Var1, Var2) =
constraint([Var1 - rat.one, Var2 - (-rat.one)], (=), rat.zero).
make_var_const_eq_constraint(Var, Constant) =
constraint([Var - rat.one], (=), Constant).
make_var_const_gte_constraint(Var, Constant) =
constraint([Var - rat.one], (>=), Constant).
true_constraint = eq([], rat.zero).
false_constraint = eq([], rat.one).
is_false(gte([], Const)) :- Const > rat.zero.
is_false(lte([], Const)) :- Const < rat.zero.
is_false(eq([], Const)) :- Const \= rat.zero.
is_true(gte([], Const)) :- Const =< rat.zero.
is_true(lte([], Const)) :- Const >= rat.zero.
is_true(eq([], Const)) :- Const = rat.zero.
% Put each constraint in the list in standard form (see below).
%
:- func lp_rational.standardize_constraints(constraints) = constraints.
standardize_constraints(Constraints) =
list.map(standardize_constraint, Constraints).
% 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. (>=) is converted to (=<) by multiplying
% through by negative one. (=) constraints should have an
% initial coefficient of (positive) 1.
%
:- func lp_rational.standardize_constraint(constraint) = 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, lp_terms::in, constant::in,
lp_terms::out, 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),
( if !.Terms = [_ - Coefficient0 | _]
then
(
AbsVal = yes,
Coefficient = rat.abs(Coefficient0)
;
AbsVal = no,
Coefficient = Coefficient0
),
( if Coefficient = rat.zero
then unexpected(this_file,
"normalize_term_and_const/5: zero coefficient.")
else true
),
DivideBy = (func(Var - Coeff) = Var - (Coeff / Coefficient)),
!:Terms = list.map(DivideBy, !.Terms),
!:Const = !.Const / Coefficient
else true
).
% 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(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(constraints) = constraints.
remove_trivial([]) = [].
remove_trivial([Constraint | Constraints]) = Result :-
( is_false(Constraint) ->
Result = [ false_constraint ]
;
Result0 = remove_trivial(Constraints),
( Result0 = [C], is_false(C) ->
Result = Result0
;
% Remove the constraint if it is trivially true or the result
% of all the variables being non-negative.
( ( is_true(Constraint) ; obvious_constraint(Constraint) ) ->
Result = Result0
;
Result = [ Constraint | Result0 ]
)
)
).
:- func remove_weaker(constraints) = constraints.
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(constraint::in, constraint::in, constraints::in,
constraints::out, bool::in, bool::out) is det.
remove_weaker_2(A, B, !Acc, !Keep) :-
( is_stronger(A, B) -> true
; is_stronger(B, A) -> list.cons(B, !Acc), !:Keep = no
; list.cons(B, !Acc)
).
:- pred is_stronger(constraint::in, 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(Old, New, Constraints0) = Constraints :-
SubstMap = map.from_corresponding_lists(Old, New),
Constraints = list.map(substitute_vars_2(SubstMap), Constraints0).
substitute_vars(SubstMap, Constraints0) = Constraints :-
Constraints = list.map(substitute_vars_2(SubstMap), Constraints0).
:- func substitute_vars_2(map(lp_var, lp_var), constraint) = constraint.
substitute_vars_2(SubstMap, lte(Terms0, Const)) = Result :-
Terms = list.map(substitute_term(SubstMap), Terms0),
Result = lte(sum_like_terms(Terms), Const).
substitute_vars_2(SubstMap, eq(Terms0, Const)) = Result :-
Terms = list.map(substitute_term(SubstMap), Terms0),
Result = eq(sum_like_terms(Terms), Const).
substitute_vars_2(_, gte(_, _)) =
unexpected(this_file, "substitute_vars_2/2: gte.").
:- func substitute_term(map(lp_var, lp_var), lp_term) = lp_term.
substitute_term(SubstMap, Var - Coeff) = SubstMap ^ det_elem(Var) - Coeff.
lp_rational.set_vars_to_zero(Vars, Constraints) =
list.map(set_vars_to_zero_2(Vars), Constraints).
:- func set_vars_to_zero_2(set(lp_var), constraint) = 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), lp_terms) = lp_terms.
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 :-
Result = lp_rational.project([Var], Varset, Constrs0),
(
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 :-
( list.member(Var, VarsToIgnore) ->
true
;
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 lp_info
---> lp(
varset :: lp_varset,
slack_vars :: lp_vars, % - slack variables.
art_vars :: lp_vars % - artificial variables.
).
lp_rational.solve(Constraints, Direction, Objective, Varset) = Result :-
Info0 = lp_info_init(Varset),
solve_2(Constraints, Direction, Objective, Result, Info0, _).
% solve_2(Eqns, Dir, Obj, Res, LPInfo0, LPInfo) takes a list
% of inequalities `Eqns', a direction for optimization `Dir', an
% objective function `Obj' and an lp_info structure `LPInfo0'.
% See inline comments for details on the algorithm.
%
:- pred solve_2(constraints::in, direction::in, objective::in,
lp_result::out, lp_info::in, lp_info::out) is det.
solve_2(!.Constraints, Direction, !.Objective, Result, !LPInfo) :-
% Simplify the inequalities and convert them to standard form by
% introducing slack/artificial variables.
Obj = !.Objective,
lp_standardize_constraints(!Constraints, !LPInfo),
% 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 = !.LPInfo ^ 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(lp_terms, lp_terms, 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(lp_terms, lp_terms, lp_vars, 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 Result = lp_res_inconsistent
else
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, _)
)
).
%-----------------------------------------------------------------------------%
:- pred lp_standardize_constraints(constraints::in, constraints::out,
lp_info::in, lp_info::out) is det.
lp_standardize_constraints(!Constraints, !LPInfo) :-
list.map_foldl(lp_standardize_constraint, !Constraints, !LPInfo).
% 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(constraint::in, constraint::out,
lp_info::in, lp_info::out) is det.
lp_standardize_constraint(Constr0 @ lte(Coeffs, Const), Constr, !LPInfo) :-
( Const < zero ->
Constr1 = negate_constraint(Constr0),
lp_standardize_constraint(Constr1, Constr, !LPInfo)
;
new_slack_var(Var, !LPInfo),
Constr = lte([Var - one | Coeffs], Const)
).
lp_standardize_constraint(Eqn0 @ eq(Coeffs, Const), Eqn, !LPInfo) :-
( Const < zero ->
Eqn1 = negate_constraint(Eqn0),
lp_standardize_constraint(Eqn1, Eqn, !LPInfo)
;
new_art_var(Var, !LPInfo),
Eqn = lte([Var - one | Coeffs], Const)
).
lp_standardize_constraint(Eqn0 @ gte(Coeffs, Const), Eqn, !LPInfo) :-
( Const < zero ->
Eqn1 = negate_constraint(Eqn0),
lp_standardize_constraint(Eqn1, Eqn, !LPInfo)
;
new_slack_var(SVar, !LPInfo),
new_art_var(AVar, !LPInfo),
Eqn = gte([AVar - one, SVar - (-one) | Coeffs], Const)
).
:- func negate_constraint(constraint) = 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(lp_terms) = lp_terms.
negate_lp_terms(Terms) = assoc_list.map_values((func(_, X) = (-X)), Terms).
:- func add_var(map(lp_var, rat), lp_var, rat) = map(lp_var, rat).
add_var(Map0, Var, Coeff) = Map :-
Acc1 = ( if Acc0 = Map0 ^ elem(Var) then Acc0 else zero ),
Acc = Acc1 + Coeff,
Map = Map0 ^ elem(Var) := Acc.
%-----------------------------------------------------------------------------%
:- func collect_vars(constraints, 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(lp_vars, int) = var_num_map.
number_vars(Vars, N) = VarNum :-
number_vars_2(Vars, N, map.init, VarNum).
:- pred number_vars_2(lp_vars::in, int::in,
var_num_map::in, var_num_map::out) is det.
number_vars_2([], _, !VarNums).
number_vars_2([Var | Vars], N, !VarNums) :-
svmap.det_insert(Var, N, !VarNums),
number_vars_2(Vars, N + 1, !VarNums).
:- pred insert_constraints(constraints::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(lp_terms::in, int::in, var_num_map::in,
tableau::in, tableau::out) is det.
insert_terms([], _, _, !Tableau).
insert_terms([Var - Const | Coeffs], Row, VarNums, !Tableau) :-
Col = VarNums ^ det_elem(Var),
set_cell(Row, Col, Const, !Tableau),
insert_terms(Coeffs, Row, VarNums, !Tableau).
%-----------------------------------------------------------------------------%
:- pred optimize(lp_vars::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(lp_vars, 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 = Map0 ^ elem(Var) := Val.
:- 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),
!:Min = ( if MinVal < zero then yes(Col - MinVal) else 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(this_file,
"simplex/3: 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 true % CellVal = 0 => multiple optimal sol'ns.
else
( if CellVal = zero
then unexpected(this_file,
"simplex/3: 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),
simplex(Result, !Tableau)
)
).
%-----------------------------------------------------------------------------%
:- pred ensure_zero_obj_coeffs(lp_vars::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),
( Val = zero ->
ensure_zero_obj_coeffs(Vars, !Tableau)
;
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(this_file,
"ensure_zero_obj_coeffs/3: zero divisor.")
else true
),
Fac = -Val / Fac0,
row_op(Fac, Row, 0, !Tableau),
ensure_zero_obj_coeffs(Vars, !Tableau)
;
Ones = [],
unexpected(this_file, "ensure_zero_obj_coeffs/3: " ++
"problem with artificial variable.")
)
).
:- pred fix_basis_and_rem_cols(lp_vars::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 = ( Val = zero -> Ones0 ; [Val - R | Ones0] )
),
solutions.aggregate(all_rows(!.Tableau), BasisAgg, [], Res),
(
Res = [one - Row]
->
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)
)
;
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(this_file, "pivot/4 - 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(this_file, "pivot/4 - 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(this_file,
"get_cell/3: attempt to address shunned row/col.")
else true
),
( if Cell0 = Tableau ^ cells ^ elem(Row - Col)
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(this_file,
"set_cell/5: 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 = !.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) = lp_vars.
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 lp_info_init(lp_varset) = lp_info.
lp_info_init(Varset) = lp(Varset, [], []).
:- pred new_slack_var(lp_var::out, lp_info::in, lp_info::out) is det.
new_slack_var(Var, !LPInfo) :-
varset.new_var(!.LPInfo ^ varset, Var, Varset),
!:LPInfo = !.LPInfo ^ varset := Varset,
Vars = !.LPInfo ^ slack_vars,
!:LPInfo = !.LPInfo ^ slack_vars := [Var | Vars].
:- pred new_art_var(lp_var::out, lp_info::in, lp_info::out) is det.
new_art_var(Var, !LPInfo) :-
varset.new_var(!.LPInfo ^ varset, Var, Varset),
!:LPInfo = !.LPInfo ^ varset := Varset,
Vars = !.LPInfo ^ art_vars,
!:LPInfo = !.LPInfo ^ 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(
label :: set(int),
% The vector's label is for redundancy checking
% during Fourier elimination - see below.
terms :: map(lp_var, coefficient),
% A map from each variable in the vector to its
% coefficient
const :: constant
).
:- type matrix == list(vector).
project(Vars, Varset, Constraints) = Result :-
project(Vars, Varset, no, Constraints, Result).
project(Vars, Varset, Constraints, Result) :-
project(Vars, Varset, no, Constraints, Result).
% For the first branch of this switch 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, pr_res_ok(Constraints)).
project(!.Vars @ [_|_], Varset, MaybeThreshold, 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(this_file, "project/5: 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(constraints) = matrix.
constraints_to_matrix(Constraints) = Matrix :-
list.foldl2(fm_standardize, Constraints, 0, _, [], Matrix).
:- pred fm_standardize(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) :-
set.singleton_set(Label, Labels).
:- func matrix_to_constraints(matrix) = constraints.
matrix_to_constraints(Matrix) = list.map(vector_to_constraint, Matrix).
:- func vector_to_constraint(vector) = 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
% (ie. 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(lp_vars::in, lp_vars::out, constraints::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),
(
eliminate_equations_2(!Vars, Equalities0, Equalities,
Inequalities0, Inequalities)
->
Result = pr_res_ok(Equalities ++ Inequalities)
;
Result = pr_res_inconsistent
).
:- pred eliminate_equations_2(lp_vars::in, lp_vars::out,
constraints::in, constraints::out, constraints::in,
constraints::out) is semidet.
eliminate_equations_2([], [], !Equations, !Inequations).
eliminate_equations_2([Var | !.Vars], !:Vars, !Equations, !Inequations) :-
eliminate_equations_2(!Vars, !Equations, !Inequations),
( find_target_equality(Var, Target, !Equations) ->
substitute_variable(Target, Var, !Equations, !Inequations,
SuccessFlag),
(
SuccessFlag = no,
list.cons(Var, !Vars),
list.cons(Target, !Equations)
;
SuccessFlag = yes
)
;
list.cons(Var, !Vars)
).
% Find an equation that constrains a variable we are trying
% to eliminate.
%
:- pred find_target_equality(lp_var::in, constraint::out, constraints::in,
constraints::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, constraints) =
maybe(pair(constraint, constraints)).
find_target_equality(Var, Eqns) = find_target_equality_2(Var, Eqns, []).
:- func find_target_equality_2(lp_var, constraints, constraints) =
maybe(pair(constraint, constraints)).
find_target_equality_2(_, [], _) = no.
find_target_equality_2(Var, [Eqn | Eqns], Acc) = MaybeTargetEqn :-
( if operator(Eqn) \= (=)
then unexpected(this_file,
"find_target_equality_2/3: inequality encountered.")
else true
),
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(constraint::in, lp_var::in,
constraints::in, constraints::out, constraints::in, constraints::out,
bool::out) is semidet.
substitute_variable(Target0, Var, !Equations, !Inequations, Flag) :-
normalize_constraint(Var, Target0, Target),
constraint(Target, TargetCoeffs, Op, TargetConst),
(
Op = (=)
;
( Op = (=<)
; Op = (>=)
),
unexpected(this_file, "substitute_variable/7: 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, lp_terms::in, constant::in,
lp_terms::out, 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),
FixedCoeffs = ( Var = Var1 -> FCoeffs0 ; [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, lp_terms::in,
constant::in, constraints::in, constraints::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, lp_terms::in,
constant::in, constraint::in, constraint::out, bool::out) is det.
substitute_into_constraint(Var, SubCoeffs, SubConst, !Constraint, Flag) :-
normalize_constraint(Var, !Constraint),
constraint(!.Constraint, TargetCoeffs, Op, TargetConst),
( list.member(Var - one, TargetCoeffs) ->
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 = constraint(FinalCoeffs, Op, FinalConst),
Flag = yes
;
Flag = no
).
%-----------------------------------------------------------------------------%
%
% Fourier elimination.
%
% Will return `no' if it aborts otherwise `yes(Matrix)', where
% `Matrix' is the result of the projection.
%
:- pred fourier_elimination(lp_vars::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.
(
PosMatrix = [_ | _],
NegMatrix = [_ | _]
->
!:Step = !.Step + 1,
(
list.foldl2(eliminate_var(!.Step, MaybeThreshold, NegMatrix),
PosMatrix, ZeroMatrix, ResultMatrix, SizeZeroMatrix, _)
->
NewMatrix = yes(ResultMatrix)
;
NewMatrix = no
)
;
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) :-
( Coefficient = Vector0 ^ terms ^ elem(Var) ->
Vector0 = vector(Label, Terms0, Const0),
normalize_vector(Var, Terms0, Const0, Terms, Const),
Vector1 = vector(Label, Terms, Const),
( if Coefficient > zero
then list.cons(Vector1, !Pos)
else list.cons(Vector1, !Neg)
)
;
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 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
->
add_vectors(TermsPos, ConstPos, TermsNeg, ConstNeg, Coeffs,
Const),
New = vector(LabelNew, Coeffs, Const),
(
(
% 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)
)
->
% If the new constraint is `true' or is
% quasi-syntactic redundant with something
% already there.
true
;
% 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),
(
list.member(Vec, !.Zeros),
label_subsumed(New, Vec)
->
% Do not add the new constraint because it is label
% subsumed by something already in the matrix.
%
true
;
filter_and_count(
(pred(Vec2::in) is semidet :-
not label_subsumed(Vec2, New)
),
!.Zeros, [], !:Zeros, 0, !:Num),
list.cons(New, !Zeros),
!:Num = !.Num + 1
)
)
;
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(lp_vars::in, matrix::in, lp_var::out,
lp_vars::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) = lp_vars.
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(this_file, "find_max/2: empty list passed as arg.").
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(lp_vars, 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) :-
( !.Map ^ elem(Var) = coeff_info(Pos0, Neg0) ->
( Coeff > zero ->
Pos = Pos0 + 1, Neg = Neg0
; Coeff < zero ->
Pos = Pos0, Neg = Neg0 + 1
;
unexpected(this_file,
"count_coeff/3: zero coefficient encountered.")
),
svmap.det_update(Var, coeff_info(Pos, Neg), !Map)
;
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(lp_vars) = 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, Const0, Terms, Const).
% Multiply the given vector by a scalar appropraite 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, coefficient)::in,
constant::in, map(lp_var, coefficient)::out, constant::out) is det.
normalize_vector(Var, !.Terms, !.Constant, !:Terms, !:Constant) :-
( Coefficient = !.Terms ^ elem(Var) ->
( if Coefficient = zero
then unexpected(this_file,
"normalize_vector/5: zero coefficient in vector.")
else true
),
DivVal = rat.abs(Coefficient),
!:Terms = map.map_values((func(_, C) = C / DivVal), !.Terms),
!:Constant = !.Constant / DivVal
;
% In this case the 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 scaler 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, constraint::in, constraint::out)
is det.
normalize_constraint(Var, Constraint0, Constraint) :-
lp_rational.constraint(Constraint0, Terms0, Op0, Constant0),
( assoc_list.search(Terms0, Var, Coefficient) ->
( if Coefficient = zero
then unexpected(this_file,
"normalize_constraint/3: zero coefficient constraint.")
else true
),
Terms = list.map((func(V - C) = V - (C / Coefficient)), Terms0),
Constant = Constant0 / Coefficient,
( if Coefficient < zero
then Op = negate_operator(Op0)
else Op = Op0
)
;
% In this case the 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 = lp_rational.unchecked_constraint(Terms, Op, Constant).
:- pred add_vectors(map(lp_var, coefficient)::in, constant::in,
map(lp_var, coefficient)::in, constant::in,
map(lp_var, coefficient)::out, 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 :-
NumA = TermsA ^ det_elem(Var),
( if Coeffs0 ^ elem(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).
%-----------------------------------------------------------------------------%
%
% 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, constraints::in,
constraints::in, constraints::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) :-
( obvious_constraint(E) ->
true
;
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).
%-----------------------------------------------------------------------------%
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(constraint::in, constraints::in, constraints::in,
constraint::out, constraints::out) is semidet.
check_for_equalities(Eqn0, [Eqn | Eqns], SoFar, NewEqn, NewEqnSet) :-
(
opposing_inequalities(Eqn0 @ lte(Coeffs, Constant), Eqn)
->
NewEqn = standardize_constraint(eq(Coeffs, Constant)),
NewEqnSet = SoFar ++ Eqns
;
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(constraint::in, constraint::in) is semidet.
opposing_inequalities(lte(TermsA, Const), lte(TermsB, -Const)) :-
TermsB = list.map((func(V - X) = V - (-X)), TermsA).
%-----------------------------------------------------------------------------%
%-----------------------------------------------------------------------------%
%
% 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(this_file, "entailed/3: inconsistent constraint set.")
;
Result = not_entailed,
fail
).
%-----------------------------------------------------------------------------%
get_vars_from_constraints(Constraints) = Vars :-
list.foldl(get_vars_from_constraint, Constraints, set.init, Vars).
:- pred get_vars_from_constraint(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(lp_terms::in, set(lp_var)::in, set(lp_var)::out)
is det.
get_vars_from_terms([], !SetVar).
get_vars_from_terms([Var - _ | Coeffs], !SetVar) :-
svset.insert(Var, !SetVar),
get_vars_from_terms(Coeffs, !SetVar).
%-----------------------------------------------------------------------------%
%
% Printing constraints.
%
% Write out a term - outputs the empty string if the term
% has a coefficient of zero.
%
:- pred write_term(lp_varset::in, lp_term::in, io::di, io::uo) is det.
write_term(Varset, Var - Coefficient, !IO) :-
( if Coefficient > zero
then io.write_char('+', !IO)
else io.write_char('-', !IO)
),
io.write_string(" (", !IO),
Num = abs(numer(Coefficient)),
io.write_string(int_to_string(Num), !IO),
( if denom(Coefficient) \= 1
then io.format("/%s", [s(int_to_string(denom(Coefficient)))], !IO)
else true
),
io.write_char(')', !IO),
io.write_string(varset.lookup_name(Varset, Var), !IO).
%-----------------------------------------------------------------------------%
%
% Debugging predicates for writing out constraints.
%
write_constraints(Constraints, Varset, !IO) :-
list.foldl(write_constraint(Varset), Constraints, !IO).
:- pred write_constraint(lp_varset::in, constraint::in, io::di, io::uo) is det.
write_constraint(Varset, Constr, !IO) :-
constraint(Constr, Coeffs, Operator, Constant),
io.write_char('\t', !IO),
list.foldl(write_constr_term(Varset), Coeffs, !IO),
io.format("%s %s\n", [s(operator_to_string(Operator)),
s(rat.to_string(Constant))], !IO).
:- pred write_constr_term(lp_varset::in, lp_term::in, io::di, io::uo) is det.
write_constr_term(Varset, Var - Coeff, !IO) :-
VarName = varset.lookup_name(Varset, Var),
io.format("%s%s ", [s(rat.to_string(Coeff)), s(VarName)], !IO).
:- func operator_to_string(operator) = string.
operator_to_string((=<)) = "=<".
operator_to_string((=) ) = "=".
operator_to_string((>=)) = ">=".
:- pred write_vars(varset::in, lp_vars::in, io::di, io::uo) is det.
write_vars(Varset, Vars, !IO) :-
io.write_string("[ ", !IO),
write_vars_2(Varset, Vars, !IO),
io.write_string(" ]", !IO).
:- pred write_vars_2(lp_varset::in, lp_vars::in, io::di, io::uo) is det.
write_vars_2(_, [], !IO).
write_vars_2(Varset, [V | Vs], !IO) :-
io.write_string(var_to_string(Varset, V), !IO),
(
Vs = []
;
Vs = [_ | _],
io.write_string(", ", !IO)
),
write_vars_2(Varset, Vs, !IO).
:- func var_to_string(lp_varset, lp_var) = string.
var_to_string(Varset, Var) = varset.lookup_name(Varset, Var, "Unnamed").
% 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(lp_varset::in, bool::in, matrix::in, io::di, io::uo)
is det.
write_matrix(Varset, Labels, Matrix, !IO) :-
io.write_list(Matrix, "\n", write_vector(Varset, Labels), !IO).
:- pred write_vector(lp_varset::in, bool::in, vector::in, io::di,
io::uo) is det.
write_vector(Varset, _WriteLabels, vector(_Label, Terms0, Constant), !IO) :-
Terms = map.to_assoc_list(Terms0),
list.foldl(write_constr_term(Varset), Terms, !IO),
io.write_string(" (=<) ", !IO),
io.write_string(rat.to_string(Constant), !IO).
%-----------------------------------------------------------------------------%
%
% Intermodule optimization stuff.
%
% The following predicates write out constraints in a form that is useful
% for (transitive) intermodule optimization.
output_constraints(OutputVar, Constraints, !IO) :-
io.write_char('[', !IO),
io.write_list(Constraints, ", ", output_constraint(OutputVar), !IO),
io.write_char(']', !IO).
:- pred output_constraint(output_var::in, constraint::in,
io::di, io::uo) is det.
output_constraint(OutputVar, lte(Terms, Constant), !IO) :-
io.write_string("le(", !IO),
output_constraint_2(OutputVar, Terms, Constant, !IO).
output_constraint(OutputVar, eq(Terms, Constant), !IO) :-
io.write_string("eq(", !IO),
output_constraint_2(OutputVar, Terms, Constant, !IO).
output_constraint(_, gte(_,_), _, _) :-
unexpected(this_file, "output_constraint/3: gte encountered.").
:- pred output_constraint_2(output_var::in, lp_terms::in,
constant::in, io::di, io::uo) is det.
output_constraint_2(OutputVar, Terms, Constant, !IO) :-
output_terms(OutputVar, Terms, !IO),
io.write_string(", ", !IO),
rat.write_rat(Constant, !IO),
io.write_char(')', !IO).
:- pred output_terms(output_var::in, lp_terms::in, io::di, io::uo)
is det.
output_terms(OutputVar, Terms, !IO) :-
io.write_char('[', !IO),
io.write_list(Terms, ", ", output_term(OutputVar), !IO),
io.write_char(']', !IO).
:- pred output_term(output_var::in, lp_term::in, io::di, io::uo)
is det.
output_term(OutputVar, Var - Coefficient, !IO) :-
io.format("term(%s, ", [s(OutputVar(Var))], !IO),
rat.write_rat(Coefficient, !IO),
io.write_char(')', !IO).
%-----------------------------------------------------------------------------%
:- func this_file = string.
this_file = "lp_rational.m".
%-----------------------------------------------------------------------------%
:- end_module libs.lp_rational.
%-----------------------------------------------------------------------------%