function [x,info] = cfp_solver_randompick(Pts,b,ColorPartition, options)
%
% Random picking algorithm to solve colorful feasibility problem.
%
% **********
% * Syntax *
% **********
% [x,info] = cfp_solver_randompick(Pts)
% [x,info] = cfp_solver_randompick(Pts, b)
% [x,info] = cfp_solver_randompick(Pts, b, ColorPartition)
% [x,info] = cfp_solver_randompick(Pts, b, ColorPartition, options)
%
% ***************
% * Description *
% ***************
% [x,info]=cfp_solver_randompick(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.
%   options is a struct that holds the addition parameters. This struct is
% for the purpose of extensibility. The field options.initT indicates the
% initial colourful simplex of the algorithm. For example, 
% options.initT=[1 5 7] corresponds to the colourful simplices generated by
% the points in the first, the fifth, and the seventh columns of Pts.
%
% ********************
% * Output Arguments *
% ********************
% x is the solution of the problem. info is a struct that tells the exiting
% state, following are its members. 
% info.iter:     the number of inner iterations to solve the problem, which
%                indicates the number of colorful simplices tested in this
%                algorithm.
% info.time:     the number of seconds spent by the solver.
% info.feasible: returns 1 if the problem is feasible, -1 if the problem
%                unknow whether problem is feasible or not because of
%                exceeding maximum loop limit.
%
% *************
% * Algorithm *
% *************
% This algorithm test a random colorful set of points in each iteration
% until a valid set is met.
%       


%%%%%%%%%%%%%%%%%%%%%%%%% Internal Comments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Written by Sui Huang, Advanced Optimization Lab, McMaster University,
% Hamilton, Ontario, Canada.
% 
% ************************
% * Modification History *
% ************************
% July   2005: First version.
% August 2005: Added comments.
% April  2007: Renamed to cfp_solver_randompick().
%              Updated interface.
% %
% ************************
% * Some Local Variables *
% ************************
% T  : A row vector containing the column indices of Pts. It indicates a 
%      colorful set of points. This set is updated in each iteration until
%      it contains b in its convex hull.
%

% 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^(-10);                        % Allowed numerical error.
LoopLimit = 100000000;                       % Allowed iterations.
if (nargin==3)||((nargin==4)&&(~isfield(options,'initT')))
    options.initT = zeros(1,NumColor);
    base = 1;
    for i=1:NumColor
        options.initT(i) = base;
        base = base + ColorPartition(i);
    end
end
if (nargin<3)
    ColorPartition = (d+1)*ones(1,NumColor); % Default is (d+1) points for
                                             % each color.
    options.initT = zeros(1,NumColor);
    base = 1;
    for i=1:NumColor
        options.initT(i) = base;
        base = base + ColorPartition(i);
    end
end
if (nargin<2)
    b=zeros(d,1);                            % Default is origin.
end
% Initialized: NumPts NumColor d TOLERANCE LoopLimit ColorPartition b

% The iterations of updating the colorful set of points until the convex
% hull of the set contains b.
T = options.initT;
LoopCnt = 0;
while (1)

    % Test whether the colorful set is valid by trying to express b as a
    % convex combination.
    [tmp, recp_cond] = linsolve([Pts(:,T); ones(1, NumColor)], [b;1]);
    if ((recp_cond < TOLERANCE) && any(tmp < -TOLERANCE))
        % When the simplex is degenerate, use linear programming to check if b
        % is in the simplex.
        [tmp,LPfval,FeasibleFlag]...
            =linprog(sparse(zeros(d+1,1)),[],[],...
            [Pts(:,T); ones(1, NumColor)],[b;1],...
            zeros(d+1,1),[],[],optimset('Display','off','LargeScale','off'));
    else
        FeasibleFlag = all(tmp >= -TOLERANCE);
    end
    if (FeasibleFlag > 0)
        % feasible solution is found.
        info.iter = LoopCnt;
        info.time = etime(clock,initTime);
        info.feasible = 1;
        x = zeros(NumPts, 1);
        x(T) = tmp;
        return;
    end

    
    % Count the number of iterations.
    LoopCnt = LoopCnt + 1;
    if (LoopCnt>LoopLimit)&&(mod(LoopCnt,LoopLimit)==1)
        % Give a warning if the problem is not solved in too many
        % iterations.
        sbuf=sprintf('Warning: over loop limit in routine %s\n', ...
                     'cfp_solver_randompick');
        disp(sbuf);
        sbuf=sprintf('Loop limit: %f   Current loop count: %f\n', ...
                     LoopLimit, LoopCnt');
        disp(sbuf);
        if (floor(LoopCnt/LoopLimit)>10)
            x = zeros(NumPts, 1);
            info.iter = LoopCnt;
            info.time = etime(clock,initTime);
            info.feasible = -1;
            return
        end
    end
    % Changed: LoopCnt
    

    % Get another random colorful set
    base = 0;
    for i = 1:length(ColorPartition)
        T(i) = rand_int(1, 1, [(base+1)   (base+ColorPartition(i))]);
        base = base + ColorPartition(i);
    end
    % Changed: T
    
end

return