function [Pts,Partition,b] = cfp_gencore_lowdensity(d, equAngle, poleAngle)
%
% A random generator of colorful feasibility problem core cases with a 
% small 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_lowdensity(d)
% [Pts,Partition,b]=cfp_gencore_lowdensity(d,equAngle)
% [Pts,Partition,b]=cfp_gencore_lowdensity(d,equAngle,poleAngle)
%
% ***************
% * Description *
% ***************
% [Pts,Partition,b]=cfp_gencore_lowdensity(d,equAngle,poleAngle) randomly
% generates a colorful feasibility problem core case with a small number of
% feasible solutions. There is significant chance to get only d^2+1 
% feasible solutions.
%
% ******************
% * Input Argument *
% ******************
% The input argument d is the dimension for Euclidean space. This routine
% requires d>=2.
%   The optional arguments equAngle and poleAngle define how close are the 
% colorful points to the equator (points with the last coordinate zero) of 
% hypersphere and the south/north pole (points with the last coordinate
% 1/-1). The smaller the values, the colser. The defaults of them are both
% 0.05*pi/d. These angles must be positive and no larger than pi/4.
%
% ********************
% * 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. 
%              Named clp_gen_sminus.m
% August 2005: Added comments.
% May    2006: Handled the case of d=2.
%              Added CoreViolate parameter.
% Mar    2007  Handle the exception of d=1
% April  2007  Renamed to cfp_gencore_lowdensity.m
%              Removed CoreViolate parameter.

%% Apply default input arguments.
if ((nargin<3)||isempty(poleAngle))
    poleAngle = 0.05*pi/d;
end
if ((nargin<2)||isempty(equAngle))
    equAngle = 0.05*pi/d;
end

%% Input error check
if (d<=1)
    error('d must be at least 2');
end

%% Generate reference points for d colors on equator.
EquRef = regular_simplex(d-1);              % General a regular simplex
c = sum(EquRef,2)/(d);                      % The center of regular simplex
n = norm(EquRef(:,1)-c, 2);                 % Distance between vertices and
                                            % center
for i=1:(d)                                 % Translate and normalize the
    EquRef(:,i) = EquRef(:,i) - c;          % regular simplex vertices
    EquRef(:,i) = EquRef(:,i)/n;
end


%% Buffer some reusable calculation results.
sin_poleAngle = sin(poleAngle);
cos_poleAngle = cos(poleAngle);
sin_equAngle = sin(equAngle);
cos_equAngle = cos(equAngle);
sin_2equAngle = sin(2*equAngle);
cos_2equAngle = cos(2*equAngle);


%% Assign some control variables over the routine.
unsymmColor = rand_int(1,1,[1 (d+1)]);      % Indicate which color is the 
                                            % last color.
baseList = 0:(d+1):(d*d+d);                 % The sum of a baseList element
                                            % and an integer from [1 d+1]
                                            % covers all point indices in
                                            % Pts.
mid_norm = norm([cos_2equAngle * sum(EquRef(:,1:(d-1)),2)/(d-1); ...
                 - sin_2equAngle]);
safeAngle=asin(sin_2equAngle/mid_norm);     % The angle with equator for d
                                            % points of the last color.
RefOrder = randperm(d);                     % The order of reference points
                                            % to be used.


%% Generate the points of the first d colors.
iRef = 0;                     % Initialize the index reference point order.
Pts = zeros(d, (d+1)*(d+1));  % Allocate memory.
for base = baseList([(1:(unsymmColor-1)) ((unsymmColor+1):(d+1))]);
    % Increment the index of reference point order.
    iRef = iRef + 1;
    % Get the index of equator point that is under equator.
    iD = rand_int(1,1,[1 (d+1)]);
    % Get the index of equator point that is above equator.
    iU = rand_int(1,1,[1 d]);
    iU = iU + (iU>=iD);
    % Get the equator points.
    Pts(:, base+iD) = [ cos_2equAngle * EquRef(:,RefOrder(iRef)); ...
                        - sin_2equAngle ];
    Pts(:, base+iU) = [ - cos_equAngle * EquRef(:,RefOrder(iRef)); ...
                        sin_equAngle ];
    % Get the other points of this color.
    if (d>2)
        if (iD<iU)
            L1 = [((base+1):(base+iD-1)) ...
                  ((base+iD+1):(base+iU-1))  ...
                  ((base+iU+1):(base+d+1)) ];
            iA = rand_int(1,1,[1 (d-1)]);
            L2 = L1([1:(iA-1)  (iA+1):(d-1)]);
            iA = iA + (iA>=iD);
            iA = iA + (iA>=iU);
        else
            L1 = [((base+1):(base+iU-1)) ...
                  ((base+iU+1):(base+iD-1))  ...
                  ((base+iD+1):(base+d+1)) ];
            iA = rand_int(1,1,[1 (d-1)]);
            L2 = L1([1:(iA-1)  (iA+1):(d-1)]);
            iA = iA + (iA>=iU);
            iA = iA + (iA>=iD);
        end
        Pts(1:(d-1),L2) = randn(d-1,d-2);
        Pts(1:(d-1),base+iA) = - Pts(1:(d-1),L2) * rand_real(d-2,1,[0 1]);
        for j = L1
            Pts(1:(d-1),j) = Pts(1:(d-1),j)/norm(Pts(1:(d-1),j),2);
        end
        Pts(1:(d-1),L1) = sin_poleAngle * Pts(1:(d-1),L1);
        Pts(d,L1) = cos_poleAngle;
    elseif (d==2)
        iA = setdiff(1:(d+1), [iU iD]);
        lati = rand_real(1,1,[-poleAngle poleAngle]);
        Pts(:, base+iA) = [sin(lati); cos(lati)];
    elseif (d==1)
        error('this routine not defined for d==1');
    end
end


%% Generate the points of the last color. 
% This section does not generate exactly the configuration of reference
% [1]. It is only an approximation.
base = baseList(unsymmColor);                       % the sum of base and
                                                    % numbers in [1 (d+1)]
                                                    % covers all point
                                                    % indices of this
                                                    % color.
iD = rand_int(1,1,[1,d+1]);
L = [(base+1):(base+iD-1)  (base+iD+1):(base+d+1)]; % List of indices of
                                                    % points distribute
                                                    % above equator.
Pts(1:(d-1), L) = randn(d-1,d);
for j = L
    Pts(1:(d-1),j) = Pts(1:(d-1),j)/norm(Pts(1:(d-1),j),2);
    Pts(1:(d-1),j) = - cos(safeAngle) * Pts(1:(d-1),j);
    Pts(d,j) = sin(safeAngle);
end
W = rand_real(d,1,[0 1]);                           % Coefficients of a
for i=1:d                                           % linear combination.
    if (rand(d,1)<(1/d))                            % By statistic the
        W(i) = d*W(i);                              % chance of only 1
    end                                             % element has large
end                                                 % coefficient is
                                                    % expected to be 1/d.
Pts(:,base+iD) = -Pts(:,L)*W/sum(W);
Pts(:,base+iD) = Pts(:,base+iD)/norm(Pts(:,base+iD),2);

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