function [x,info]=cfp_solver_maxvolume(Pts,b,ColorPartition,options)
%
% Maximum Volume heuristic algorithm to solve the colorful feasibility
% problem core cases.
%
% **********
% * Syntax *
% **********
% [x,info] = cfp_solver_maxvolume(Pts)
% [x,info] = cfp_solver_maxvolume(Pts, b)
% [x,info] = cfp_solver_maxvolume(Pts, b, ColorPartition)
% [x,info] = cfp_solver_maxvolume(Pts, b, ColorPartition, options)
%
% ***************
% * Description *
% ***************
% [x,info] = cfp_solver_maxvolume(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. The points in each color of
% Pts must contain b in their convex hulls.
%   b is a column vector representing a point in the d-Euclidean space.
% This argument is optional. In default it will be 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 iterations to solve the problem.
% info.time:   the number of seconds spent by the solver.
% info.feasible: 1 if feasible, -1 if exceeding maximum loop limit, -3 if
%                degenerated colorful simplex is encountered.
%
% *************
% * Algorithm *
% *************
% This algorithm keeps a colorful set of points, and update the set at each
% iteration. At each iteration, there is some point in the set separated
% from b by its opposite facet. These seperated points indicates potential
% replacements. Each potential replacement replaces a point in the set by a
% point of the same color and not seperated from b by the opposite facet.
% The algorithm performs one of these potential replacement such that the
% hyper volume of the convex hull of the set is maximized.
%   An important point for the implementation and outcome of this algorithm
% is: since the algorithm tries to keep the hyper volume of the simplex
% large, the colorful set of points will not be affinely dependent.


%%%%%%%%%%%%%%%%%%%%%%%%% Internal Comments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Written by Sui Huang, Advanced Optimization Lab, McMaster University,
% Hamilton, Ontario, Canada.
% 
% ************************
% * Modification History *
% ************************
% May    2005: First version.
% July   2005: Added comments.
% August 2005: Used function linsolve() to solve linear equation system.
% May    2006: Fixed a bug in generating data structure ColorMap.
%              Made vector M compatible abritary size of point sets.
% April  2007: Renamed to cfp_solver_maxvolume().
%
% ************************
% * 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. In the comment of code, T can be
%      used to denote a set of points according to the context. For
%      example, aff(T\{T{i}}) denotes the affine hyperplane generated by
%      the points indexed by T except the one of color i.
% LUM: A row vector containing coefficients for points in the colorful set
%      indicated by T. LUM represents a linear combination mapping the
%      colorful set to b. Discarding numerical error, sum(LUM)=1 should be
%      true after each iteration, and all(LUM>=0) should also be true after
%      the last iteration.
% M  : A row vector recording the relative distance of all the points to
%      their corresponding affine hyperplanes. These conditions holds after
%      each iteration: among the points of the same color i, the distance
%      between point j and aff(T\{T(j)}) is proportional to M(j); the sign 
%      of M(j) indicates which side of the hyperplane point j locates.
%

% 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^(-12);                        % Allowed numerical error.
LoopLimit = 1000000;                         % 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


% Preprocess the inputs:
%   _
%  |  b is translated to the origin
%  |  The coordinates after translation are kept in OriginPts
%  |_ The coordinates after translation and normalization are kept in Pts
OriginPts = Pts;
for cnt = 1:NumPts
    OriginPts(:,cnt) = OriginPts(:,cnt) - b;
    Pts(:,cnt) = Pts(:,cnt) - b;
    Pts(:,cnt) = Pts(:,cnt)/norm(Pts(:,cnt),2);
end
b = zeros(d,1);
% Initialized: OriginPts
% Changed: Pts b


% construct ColorMap. ColorMap is an alternative data structure to present
% which point is in which color to ease the programming.
ColorMap(NumColor).list = [];
base = 0;
for cnt=1:NumColor
    ColorMap(cnt).list = (base+1) : (base+ColorPartition(cnt));
    base = base + ColorPartition(cnt);
end
% Keep: ColorMap


% Initialize the state variables according to the initial colorful set.
% After this section, it can be proved that the state variables satisfy the
% following conditions:
%   Among the points of the same color i, the distance between point j and
% aff(T\{T(j)}) is proportional to M(j). The sign of M(j) indicates which
% side of the hyperplane point j locates.
%   NegColorList records the indices of elements in LUM that has negative
% value. NumNegColor is the length of NegColorList.
T = options.initT;
A = [Pts(:, T); ones(1,NumColor)];
trail = 0;
while (rank(A)<(d+1))
    % if T is degenerated, take some chance to find a non-degenerated one
    % by random process
    if (trail>10)
        info.iter = 0;
        info.time = etime(clock,initTime);
        info.feasible = -3;
        x = zeros(NumPts, 1);
        return
    end
    base = 0;
    for cnt = 1:NumColor
        T(cnt) = base+rand_int(1,1,[1 ColorPartition(cnt)]);
        base = base + ColorPartition(cnt);
    end
    trail = trial + 1;
end
Y = inv(A);
M = zeros(1, NumPts);
for i = 1:NumColor
    M(ColorMap(i).list) = ...
        Y(i,:) * [Pts(:,ColorMap(i).list);ones(1,ColorPartition(i))];
end
M(T) = 1;
LUM = Y * [b;1];
NegColorList = find(LUM<(-TOLERANCE))';
NumNegColor = length(NegColorList);
% Initialized: T LUM M NegColorList NumNegColor


% Initialize the counter of iterations.
LoopCnt = 0;
% Initialized: LoopCnt


% The iterations of updating the colorful set of points until the convex
% hull of the set contains b. When all elements in LUM are non-negative,
% the convex hull contains b.
while (NumNegColor > 0)
    
    % 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_maxvolume');
        disp(sbuf);
        sbuf=sprintf('Loop limit: %f   Current loop count: %f\n', ...
                     LoopLimit, LoopCnt');
        disp(sbuf);
        if (floor(LoopCnt/LoopLimit)>10)
            info.iter = LoopCnt;
            info.time = etime(clock,initTime);
            info.feasible = -1;
            x = zeros(NumPts, 1);
            return;
        end
    end
    % Changed: LoopCnt
    
    % Find which color of the colorful set to update, and which point of
    % that color to bring into the set. The general rule of updating the
    % colorful set is described at the begining of this file. In the
    % implementation aspect, two facts are utilized:
    %   (1) The points in T that have negative coefficients in LUM is
    %       seperated from b by their corresponding affine hyperplanes,
    %       which is generated by the points in T except themselves.
    %   (2) If the bottom facet of a hyper simplex does not change, its
    %       hyper volume is proportional to the distance between its top
    %       vertex and the hyperplane that contains the bottom facet.
    if ( NumNegColor > 1 )
      % If there are multiple candidates colors to be chosen from,
      % calculate and compare the colorful simplex volumes, and choose the
      % color that leading to the maximum volume.
      tmpT = T;
      maxVol = 0;
      for i = NegColorList
          % Find the local maximum volume replacement for the current 
          % color. The point indices of this color are buffered in
          % ColorMap. The search of local maximum volume replacement is
          % implemented by the search of minimum M element corresponding to
          % the points of this color, since the bottom facet of the hyper
          % simplex does not change.
          [dum, tmp] = min(M(ColorMap(i).list));
          j = ColorMap(i).list(tmp);
          % Get the volume of the local maximum replacement. This section
          % may not need to get the exact value of volume, because any
          % value proportional to the volume can be used to find the best
          % replacement by comparison.
          tmpT(i) = j;                              % The colorful set this 
                                                    % replacement leads to.
          selected_pts = Pts(1:d, tmpT);
          for cnt=2:(d+1)
              selected_pts(:,cnt) =   selected_pts(:,cnt) ...
                                    - selected_pts(:,1);
          end
          Vol = abs(det(selected_pts(:, 2:(d+1)))); % Vol is d! times of
                                                    % the actual volume.
          tmpT(i) = T(i);                           % Restore the temporary 
                                                    % variable.
          % Record the global maximum volume replacement by comparison.
          if (Vol>maxVol)
              maxVol = Vol;
              iMaxVol = i;
              jMaxVol = j;
          end 
      end
    else
      % If there is only one candidate color can be updated, decide which
      % point to bring into the colorful set, without calculating or
      % comparing the resulted colorful simplex volumes.
      iMaxVol = NegColorList;
      [dum, tmp] = min(M(ColorMap(iMaxVol).list));
      jMaxVol = ColorMap(iMaxVol).list(tmp);
    end
    T(iMaxVol) = jMaxVol;
    % Initialized: iMaxVol jMaxVol
    % Changed: T
    
    % Update the state variables according to the initial colorful set.
    % After this section, it can be proved that the state variables satisfy
    % following conditions:
    %   Among the points of the same color i, the distance between point j
    % and aff(T\{T(j)}) is proportional to M(j). The sign of M(j) indicates
    % which side of the hyperplane point j locates.
    A = [Pts(:, T); ones(1,NumColor)];
    Y = inv(A);
    M = zeros(1, NumPts);
    for i = 1:NumColor
        M(ColorMap(i).list) = ...
            Y(i,:)*[Pts(:,ColorMap(i).list); ones(1,ColorPartition(i))];
    end
    M(T) = 1;
    LUM = Y * [b;1];
    NegColorList = find(LUM<(-TOLERANCE))';
    NumNegColor = length(NegColorList);
    % Changed: LUM M NegColorList NumNegColor
    
end
% Changed: T

% Prepare the outputs. The solution x needs to be calculated from
% OriginPts, whose point coordinates are not normalized. The calculation
% does not need to handle the case of affinely dependent colorful set,
% due to the update rule in this algorithm.
info.iter = LoopCnt;
info.time = etime(clock,initTime);
info.feasible = 1;
x = zeros(NumPts, 1);
x(T) = linsolve([OriginPts(:,T); ones(1, NumColor)], [b; 1]);

return
