function [x,info] = cfp_solver_sdprelax(Pts, b, ColorPartition)
%
% The semidefinite relaxation solver to solve colorful feasibility problem.
% This function relaxes the colorful linear feasibility problem to a
% semidefinite programming problem and tries to obtain a feasible solution
% of colorful feasibility problem.
%
% **********
% * Syntax *
% **********
% [x,info] = cfp_solver_sdprelax(Pts)
% [x,info] = cfp_solver_sdprelax(Pts, b)
% [x,info] = cfp_solver_sdprelax(Pts, b, ColorPartition)
%
% ***************
% * Description *
% ***************
% [x,info]=clp_solver_sdprelax(Pts,b,ColorPartition) 
% solves for x satisfying the constraints:
%        _
%       |  Pts*x = b
%       |  if x(i) and x(j) have the same color, then x(i)*x(j)=0
%       |  sum(x)=1
%       |_ x>=0
%
% *******************
% * Input 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.
%   b is a column vector representing a point in the d-Euclidean space.
% This argument is optional. In default it is the origin.
%   ColorPartition is a row vector of length (d+1). Each element is an
% integer, 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 argument is
% optional. In default it assumes (d+1) points in each of the (d+1) colors.
% If the problem is not in the default case, user must provide this
% argument.
%
% ********************
% * Output Arguments *
% ********************
% x is the solution of the problem.
% info is a struct that tells the exiting state, following are its members.
% info.feasible:   1 if x is feasible for the input colorful feasibility
%                  problem, 0 if infeasible, -1 if unknown. If this member
%                  is 1, then the info.relaxfeasi relax must be 1.
% info.relaxfeasi: 1 if the relaxed problem is feasible, 0 if not. If this
%                  member is 0, then info.feasible must be 0.
% info.time:       the number of seconds spent by the solver.
%
% *************
% * Algorithm *
% *************
% We use Q to notate a square matrix such that Q(i,j)=1 if x(i) and x(j)
% are for the same color, and Q(i,j)=0 otherwise, X to notate a square
% matrix of proper size. Then the colorful feasibility problem can be
% formulated as the nondefinite QP:
%
%         min x'Qx
%   such that     Px = 0
%             sum(x) = 0
%                 x >= 0
%
% and it can be relaxed to the SDP:
%
%         min QX
%   such that     Px = 0
%             sum(x) = 0
%                 x >= 0
%       [1 x'; x X] is PSD and symmetric
%                 X >= 0
%
% If the SDP optimal objective value is grater than zero, then the colorful
% feasibility problem is not feasible. If the SDP optimal objective value
% is zero, then the colorful feasibility may or may not be feasible.
%
% ****************
% * Dependencies *
% ****************
% This routine following toolboxes or libraries:
%     yalmip library at http://control.ee.ethz.ch/~joloef/yalmip.php
%     sedumi at http://sedumi.mcmaster.ca
% 

%%%%%%%%%%%%%%%%%%%%%%%%% Internal Comments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Written by Sui Huang, Advanced Optimization Lab, McMaster University,
% Hamilton, Ontario, Canada.
% 
% ************************
% * Modification History *
% ************************
% April 2007: First version by Sui Huang
% 
% ************************
% * Some Local Variables *
% ************************
% F      : an object obtained by yalmip toolbox to hold all the constraints
%          of an optimization problem.
% xvar   : An object obtained by yalmip toolbox to hold optimization
%          variables.
% bigXvar: an object obtained by yalmip toolbox to hold optimization
%          variables.
%

% Initialize the time counter.
initTime = clock;
% Initialized: initTime

% Assign some global control numbers and apply the default arguments.
[d NumPts]=size(Pts);                        % Space dimension and number 
                                             % of points.
NumColor = d+1;                              % Number of colors.
TOLERANCE = 10^(-8);                         % Allowed numerical error.
if (nargin<3)
    ColorPartition = (d+1)*ones(1,NumColor); % Default is (d+1) points for 
                                             % each color.
end
if (nargin<2)
    b=zeros(d,1);                            % Default is origin.
end
% Initialized: NumPts NumColor d TOLERANCE ColorPartition b

% Initialize the info struct for output.
info.feasible = 1;
info.relaxfeasi = 1;
info.time = 0;
% Initialized: initTime

% Formulate and solve the SDP.
Q=[];
for i=1:length(ColorPartition)
    Q=blkdiag(Q,~eye(ColorPartition(i)));
end
xvar = sdpvar(NumPts,1);
bigXvar = sdpvar(NumPts, NumPts, 'symmetric');
F1 = set(Pts*xvar==b, 'x elements are cofficients to get b');
F2 = set(sum(xvar)==1, 'x elements sum up to 1');
F3 = set(xvar>=0, 'elements in x are non-negative');
F4 = set([1 xvar'; xvar bigXvar]>=0, '[1 x; x bigX] is PSD');
F5 = set(sum(sum(bigXvar))==1, 'bigX elements sum up to 1');
F6 = set(bigXvar(:)>=0, 'bigX elements are non-negatives');
F = F1+F2+F3+F4+F5+F6;
solveroption = sdpsettings('solver','sedumi', 'verbose',0, ...
                           'sedumi.stepdif',2, 'sedumi.sdp',1, ...
                           'sedumi.cg.qprec', 1);
diagnostics = solvesdp(F, sum(sum(Q.*bigXvar)), solveroption);
x = double(xvar);
bigX = double(bigXvar);
% initialized: diagnostics x

if (diagnostics.problem==1)
    % SDP is infeasible
    info.relaxfeasi = 0;
    info.feasible = 0;
    info.time = etime(clock,initTime);
    return;
else
    % SDP is feasible
    if (bigX(:)'*Q(:)>TOLERANCE)
        % colorful feasibility problem is infeasible
        info.relaxfeasi = 1;
        info.feasible = 0;
        info.time = etime(clock,initTime);
    elseif iscolorful(x,ColorPartition,TOLERANCE)
        % fortunately the solution of SDP is also a solution of colorful
        % feasibility problem
        info.relaxfeasi = 1;
        info.feasible = 1;
        info.time = etime(clock,initTime);
    else
        % the solution of SDP is not a solution of colorful feasibility
        % problem
        info.relaxfeasi = 1;
        info.feasible = -1;
        info.time = etime(clock,initTime);
    end
end

return

function flag = iscolorful(x, ColorPartition, TOLERANCE)
base = 0;
posFlag = x>TOLERANCE;
for i=1:length(ColorPartition)
    if (sum(posFlag((base+1):(base+ColorPartition(i))))>1)
        flag = 0;
        return;
    end
    base = base + ColorPartition(i);
end
flag = 1;
return