mercury/compiler/polyhedron.m

%------------------------------------------------------------------------------%
% Copyright (C) 2003, 2005 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: polyhedron.m
% main author: juliensf

% Provides closed convex polyhedra over Q^n.
% These are useful as an abstract domain for describing numerical 
% relational information.  

% The set of closed convex polyhedra is partially ordered by subset
% inclusion.  It forms a lattice with intersection as the binary meet
% operation and convex hull as the binary join operation.

% This module includes procedures for:
% 	- computing the intersection of two convex polyhedra
% 	- computing the convex hull of two convex polyhedra
% 	- approximating the convex union of two convex polyhedra by means
% 	  other than the convex hull when that becomes too computationally
% 	  expensive.
% 	- converting a convex polyhedron to and from an equivalent system
% 	  of linear constraints.

% It also includes an implementation of widening for convex polyhedra.

% NOTE: many of the operations in this module require you to pass in 
% the varset that the variables in the constraints that define the polyhedron
% were allocated from.  This because the code for computing the convex hull
% and the linear solver in lp_rational.m need to allocate fresh temporary
% variables.  
%
% XXX We could avoid this with some extra work.

% TODO:
%	* See if using the double description method is any faster.

%------------------------------------------------------------------------------%

:- module libs.polyhedron.

:- interface.

:- import_module libs.lp_rational.

:- import_module io.
:- import_module list.
:- import_module map.
:- import_module set.
:- import_module std_util.

%------------------------------------------------------------------------------%

:- type polyhedron. 

:- type polyhedra == list(polyhedron).

%------------------------------------------------------------------------------%

	% The `empty' polyhedron.  Equivalent to the constraint `false'.
	%
:- func polyhedron.empty = polyhedron.

	% The `universe' polyhedron.  Equivalent to the constraint `true'.
	%   
:- func polyhedron.universe = polyhedron.

	% Constructs a convex polyhedron from a system of linear constraints.
	%
:- func polyhedron.from_constraints(constraints) = polyhedron.

	% Returns a system of constraints whose solution space defines
	% the given polyhedron.
	%
:- func polyhedron.constraints(polyhedron) = constraints.
	
	% As above but throws an exception if the given polyhedron is empty.
	%
:- func polyhedron.non_false_constraints(polyhedron) = constraints.

	% Succeeds iff the given polyhedron is the empty polyhedron.
	%
	% NOTE: this only succeeds if the polyhedron is actually *known* 
	% to be empty.  It might fail even when the constraint set is
	% is inconsistent.  You currently need to call polyhedron.optimize
	% to force this to always work. 
	%	
:- pred polyhedron.is_empty(polyhedron::in) is semidet.

	% Succeeds iff the given polyhedron is the `universe' polyhedron,
	% that is the one whose constraint representation corresponds to `true'.
	% (ie. it is unbounded in all dimensions).
	%
:- pred polyhedron.is_universe(polyhedron::in) is semidet.

	% Optimizes the representation of a polyhedron.
	% At the moment this performs a consistency check and then
	% replaces the polyhedron by empty if necessary.
	%
:- pred polyhedron.optimize(lp_varset::in, polyhedron::in, polyhedron::out)
	is det.

	% polyhedron.intersection(A, B, C).
	% The polyhedron `C' is the intersection of the the
	% polyhedra `A' and `B'. 
	%
:- func polyhedron.intersection(polyhedron, polyhedron) = polyhedron.
:- pred polyhedron.intersection(polyhedron::in, polyhedron::in,
	polyhedron::out) is det.

	% Returns a polyhedron that is a closed convex approximation of 
	% union of the two polyhedra.
	%
:- func polyhedron.convex_union(lp_varset, polyhedron, polyhedron) = 
	polyhedron.
:- pred polyhedron.convex_union(lp_varset::in, polyhedron::in, polyhedron::in,
	polyhedron::out) is det.

	% As above but takes an extra argument that weakens the approximation
	% even further if the size of the internal matrices exceeds the 
	% supplied threshold
	% 
