[m-rev.] for review: new termination analyser (part 2 of 6)

Julien Fischer juliensf at cs.mu.OZ.AU
Thu Mar 24 15:44:40 AEDT 2005


===================================================================
RCS file: compiler/lp_rational.m
diff -N compiler/lp_rational.m
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ compiler/lp_rational.m	23 Mar 2005 22:47:32 -0000
@@ -0,0 +1,2366 @@
+%-----------------------------------------------------------------------------%
+% vim: ft=mercury ts=4 sw=4 et
+%-----------------------------------------------------------------------------%
+% Copyright (C) 1997-2002 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 above.)
+%
+%-----------------------------------------------------------------------------%
+
+:- module libs.lp_rational.
+
+:- interface.
+
+:- import_module libs.rat.
+
+:- import_module io.
+:- import_module list.
+:- import_module map.
+:- import_module set.
+:- import_module std_util.
+:- 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.
+
+:- func lp_rational.simplify_constraints(constraints) = constraints.
+
+    % substitute_vars(VarsA, VarsB, Constraints0) = Constraints.
+    % 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
+    --->    unbounded
+    ;       inconsistent
+    ;       satisfiable(rat, map(lp_var, rat)).
+                % 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
+    --->    ok(constraints) % projection succeeded.
+    ;       inconsistent    % matrix was inconsistent.
+    ;       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 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(Vars , Vars):
+    % 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'.
+    % Then tests if `P' satisfies `C'.
+    % This assumes that all the variables are non-negative.
+    % Throws an exception if `Cs' is inconsistent.
+    %
+:- pred lp_rational.entailed(lp_varset::in, constraints::in,
+    constraint::in) is semidet.
+
+:- func lp_rational.entailed(lp_varset, constraints, constraint) =
+    entailment_result.
+
+%------------------------------------------------------------------------------%
+%
+% Stuff for intermodule optimization.
+%
+
+    % A function that converts an lp_var into a string.
+    %
+:- type output_var == (func(lp_var) = string).
+:- inst output_var == (func(in) = out is det).
+
+    % 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(output_var),
+    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 assoc_list.
+:- import_module bool.
+:- import_module exception.
+:- import_module int.
+:- 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  throw("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  throw("non_false_constraint/4: false_constraint.")
+      else  true
+    ),
+    (
+        Constraint = lte(Terms, Constant),
+        Operator   = (=<)
+    ;
+        Constraint = eq(Terms, Constant),
+        Operator   = (=)
+    ;
+        Constraint = gte(_, _),
+        throw("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(_,_))  = throw("operator/1: gte.").
+
+:- func negate_operator(operator) = operator.
+
+negate_operator((=<)) = (>=).
+negate_operator((=))  = (=).
+negate_operator((>=)) = (=<).
+
+nonneg_constr(lte([_ - (-rat.one)], rat.zero)).
+nonneg_constr(gte(_, _)) :-
+    throw("nonneg_constr/1: get.").
+
+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 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  throw("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],
+        Result = lp_rational.solve(Constraints, max, DummyObjective, Vars),
+        Result = 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(_, _)) = throw("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)),
+    BoundingBox = list.foldl((func(Var, Constrs0) = Constrs :-
+            Result = lp_rational.project([Var], Varset, Constrs0),
+            (
+                Result = inconsistent,
+                Constrs = [false_constraint]
+            ;
+                % If we needed to abort this computation
+                % we will just approximate the whole lot
+                % by `true'.
+                Result = aborted,
+                Constrs = []
+            ;
+                Result = ok(Constrs)
+            )
+        ), 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 = unbounded,
+            Result = Result0
+        ;
+            Result0 = inconsistent,
+            Result = Result0
+        ;
+            Result0 = satisfiable(NOptVal, OptCoffs),
+            OptVal = -NOptVal,
+            Result = 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 = unbounded,
+        Result = unbounded
+    ;
+        Result0 = inconsistent,
+        Result = inconsistent
+    ;
+        Result0 = satisfiable(Val, _ArtRes),
+        ( if    Val \= zero
+          then  Result = 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)
+    ),
+    std_util.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  = unbounded
+    ;
+        Result0 = yes,
+        ObjVal  = !.Tableau ^ elem(0, !.Tableau ^ cols),
+        ObjMap  = extract_objective(ObjVars, !.Tableau),
+        Result  = 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)
+    ),
+    std_util.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 )
+        )
+    ),
+    std_util.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  throw("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  throw("simplex/3: zero divisor.")
+                          else  true
+                        ),
+                        MaxVal1 = MVal / CellVal,
+                        ( if    MaxVal1 =< MaxVal0
+                          then  !:Max = yes(Row - MaxVal1)
+                          else  true
+                        )
+                )
+            )
+        ),
+        std_util.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
+        ),
+        std_util.solutions(FindOne, Ones),
+        (
+            Ones = [Row - Fac0 | _],
+            ( if    Fac0 = zero
+              then  throw("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 = [],
+            throw("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] )
+    ),
+    std_util.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
+        ),
+        std_util.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  throw("pivot/4 - ScaleCell: zero divisor.")
+          else  true
+        ),
+        T = T0 ^ elem(J, K) := Ajk - Apk * Ajq / Apq
+    ),
+    std_util.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
+    ),
+    std_util.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  throw("pivot/4 - ScaleRow: zero divisor.")
+          else  true
+        ),
+        T = T0 ^ elem(P, K) := Apk / Apq
+    ),
+    std_util.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
+    ),
+    std_util.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  throw("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  throw("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
+        ),
+        std_util.solutions(NonZeroGoal, Solns),
+        Solns = [_ - one]
+    ),
+    std_util.solutions(BasisCol, Cols),
+    BasisVars = (pred(V::out) is nondet :-
+        list.member(Col, Cols),
+        map.member(Tableau ^ var_nums, V, Col)
+    ),
+    std_util.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, ok(Constraints)).
+project(!.Vars @ [_|_], Varset, MaybeThreshold, Constraints0, Result) :-
+    eliminate_equations(!Vars, Constraints0, EqlResult),
+    (
+        EqlResult = inconsistent,
+        Result = inconsistent
+    ;
+        % Elimination of equations should not cause an abort
+        % since we are always making the matrix smaller.
+        %
+        EqlResult = aborted,
+        throw("project/5: abort from eliminate_equations.")
+    ;
+        EqlResult = 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 = ok(Constraints)
+            ;
+                FourierResult = no,
+                Result = 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.
+            Result = 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 then substitute the value of that
+% variable into the other constraints.
+
+:- 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 = ok(Equalities ++ Inequalities)
+    ;
+        Result = 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)
+        ;
+            true
+        )
+    ;
+        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  throw("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),
+    ( if    Op \= (=)
+      then  throw("substitute_variable/7: inequality encountered.")
+      else  true
+    ),
+    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 a one or more substitutions was 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 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)
+    ),
+    ( if    NewMatrix = yes(Matrix)
+      then  fourier_elimination(OtherVars, Varset, MaybeThreshold, !.Step,
+            Matrix, Result)
+      else  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), matrix, matrix, matrix, int, int) is det.
+:- mode filter_and_count(pred(in) is semidet, in, in, out, in, 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([]) = throw("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
+        ;
+            throw("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).
+    % Converts a vector, broken up into `Terms0' and `Const0', into a
+    % parallel one where the coefficient of `Var' is +/- one.
+    % Throws an exception if the map contains a zero coefficient for `Var'.
+    %
+:- 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  throw(
+            "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  throw("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)
+        )
+    ),
+    std_util.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(C, Cs, Vars).
+    % 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 (or minimizes if `C' contains `>=' )
+    % a function parallel to `C'.  Then tests whether `P' satisfies `C'.
+    % This assumes that all the variables are non-negative.
+    %
+entailed(Varset, Constraints, lte(Objective, Constant)) = Result :-
+    SolverResult = lp_rational.solve(Constraints, max, Objective, Varset),
+    (
+        SolverResult = satisfiable(MaxVal, _),
+        Result = ( if MaxVal =< Constant then entailed else not_entailed )
+    ;
+        SolverResult = unbounded,
+        Result = not_entailed
+    ;
+        SolverResult = 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 = Result
+    ).
+entailed(Varset, Constraints, gte(Objective, Constant)) = Result :-
+    SolverResult = lp_rational.solve(Constraints, min, Objective, Varset),
+    (
+        SolverResult = satisfiable(MinVal, _),
+        Result = ( if MinVal >= Constant then entailed else not_entailed )
+    ;
+        SolverResult = unbounded,
+        Result = not_entailed
+    ;
+        SolverResult = inconsistent,
+        Result = inconsistent
+    ).
+
+entailed(Varset, Constraints, Constraint) :-
+    Result = entailed(Varset, Constraints, Constraint),
+    (
+        Result = entailed
+    ;
+        Result = inconsistent,
+        throw("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),
+    ( if Vs = [] then true else 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(output_var), 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(_,_), _, _) :-
+    throw("output_constraint/3: gte encountered.").
+
+:- pred output_constraint_2(output_var::in(output_var), 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(output_var), 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(output_var), 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.
+%------------------------------------------------------------------------------%

--------------------------------------------------------------------------
mercury-reviews mailing list
post:  mercury-reviews at cs.mu.oz.au
administrative address: owner-mercury-reviews at cs.mu.oz.au
unsubscribe: Address: mercury-reviews-request at cs.mu.oz.au Message: unsubscribe
subscribe:   Address: mercury-reviews-request at cs.mu.oz.au Message: subscribe
--------------------------------------------------------------------------



More information about the reviews mailing list