function x = distribpoints(P, n, yalmipoptions)
%DISTRIBPOINTS Distribute n points such that distances between them are maximized
%
% x = distribpoints(P, n)
% x = distribpoints(P, n, yalmipoptions)
%
% ---------------------------------------------------------------------------
% DESCRIPTION
% ---------------------------------------------------------------------------
% Distributes "n" points inside of a given polytope "P" such that the minimal
% distance between any of the points is maximized.
%
% NOTE! This procedure is based on mixed-integer programming and can be VERY
% slow for higher dimensions / number of points!
%
% ---------------------------------------------------------------------------
% INPUT
% ---------------------------------------------------------------------------
% P              - polytope object
% n              - number of points to distribute
% yalmipoptions  - YALMIP options as generated by sdpsettings
%
% ---------------------------------------------------------------------------
% OUTPUT                                                                                                    
% ---------------------------------------------------------------------------
% x              - coordinates of n points, stored column-wise!
%

% Copyright is with the following author(s):
%
% (C) 2005 Michal Kvasnica, Automatic Control Laboratory, ETH Zurich,
%          kvasnica@control.ee.ethz.ch
% (C) 2005 Johan Loefberg, Automatic Control Laboratory, ETH Zurich,
%          loefberg@control.ee.ethz.ch

% ---------------------------------------------------------------------------
% Legal note:
%          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 2.1 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 library; if not, write to the 
%          Free Software Foundation, Inc., 
%          59 Temple Place, Suite 330, 
%          Boston, MA  02111-1307  USA
%
% ---------------------------------------------------------------------------

error(nargchk(2,3,nargin));

global mptOptions;
if ~isstruct(mptOptions),
    mpt_error;
end

if ~isa(P,'polytope')
    error('DISTRIBPOINTS: First input argument must be a polytope.');
end
if ~isempty(P.Array),
    error('DISTRIBPOINTS: This function does not work with polytope arrays.');
end
if ~isa(n, 'double'),
    error('DISTRIBPOINTS: Second input argument must be a scalar.');
end
if any(size(n)~=1),
    error('DISTRIBPOINTS: Second input argument must be a scalar.');
end

if ~isfulldim(P)
    error('DISTRIBPOINTS: Polytope must be fully dimensional.');
end
if ~isbounded(P)
    error('DISTRIBPOINTS: Polytope must be bounded.');
end
if ~isminrep(P)
    P = reduce(P);
end
if ~isnormal(P)
    P = normalize(P);
end

dimP = dimension(P);
[A, b] = double(P);
x = sdpvar(dimP, n);

% make tighter bounds on variables
[dummy,l,u] = bounding_box(P, struct('noPolyOutput', 1));
for i = 1:size(x,2)
    bounds(x(:,i),l,u);
end

% all points have to be inside of the polytope, i.e. A*x<=b \forall x
F = set(A*x <= repmat(b, 1, n));

% define slack variables for each l_1 norm distance
% i.e. t12 <= sum(abs(x1-x2))
%      t13 <= sum(abs(x1-x3))
%      t23 <= sum(abs(x2-x3))
t = sdpvar(dimP, nchoosek(n, 2));
c = 0;
for ii = 1:n-1,
    for jj = ii+1:n,
        c = c + 1;
        F = F + set(sum(abs(x(:, ii) - x(:, jj))) > t(:, c));
    end
end

if nargin < 3,
    % use default sdpsettings
    if ~isempty(mptOptions.sdpsettings)
        yalmipoptions = mptOptions.sdpsettings;
    else
        yalmipoptions=sdpsettings('Verbose',1);
    end
end

% maximize the minimum distance between any pair of points
sol = solvesdp(F, -min(t(:)), yalmipoptions);
if sol.problem ~= 0,
    error(yalmiperror(sol.problem));
end
x = double(x);