:- func polyhedron.convex_union(lp_varset, maybe(int), polyhedron, 
	polyhedron) = polyhedron.
:- pred polyhedron.convex_union(lp_varset::in, maybe(int)::in, polyhedron::in,
	polyhedron::in, polyhedron::out) is det.

	% Approximate a (convex) polyhedron by a rectangular region
	% whose sides are parallel to the axes.
	%
:- func polyhedron.bounding_box(polyhedron, lp_varset) = polyhedron.

	% polyhedron.widen(A, B, Varset) = C.
	% Remove faces from the polyhedron `A' to form the polyhedron `C' 
	% according to the rules that the smallest number of faces
	% should be removed and that `C' must be a superset of `B'.
	% This operation is not commutative.
	%
:- func polyhedron.widen(polyhedron, polyhedron, lp_varset) = polyhedron. 

	% project_all(Varset, Variables, Polyhedra) returns a list
	% of polyhedra in which the variables listed have been eliminated
	% from each polyhedron. 
	%
:- func polyhedron.project_all(lp_varset, lp_vars, polyhedra) = polyhedra.

:- func polyhedron.project(lp_vars, lp_varset, polyhedron) = polyhedron.
:- pred polyhedron.project(lp_vars::in, lp_varset::in, polyhedron::in,
	polyhedron::out) is det.

	% XXX It might be nicer to think of this as relabelling the axes.
	% Conceptually it alters the names of the variables in the polyhedron
	% without explicitly converting it back into constraint form - this
	% is easy to do (at the moment) as the polyhedra are represented
	% as constraints anyway.
	%
:- func polyhedron.substitute_vars(lp_vars, lp_vars, polyhedron) = polyhedron.
:- func polyhedron.substitute_vars(map(lp_var, lp_var), polyhedron) = 
	polyhedron.

	% polyhedron.zero_vars(Set, Polyhedron0) = Polyhedron <=>
	%
	% 	Constraints0 = polyhedron.constraints(Polyhedron0),
	% 	Constraints  = lp_rational.set_vars_to_zero(
	% 		Set, Constraints0)
	% 	Polyhedron   = polyhedron.from_constraints(Constraints)
	%
	% This is a little more efficient than the above because
	% we don't end up traversing the list of constraints as much.
	% 	
:- func polyhedron.zero_vars(set(lp_var), polyhedron) = polyhedron.

	% Print out the polyhedron using the names of the variables in the 
	% varset.
	%
:- pred polyhedron.write_polyhedron(polyhedron::in, lp_varset::in, io::di,
	io::uo) is det.

%------------------------------------------------------------------------------%
%------------------------------------------------------------------------------%

:- implementation.

:- import_module libs.rat.

:- import_module bool.
:- import_module exception.
:- import_module int.
:- import_module string.
:- import_module svmap.
:- import_module svvarset.
:- import_module varset.

%------------------------------------------------------------------------------%

	% XXX The constructor eqns/1 should really be called something
	% more meaningful.
	%
:- type polyhedron 
	--->	eqns(constraints)
	;	empty_poly.

%------------------------------------------------------------------------------%
%
% Creation of polyhedra.
%

polyhedron.empty = empty_poly.

polyhedron.universe = eqns([]).
	
	% This does the following:
	% 	- checks if the constraint is false.
	% 	- simplifies the representation of the constraint.
	%	 - calls intersection/3 (which does further simplifications).
	%
polyhedron.from_constraints([]) = eqns([]).
polyhedron.from_constraints([C | Cs]) = Polyhedron :-
	( lp_rational.is_false(C) ->
	  	Polyhedron = empty
	;
	  	Polyhedron0 = polyhedron.from_constraints(Cs),
		Polyhedron = polyhedron.intersection(eqns([C]), Polyhedron0)
	).

