function [Pts, Partition, b] = cfp_gencore_middensity(d, r)
%
% A random generator of colorful feasibility problem core cases with medium
% number of feasible solutions. These random problems are essentially some
% particular structured configurations up to perturbation and shuffling
% of the order of points.
%
% **********
% * Syntax *
% **********
% [Pts, Partition, b] = cfp_gencore_middensity(d)
% [Pts, Partition, b] = cfp_gencore_middensity(d, r)
%
% ***************
% * Description *
% ***************
% [Pts,Partition,b]=cfp_gencore_middensity(d,r) randomly generates a
% colorful feasibility problem core case with medium number of feasible
% solutions. It must be r<1/d such that the problem has (d+1)! feasible
% solutions.
%
% ******************
% * Input Argument *
% ******************
% The input argument d is the dimension for Euclidean space.
% The optional argument r defines how close do d points of the different
% colors cluster. The default is 1/d to garantee there exist (d+1)!
% feasible solutions and the problem is solvable.
%
% ********************
% * Output Arguments *
% ********************
% Pts is a matrix storing the coordinates of points. Each column of Pts
% stores the coordinate of one point. The number of rows is d, which is
% the number of dimensions of the Euclidean space. The points in each color
% of Pts must contain b in their convex hulls. All the colorful points are
% unit vectors. Scaling the colorful points does not affect whether a
% convex hull generated by them contain the origin, and b is the origin.
% Hence, all the colorful points can be scaled to unit vectors.
%   b is a column vector representing a point, and it is always generated
% as the origin. This output argument is for future potential update.
%   Partition is a row vector of length (d+1). Each element is an
% integer no less than (d+1), specifying the number of points in a color.
% For example [3 4 3] tells that the first 3 points in Pts are in the first
% color, the following 4 points are in the second color, and so on. This
% function only assigns (d+1) points to each color, since in practice a
% subset of size (d+1) containing b in the convex hull can always be found
% when there are more than (d+1) points. Partition is an output
% argument for future potential update.
%
% **************
% * References *
% **************
% [1] A Deza, S Huang, T Stephen and T Terlaky, "Colourful Simplicial
% Depth", Discrete & Computational Geometry, Volume 35, 597-615.

%%%%%%%%%%%%%%%%%%%%%%%%% Internal Comments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Written by Sui Huang, Advanced Optimization Lab, McMaster University,
% Hamilton, Ontario, Canada.
% 
% ************************
% * Modification History *
% ************************
% May    2005: First version.
% August 2005: Added comments.
% May    2006: Added CoreViolate parameter.
% April  2007: Removed CoreViolate parameter.
%              Renamed to cfp_gencore_middensity().

%% Apply the default input arguments.
if ((nargin<2)||isempty(r))
    r = 1/d;
end

%% Generate regular simplex points
% to guide the distribution of colorful points.
RegularPts = regular_simplex(d);
c = sum(RegularPts,2)/(d+1);
n = norm(RegularPts(:,1)-c, 2);
for i=1:(d+1)
    RegularPts(:,i) = RegularPts(:,i) - c;
    RegularPts(:,i) = RegularPts(:,i)/n;
end

%% Allocate memory for Pts.
Pts = zeros(d, (d+1)*(d+1));

%% Decide the ordering of points
% OrderMatrix is use to randomly shuffle the order of points in each color.
% Each column of OrderMatrix stores a random permutation of points inside
% one color.
OrderMatrix = zeros(d+1, d+1);
for i = 0:d
    OrderMatrix(:,i+1) = randperm(d+1)' + i*(d+1);
end

%% Generate the points of each color.
for i = 1:(d+1)
    % Buffer a vertex of the regular simplex.
    c = RegularPts(:,i);
    % Generate d points of different colors around a vertex of the regular
    % simplex.
    for j = OrderMatrix(i,:)
        t = randn(d,1);
        t = c + r*t/norm(t,2);
        Pts(:,j) = t/norm(t,2);
    end
end

%% Generate b as the origin.
b = zeros(d,1);

%% Partition tells that each color has (d+1) points.
Partition = (d+1)*ones(1,(d+1));

return
