function [hrep rot x0] = reduceHrepDimension(hrep, zerotol)
% project lower dimensional (flat) polytope into embedding space
%
% Syntax: hrep = reduceVrepDimension(hrep, zerotol)
%   or    [hrep, rot, x0] = reduceVrepDimension(hrep, zerotol)
% 
% The input polytope is reduced in its dimension by first finding the
%   linear constraints that define the embedding linear subspace. In a
%   second step the polytope is rotated such that this subspace spans the
%   last coordinate axes. Note that herewith, the position of the polytope
%   in its original space is lost so that this function returns a rotatoin
%   matrix and an offset. Note that an empty polytope returns:
%     hrep = [], rot = eye(nDim), x0 = []
%   while a zero-dimensional polytope returns:
%     hrep.A = [], hrep.h = []
%   and a proper offset vector x0.
%
% Use liftHrepDimension to restore the original polytope. Though, according
%   to the allowed tolerances, the result might have an position offset in
%   the range of zerotol.
%
% See also: liftHrepDimension, reduceVrepDimension, obtainPolytope,
%   rotatePolytope, movePolytope

% If the polytope vrep has a lower dimension r<n than the space it is
%   embedded in, this fucntion moves and rotates the polytope such that the
%   values for the last (n-r) coordinate axes are zero for all points in
%   vrep.
%
% The zerotol determines how flat a polytope must be to be recognized.
%   Points can be approximately a distance of zerotol away from the flat
%   polytope. Default: see elkZerotol.
%
% With the second syntax, the rotation matrix rot and the offset x0 is
%   returned so that the original polytope can be obtained by:
%     removedCoordinates = zeros(size(vrep, 1), size(V, 1) - size(vrep, 2));
%     vrepOriginal = [vrep, removedCoordinates]*rot';
%     vrepOriginal = moveVrep(vrepOriginal, x0);
%
% See also: reduceVrep, scaleVrep

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2012 Alexander Reinhold
%
% This program is free software: you can redistribute it and/or modify it 
%   under the terms of the GNU General Public License as published by the 
%   Free Software Foundation, either version 3 of the License, or (at your 
%   option) any later version.
%
% This program is distributed in the hope that it will be useful, but 
%   WITHOUT ANY WARRANTY; without even the implied warranty of 
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
%   General Public License for more details.
%
% You should have received a copy of the GNU General Public License along 
%   with this program.  If not, see <http://www.gnu.org/licenses/>.

%% input
if ~exist('zerotol', 'var')
    zerotol = elkZerotol;
end
% reject empty polytopes
if isempty(hrep)
    rot = [];
    x0 = [];
    return
end
% reject unbounded or zero-dim polytopes
if isempty(hrep.A)
    rot = [];
    x0 = [];
    return
end

%% algorithm
% The algorithm idea is based on some lecture notes by Komei Fukuda, Zürich 
%   2011 obtained, here:
%     http://www.ifor.math.ethz.ch/teaching/lectures/poly_comp_ss11/lecture8
%   However, the algorithm is altered quite a bit, mainly to allow for some
%   tolerance that is required in the scope of this library. 
% First, we use the lagrangian multipliers that are also returned by the 
%   constructed optimization problems. These lagrangian multipliers
%   indicate which constraints are limiting the result. How to obtain
%   these variables was not clear in the given source.
% Secondly, we do not use equality constraints. Instead, we add the
%   identified equality constraints as inequalities. This is necessary if
%   you (for example) have two opposing facets that are some distance away
%   - even though this distance is smaller than zerotol. The algorithm
%   would then select ONE of these constraints while the two constraints
%   would not exactly give the same subdomain.
%
% To make the desciption general for all steps, assume that a set of
%   constraints j from J are already found that limit the dimension by:
%      a_j * x = h_j
%   while the remainder of constraints i still gives:
%      a_i * x <= h_i
%   The optimization problem shall maximize a radius or general (possible)
%   decrease for all h_i values r so that the constraints:
%      a_i * x + r <= h_i    equaling   a_i * x <= h_i - r
%   is fulfilled together with the equality constraints.
% In the initial run where the set of equations is emtpy, r can be negative
%   which implies an empty polytope. In subsequent runs, this is not
%   possible.
% If r is positive (here, larger than zerotol) than the polytope is
%   fulldimensional when it is embedded in the subspace indicated by the
%   equality constraints.
% If r is zero (here, abs(r) < zerotol) then the polytope is degenerated
%   and the set of equations must be extended. Since only a minimal number
%   of inequalities will be marked limiting by the Lagrangian multipliers,
%   the full set should be identified by finding the linearly dependent
%   rows from [A h].
% The optimization problem might be unbounded, indicating a cone. In this
%   case, r is treated simply as positive. 
% The optimization problem cannot be infeasible since some negative value
%   for r always exists to fulfill Ax + r <= h.