polyhedron.constraints(eqns(Constraints)) = Constraints.
polyhedron.constraints(empty_poly) = [lp_rational.false_constraint].

polyhedron.non_false_constraints(eqns(Constraints)) = Constraints.
polyhedron.non_false_constraints(empty_poly) = 
	throw("non_false_constraints/1: empty polyhedron.").

polyhedron.is_empty(empty_poly).

polyhedron.is_universe(eqns(Constraints)) :-
	list.all_true(lp_rational.nonneg_constr, Constraints).

polyhedron.optimize(_, empty_poly, empty_poly).
polyhedron.optimize(Varset, eqns(Constraints0), Result) :- 
	Constraints = simplify_constraints(Constraints0),
	( inconsistent(Varset, Constraints) ->
		Result = empty_poly
	;	
		Result = eqns(Constraints)
	).	

%------------------------------------------------------------------------------%
%
% Intersection
%

polyhedron.intersection(PolyA, PolyB, polyhedron.intersection(PolyA, PolyB)).

polyhedron.intersection(empty_poly, _) = empty_poly.
polyhedron.intersection(eqns(_), empty_poly) = empty_poly.
polyhedron.intersection(eqns(MatrixA), eqns(MatrixB)) = eqns(Constraints) :-
	Constraints0 = MatrixA ++ MatrixB,	
	restore_equalities(Constraints0, Constraints1),
	Constraints = simplify_constraints(Constraints1).

%------------------------------------------------------------------------------%
%
% Convex union.
%
% XXX At the moment this just calls convex_hull; it should actually back
% out of the convex_hull calculation if it gets too expensive (we can
% keep track of the size of the matrices during projection) and use a 
% bounding box approximation instead.
%

polyhedron.convex_union(Varset, PolyhedronA, PolyhedronB) = Polyhedron :-
	convex_union(Varset, no, PolyhedronA, PolyhedronB, Polyhedron).

polyhedron.convex_union(Varset, PolyhedronA, PolyhedronB, Polyhedron) :-
	convex_union(Varset, no, PolyhedronA, PolyhedronB, Polyhedron).

polyhedron.convex_union(Varset, MaxMatrixSize, PolyhedronA, PolyhedronB) 
		= Polyhedron :-
	convex_union(Varset, MaxMatrixSize, PolyhedronA, PolyhedronB, 
		Polyhedron).

polyhedron.convex_union(_, _, empty_poly, empty_poly, empty_poly).
polyhedron.convex_union(_, _, eqns(Constraints), empty_poly, eqns(Constraints)).
polyhedron.convex_union(_, _, empty_poly, eqns(Constraints), eqns(Constraints)).
polyhedron.convex_union(Varset, MaybeMaxSize, eqns(ConstraintsA),
		eqns(ConstraintsB), Hull) :-
	convex_hull([ConstraintsA, ConstraintsB], Hull, MaybeMaxSize, Varset).

%------------------------------------------------------------------------------%
%
% Convex hull calculation.
%

% The following transformation for computing the convex hull is described in
% the paper:
%
% F. Benoy and A. King.  Inferring Argument Size Relationships with CLPR(R).
% Logic-based Program Synthesis and Transformation,
% LNCS 1207: pp. 204-223, 1997.
%
% Further details can be found in: 
%
% F. Benoy, A. King, and F. Mesnard.
% Computing Convex Hulls with a Linear Solver
% Theory and Practice of Logic Programming 5(1&2):259-271, 2005.

:- type convex_hull_result
	--->	ok(polyhedron)
	;	aborted.	

:- type var_map == map(lp_var, lp_var).

:- type var_maps == list(var_map).

	% We introduce sigma variables into the constraints as
	% part of the transformation (See the above papers for details).
	%
:- type sigma_var == lp_var.

:- type sigma_vars == list(sigma_var).

