function [hrep, scale, x0] = scaleHrep(hrep, zerotol)
% scale polytope to reasonable size and position
%
% Syntax: [hrep, scale, x0] = scaleHrep(hrep)
%   or    [hrep, scale, x0] = scaleHrep(hrep, zerotol)
%
% The Chebychev ball is the largest ball that fits into a given polytope.
%   The input polytope is scaled so that the center of the Chebychev ball
%   is the origin of the coordinate system and its radius is 1.
% The polytope remains unscaled whenever the radius of the Chebychev ball
%   is below zerotol or above 1/zerotol. For abs(r)<=zerotol the output
%   is assigned to scale=0 (flat polytope). For r=inf the output is 
%   assigned to scale=inf (unbounded polytope) and for r<(-1*zerotol) the
%   output is assigned to scale=-inf (empty polytope). Default: elkZerotol.
%   x0 is empty in unscaled cases and hrep is left unchanged.
%
% The polytope is scaled back to original by:
%   originalHrep = stretchHrep(scaledHrep, scale);
%   originalHrep = moveHrep(originalHrep, x0);
%
% See also: convertHrepToVrep, scaleHrep, reduceVrepDimension, reduceVrep

% 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/>.

%! DEV: This function could in principle use a non-uniform scaling, given
%       the procedures in reduceHrepDimension. Analogue to this procedure,
%       we could identify a limiting set of constraints (and therefore a
%       subspace for scaling) at a time.
% However, for most cases this simple size-scaling should be sufficient.
%/

%! DEV: The radius of the Chebychev ball can be limited even though the
%       polytope is not bounded. See the simple 2D example:
%          /-----------
%          |   
%          |         --> unbounded    
%          |   
%          \-----------
%       The resulting scaling does make sense, but the scaling no longer
%       guarantees to find unbounded polytopes.
%/

%% Input
if ~exist('zerotol', 'var')
    zerotol = elkZerotol;
end
nDim = size(hrep.A, 2);
nFacet = size(hrep.A, 1);
x0 = [];

%% 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);

%% find Chebychev ball
% f = [1 0 0 0]
% b = [0 h]
% A = [1 A]
problem.solver = 'linprog';
problem.options = optimset('Display', 'off', 'largeScale', 'off', ...
                           'TolFun', 0.01*zerotol, 'simplex', 'on');
problem.f = [-1 zeros(1, nDim)];
problem.Aineq = [ones(nFacet, 1) hrep.A];
problem.bineq = hrepCentered.h;
% problem.lb = [-1 -1*sqrt(1/zerotol)*ones(1, nDim)];
% problem.ub = [1  sqrt(1/zerotol)*ones(1, nDim)];
% have to add -inf to match prior implementation but avoid errors in MATLAB
%   R2012b
problem.lb = [-1*sqrt(1/zerotol) -inf(1, nDim)];
problem.up = [sqrt(1/zerotol) inf(1,nDim)];
[x, ~, exitFlag, ~, lambda] = linprog(problem);
if lambda.upper(1) ~= 0
    % this also means unbounded
    exitFlag = -3;
end
% This optimization problem should always find a feasible point since the
%   Chebychev radius is allowed to be negative to ensure a feasible point.
%   Note also that the simplex algorithm should always find this feasible
%   point automatically.

switch exitFlag
    case 1
        % solution found we have to distinguish the different radii cases
        %   while an internal point is provided:
        if x(1) < -zerotol
            % empty polytope
            hrep = [];
            scale = -inf;
            return
        elseif x(1) > 0.9/zerotol
            % we consider these cases as unbounded even if the probloem in
            %   itself is bounded. But we still should not scale this
            %   result.
            % Note that it happens that it exits with flag 1 even though
            %   the problem is quite unbounded. In this case it has hit the
            %   upper bounds.
            scale = inf;
            return
        elseif x(1) > zerotol
            % properly scalable case
            scale = x(1);
            x0 = x(2:end) + extraX0;
        else
            % this leaves the case of a flat polytope (abs(x) <= zerotol)
            %   which should not be scaled
            scale = 0;
            return
        end
    case 0
        % no solution / maximum iterations
        error('elk:polytope:numericalError', ['Scaling failed by ' ...
              'reaching maximum iterations of linprog. This should not ' ...
              'happen.']);
    case -2
        % no feasible point / polytope empty - see above comments for the
        %   formulated optimization problem
        error('elk:polytope:numericalError', ['Scaling failed by ' ...
              'not finding any feasible solution to the optimization ' ...
              'problem. This should not happen.']);
    case -3
        % polytope unbounded
        scale = inf;
        return
    otherwise
        % other errors
        error('elk:polytope:numericalError', ...
            'Scaling failed for an unknown reason. This should not happen.');
end

%% Scale polytope
hrep = moveHrep(hrep, -1*x0);
hrep = stretchHrep(hrep, 1/scale, zeros(1, nDim));
hrep = normalizeHrep(hrep);