%-----------------------------------------------------------------------------%
% vim: ft=mercury ts=4 sw=4 et
%-----------------------------------------------------------------------------%
% Copyright (C) 2003, 2005-2006 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 maybe.
:- import_module set.

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

:- 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.compiler_util.
:- import_module libs.rat.

:- import_module pair.
:- 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) = 
    unexpected(this_file, "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([], _, _, _) :- unexpected(this_file, "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.
            %     It seems unnecessary.
            !: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, _) =
    unexpected(this_file, "widen/2: empty polyhedron").
polyhedron.widen(empty_poly, eqns(_), _) =
    unexpected(this_file, "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,
                unexpected(this_file, "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,
        unexpected(this_file, "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).    

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

:- func this_file = string.

this_file = "polyhedron.m".

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