:- type polyhedra_info 
	---> 	polyhedra_info(
			var_maps :: var_maps, 
				% There is one of these for each polyhedron.
				% It maps the original variables in the
				% constraints to the temporary variables
				% introduced by the the transformation.
				% A variable that occurs in more than one
				% polyhedron is mapped to a separate
				% temporary variable for each one.
			
			sigmas :: sigma_vars,
				% The sigma variables introduced by the
				% transformation. 
			
			poly_varset :: lp_varset
				% The varset the variables are allocated.
				% The temporary and sigma variables need
				% to be allocated from this as well in order
				% to prevent clashes when using the solver.
	).

:- type constr_info
	---> 	constr_info(
			var_map :: var_map,
				% Map from original variables
				% to new (temporary) ones.
				% There is one of these for
				% each constraint. 
			
			constr_varset :: lp_varset
  	 ).

:- pred convex_hull(list(constraints)::in, polyhedron::out, maybe(int)::in,
	lp_varset::in) is det.

convex_hull([], _, _, _) :- throw("convex_hull/4: empty list").
convex_hull([Poly], eqns(Poly), _, _).
convex_hull(Polys @ [_,_|_], ConvexHull, MaybeMaxSize, Varset0) :-
	%
	% Perform the matrix transformation from the paper
	% by Benoy and King.  Rename variables and add sigma
	% constraints as necessary.
	%
	PolyInfo0 = polyhedra_info([], [], Varset0),
	transform_polyhedra(Polys, Matrix0, PolyInfo0, PolyInfo), 
	PolyInfo = polyhedra_info(VarMaps, Sigmas, Varset),
	add_sigma_constraints(Sigmas, Matrix0, Matrix1),
	Matrix   = add_last_constraints(Matrix1, VarMaps),
	AppendValues = (func(Map, Varlist0) = Varlist :-
		Varlist = Varlist0 ++ map.values(Map)
	),
	VarsToEliminate = Sigmas ++ list.foldl(AppendValues, VarMaps, []),	
	%
	% Calculate the closure of the convex hull of the original
	% polyhedra by projecting the constraints in the transformed
	% matrix onto the original variables (ie. eliminate all the
	% sigma and temporary variables).  Since the resulting matrix
	% tends to contain a large number of redundant constraints,
	% we need to do a redundancy check after this.
	%
	lp_rational.project(VarsToEliminate, Varset, MaybeMaxSize, Matrix,
		ProjectionResult),
	(
		% XXX We should try using a bounding box first.
		ProjectionResult = aborted,
		ConvexHull = eqns(lp_rational.nonneg_box(VarsToEliminate, Matrix))
	;
		ProjectionResult = inconsistent,
		ConvexHull = empty_poly
	;
		some [!Hull] (
			ProjectionResult = ok(!:Hull),
			restore_equalities(!Hull),
			% XXX We should try removing this call to
			% simplify constraints.
			!:Hull = simplify_constraints(!.Hull),
			( remove_some_entailed_constraints(Varset, !Hull) ->
				ConvexHull = eqns(!.Hull)
			;
				ConvexHull = empty_poly	
			)	
		)
	).

:- pred transform_polyhedra(list(constraints)::in, constraints::out, 
	polyhedra_info::in, polyhedra_info::out) is det.

transform_polyhedra(Polys, Eqns, !PolyInfo) :-
	list.foldl2(transform_polyhedron, Polys, [], Eqns, !PolyInfo).

:- pred transform_polyhedron(constraints::in, constraints::in, 
	constraints::out, polyhedra_info::in, polyhedra_info::out) is det.

transform_polyhedron(Poly, Polys0, Polys, !PolyInfo) :-
	some [!Varset] (
		!.PolyInfo = polyhedra_info(VarMaps, Sigmas, !:Varset), 
		svvarset.new_var(Sigma, !Varset),
		list.map_foldl2(transform_constraint(Sigma), Poly, NewEqns,
			map.init, VarMap, !Varset),
		Polys = NewEqns ++ Polys0,
		!:PolyInfo = polyhedra_info([VarMap | VarMaps],
			[Sigma | Sigmas], !.Varset)
	).
 
	% transform_constraint: takes a constraint (with original variables) and the 
	% sigma variable to add, and returns the constraint where the original 
	% variables are substituted for new ones and where the sigma variable 
	% is included.  The map of old to new variables is updated if necessary.
	%