%% additional scaling
% the success of the linear program is dependent on scaling and centering.
%   While polytopes are, for this task, somehow necessarily badly scaled
%   (degenerate polytopes that are reduced in dimension), there position
%   might be adopted.
extraX0 = 0.5*hrep.A'*(hrep.h./sum(hrep.A.^2, 2));
hrepCentered = moveHrep(hrep, -1*extraX0);
% (hrep.A*[0 0 0]' - hrep.h)'
% (hrepCentered.A*[0 0 0]' - hrepCentered.h)'

%% initialize loop
nDim = size(hrep.A, 2);
nFacet = size(hrep.A, 1);
inequalityFilter = true(1, nFacet);
isFirstRun = true;
keepRunning = true;
dimension = [];

problem.solver = 'linprog';
problem.options = optimset('Display', 'off', 'diagnostics', 'off', ...
    'largeScale', 'off', ...
    'TolFun', 0.01*zerotol, 'simplex', 'on');
% simplex 'on' is choosen to automatically find a feasible interior point
% simplex 'on' produces problems for decomposition validation of generic
%   lactose 7
problem.f = [-1; zeros(nDim,1)];
% problem.lb = [-1 -1*sqrt(1/zerotol)*ones(1, nDim)];
% problem.up = [1  sqrt(1/zerotol)*ones(1, nDim)];
problem.lb = [-1*sqrt(1/zerotol) -inf(1, nDim)];
problem.up = [sqrt(1/zerotol) inf(1, nDim)];
% the euqaion matrix reflects entries of the type [hrep.A hrep.h]:
equationMatrix = zeros(0, nDim+1);
nCycle = 0;

while keepRunning    
    % count the cycles: we need to know the first clycle, at least
    nCycle = nCycle + 1;
    
    %% solve optimization problem
%     inequalityFilter
%     figure,viewPolytope(struct('A', hrep.A(inequalityFilter, :), 'h', hrep.h(inequalityFilter)));
    % the ineqaulity constraints contain the original inequality
    %   constrainty, but additionally the equations are also formulated as
    %   (tight) inequality constraints to model a tolerance on the original
    %   equality constraints. This is required since from two limiting and
    %   opposing facets, only one is usually selected to reflect the
    %   constraint while x-vectors living on this constraint might (for
    %   numerical reasons) not intersect anymore with the original
    %   polytope.
    nEq = size(equationMatrix, 1);
    if nCycle == 1
        problem.Aineq = [ones(nFacet, 1), hrep.A];
        problem.bineq = [hrepCentered.h];
    else
        problem.Aineq = [ones(sum(inequalityFilter), 1), ...
                     hrep.A(inequalityFilter, :); ...
                     zeros(nEq, 1)    equationMatrix(:, 1:(end-1)); ...
                     zeros(nEq, 1) -1*equationMatrix(:, 1:(end-1))];
        problem.bineq =  [hrep.h(inequalityFilter); ...
                      equationMatrix(:, end) + zerotol; ...
                   -1*equationMatrix(:, end) + zerotol];
    end
%     max(problem.Aineq * [-0.1 0 0 0]' - problem.bineq)

%     inequalityFilter
    if any(inequalityFilter)
        [x, ~, exitFlag, ~, lambda] = linprog(problem);
%         inequalityFilter
%         x'
%         exitFlag
%         lambda.ineqlin'
%         (problem.Aineq*x - problem.bineq)'
        lambdaIneq = lambda.ineqlin(1:(end-2*nEq));
    else
        % there are no inequalities left - this happens when the polytope
        %   is zero-dimensional. The safe result here is to state that the
        %   result of the optimization problem is unbounded.
        x = [];
        exitFlag = -3;
    end
    % only the first run shall use the simplex approach, all subsequent
    %   runs shall use active-set so that we can provide an initial
    %   feasible point.
    % Note that the simplex algorithm might not find a feasible point
    %   anymore when the feasible region is too narrow.
    % x is only empty for an initially infeasible problem
    if nCycle == 1 && ~isempty(x)
        problem.options.Simplex = 'off';
        problem.x0 = [0; x(2:end)+extraX0];
        x0 = x(2:end) + extraX0;
    elseif nCycle == 1
        x0 = [];
    end
    
    % handle the exitFlag
    switch exitFlag
        case 1
            % solution found
            if x(1) < -zerotol
                % empty polytope
                reductionFlag = -1;
                keepRunning = false;
            elseif x(1) > zerotol
                % (now) full-dimensional polytope
%                 disp('full')
                reductionFlag = false;
                keepRunning = false;
            else
                % not full-dimensional, more processing required
%                 disp('reduce')
                reductionFlag = 1;
            end
        case 0
            % no solution / maximum iterations
            error('elk:polytope:numericalError', ['Dimension reduction ' ...
                'failed by reaching maximum iterations for linprog. '...
                'This should not happen.']);
        case -2
            % no feasible point this should not happen since x(1) can be used
            %   to enforce the existence of a solution
            error('elk:polytope:numericalError', ['Dimension reduction ' ...
                'failed by finding no feasible solution. '...
                'This should not happen.']);
        case -3
            % problem is unbounded - well in this case, the polytope is
            %   obviously full-dimensional even though we cannot return the
            %   inner point
            reductionFlag = 0;
            keepRunning = 0;
        otherwise
            % other errors
            error('elk:polytope:numericalError', ['Dimension reduction ' ...
                'failed for an unknown reason. '...
                'This should not happen.']);
    end
    
    if reductionFlag
%         disp('reducing')
%         lambda.ineqlin
        % index values of current inequality constraints
        inequalityIndexVector = find(inequalityFilter);
        % determine the rows in hrep.A that become equations
%         inequalityIndexVector
%         lambdaIneq
        newEqualityIndexVector = inequalityIndexVector(lambdaIneq > zerotol);
        % remove the same rows from the inequalities
        inequalityIndexVector(lambdaIneq > zerotol) = [];
        inequalityFilter(newEqualityIndexVector) = false;
        
        % get (reduced) row echolon form of the new equation rows
        [R, jb] = rref([equationMatrix; ...
                        hrep.A(newEqualityIndexVector, :) ...
                        hrep.h(newEqualityIndexVector)], 10*zerotol);
        % remove the rows that contain zero-rows to obtain new equation
        % matrix
        equationMatrix = R(1:length(jb), :);
        
        % find linearly dependent rows in the remaining [hrep.A hrep.h]
        for iRow = inequalityIndexVector
            thisRow = [hrep.A(iRow, :) hrep.h(iRow)];
            thisRow = thisRow - thisRow(jb)*equationMatrix;
            if all(abs(thisRow) < 10*zerotol)
                inequalityFilter(iRow) = false;
            end
        end
    end
end

%% compute output
% this is based on the last reductionFlag
switch reductionFlag
    case -1
        % polytope is empty
        hrep = [];
        rot = eye(nDim);
    case 0
%         disp('properly full in some dim')
        % polytope is properly full dimensional in some subspace
        if all(inequalityFilter)
            % polytope is full-dimensional as provided
            % hrep remains
            rot = eye(nDim);
        else
%             disp('full in subspace, obtaining rot')
            % polytope is full-dimensional in some subspace for which the
            %   rotation matrix is obtained by singular value
            %   decomposition:
            [u,s,v] = svd(hrep.A(~inequalityFilter, :));
            dimensionMinus = sum(diag(s) > zerotol);
            rot = v';
            % we flip the dimensions
            rot = rot(end:-1:1, :);
            % take rotated versions from A
            hrep.A = hrep.A(inequalityFilter, :)*rot';
            % cut off dimensions
            hrep.A = hrep.A(:, 1:(end-dimensionMinus));
            hrep.h = hrep.h(inequalityFilter);
            % remove zero rows
            eliminationFilter = all(abs(hrep.A) < zerotol, 2);
            hrep.A(eliminationFilter, :) = [];
            hrep.h(eliminationFilter) = [];
        end
    otherwise
end