:- pred transform_constraint(lp_var::in, constraint::in, constraint::out, 
	var_map::in, var_map::out, lp_varset::in, lp_varset::out) is det.

transform_constraint(Sigma, !Constraint, !VarMap, !Varset) :-
	some [!Terms] (
		constraint(!.Constraint, !:Terms, Op, Const),
		list.map_foldl2(change_var, !Terms, !VarMap, !Varset),
		list.cons( Sigma - (-Const), !Terms),
		!:Constraint = constraint(!.Terms, Op, zero)
	).

	% change_var: takes a Var-Num pair with an old variable and 
	% returns one with a new variable which the old variable maps 
	% to.  Updates the map of old to new variables if necessary.
	%
:- pred change_var(lp_term::in, lp_term::out, var_map::in, var_map::out,
	lp_varset::in, lp_varset::out) is det. 

change_var(!Term, !VarMap, !Varset) :-
	some [!Var] (
		!.Term = !:Var - Coefficient,
		%
		% Have we already mapped this original variable to a new one?
		%
		( !:Var = !.VarMap ^ elem(!.Var) -> 
			true
		;
			svvarset.new_var(NewVar, !Varset),
			svmap.det_insert(!.Var, NewVar, !VarMap),
			!:Var = NewVar
		),
		!:Term = !.Var - Coefficient
	).

:- pred add_sigma_constraints(sigma_vars::in,
	constraints::in, constraints::out) is det.

add_sigma_constraints(Sigmas, !Constraints) :- 
	%
	% Add non-negativity constraints for each sigma variable.
	%
	SigmaConstraints = list.map(make_nonneg_constr, Sigmas),
	list.append(SigmaConstraints, !Constraints),
	%
	% The sum of all the sigma variables is one.
	%	
	SigmaTerms = list.map(lp_term, Sigmas),
	list.cons(constraint(SigmaTerms, (=), one), !Constraints).
		 
	% Add a constraint specifying that each variable is the sum of the
	% temporary variables to which it has been mapped.
	%
:- func add_last_constraints(constraints, var_maps) = constraints.

add_last_constraints(!.Constraints, VarMaps) = !:Constraints :-
	Keys = get_keys_from_maps(VarMaps),
	NewLastConstraints = set.filter_map(make_last_constraint(VarMaps),
		Keys),
	list.append(set.to_sorted_list(NewLastConstraints), !Constraints).
	
	% Return the set of keys in the given list of maps.
	%
:- func get_keys_from_maps(var_maps) = set(lp_var).

get_keys_from_maps(Maps) = list.foldl(get_keys_from_map, Maps, set.init).

:- func get_keys_from_map(var_map, set(lp_var)) = set(lp_var).

get_keys_from_map(Map, KeySet) = set.insert_list(KeySet, map.keys(Map)).

:- func make_last_constraint(var_maps, lp_var) = constraint is semidet.

make_last_constraint(VarMaps, OriginalVar) = Constraint :- 
	list.foldl(make_last_terms(OriginalVar), VarMaps, [], LastTerms),
	Constraint = constraint([OriginalVar - one | LastTerms], (=), zero).

:- pred make_last_terms(lp_var::in, var_map::in, lp_terms::in, lp_terms::out)
	is semidet.

make_last_terms(OriginalVar, VarMap, !Terms) :-
	NewVar = VarMap ^ elem(OriginalVar),
	list.cons(NewVar - (-one), !Terms).

%------------------------------------------------------------------------------%
%
% Approximation of polyhedron by a bounding box.
%

polyhedron.bounding_box(empty_poly, _) = empty_poly.
polyhedron.bounding_box(eqns(Constraints), Varset) = 
	eqns(lp_rational.bounding_box(Varset, Constraints)).	

%------------------------------------------------------------------------------%
%
% Widening
%

polyhedron.widen(empty_poly, empty_poly, _) = empty_poly.
polyhedron.widen(eqns(_), empty_poly, _) = throw("widen/2: empty polyhedron").
polyhedron.widen(empty_poly, eqns(_), _) = throw("widen/2: empty polyhedron"). 
polyhedron.widen(eqns(Poly1), eqns(Poly2), Varset) = eqns(WidenedEqns) :-
	WidenedEqns = list.filter(entailed(Varset, Poly2), Poly1).

%------------------------------------------------------------------------------%
%
% Projection.
%

polyhedron.project_all(Varset, Locals, Polyhedra) = 
	list.map((func(Poly0) = Poly :-
		(
			Poly0 = eqns(Constraints0),
			lp_rational.project(Locals, Varset, Constraints0,
				ProjectionResult),
			(
				ProjectionResult = aborted,
				throw("project_all/4: abort from project.")
			;
				ProjectionResult = inconsistent,
				Poly = empty_poly
			;
				ProjectionResult = ok(Constraints1),
				restore_equalities(Constraints1, Constraints),
				Poly = eqns(Constraints)
			)
		;
			Poly0 = empty_poly,
			Poly  = empty_poly
		)
	), Polyhedra).

polyhedron.project(Vars, Varset, Polyhedron0) = Polyhedron :-
	polyhedron.project(Vars, Varset, Polyhedron0, Polyhedron).

polyhedron.project(_, _, empty_poly, empty_poly).
polyhedron.project(Vars, Varset, eqns(Constraints0), Result) :-
	lp_rational.project(Vars, Varset, Constraints0, ProjectionResult),
	(
		ProjectionResult = aborted,
		throw("project/4: abort from project")
	;
		ProjectionResult = inconsistent,
		Result = empty_poly
	;
		ProjectionResult = ok(Constraints1),
		restore_equalities(Constraints1, Constraints),
		Result = eqns(Constraints)
	).
%------------------------------------------------------------------------------%
%
% Variable substitution.
%

polyhedron.substitute_vars(OldVars, NewVars, Polyhedron0) = Polyhedron :-
	Constraints0 = polyhedron.non_false_constraints(Polyhedron0),
	Constraints = lp_rational.substitute_vars(OldVars, NewVars, 
		Constraints0),	
	Polyhedron = polyhedron.from_constraints(Constraints).

polyhedron.substitute_vars(SubstMap, Polyhedron0) = Polyhedron :-
	Constraints0 = polyhedron.non_false_constraints(Polyhedron0),
	Constraints = lp_rational.substitute_vars(SubstMap, Constraints0),	
	Polyhedron = polyhedron.from_constraints(Constraints).

%------------------------------------------------------------------------------%
%
% Zeroing out variables.
%

polyhedron.zero_vars(_, empty_poly) = empty_poly.
polyhedron.zero_vars(Vars, eqns(Constraints0)) = eqns(Constraints) :-
	Constraints = lp_rational.set_vars_to_zero(Vars, Constraints0).

%------------------------------------------------------------------------------%
%
% Printing. 
%

polyhedron.write_polyhedron(empty_poly, _, !IO) :-
	io.write_string("\tEmpty\n", !IO).
polyhedron.write_polyhedron(eqns([]),   _, !IO) :-
	io.write_string("\tUniverse\n", !IO).
polyhedron.write_polyhedron(eqns(Constraints @ [_|_]), Varset, !IO) :-
	lp_rational.write_constraints(Constraints, Varset, !IO).	

%------------------------------------------------------------------------------%
:- end_module libs.polyhedron.
%------------------------------------------------------------------------------%