% **************************************************************
% * MATLAB FRAMEWORK CAP DISCREPANCY:                          *
% *                                                            *
% *   AUTHOR    : Holger Heitsch                               *
% *   INSTITUTE : WIAS Berlin                                  *
% *   (C)       : 2020                                         *
% *                                                            *
% **************************************************************

% **************************************************************
%
%   CAP DISCREPANCY MAIN FUNCTION
%
% **************************************************************

% **************************************************************
%
%   Implementation of the exact enumeration formula for the cap
%   discrepancy by H. Heitsch and R. Henrion
%
%   Input: X - array n x N of spherical points
%              (n space dimension, N number of samples)
%
%   Return: Discr - cap discrepancy
%               w - normal direction of hyperplane
%               t - shift of hyperplane
%            time - computation time
%
%   Required subroutines: cap_measure.m
%                         emp_measure.m
% 
% **************************************************************


function [Discr, w, t, time] = cdiscr (X)

    % problem size
    n = size (X,1);
    N = size (X,2);
    
    % affine extension
    Y = X;
    Y(n+1, :)  = -1 * ones (1, N);
    
    RankY = min (rank (Y), n);
     
    tm_start = clock;
    
    % initialize iteration
    pos = 1;   % index number ( 1 <= pos <= n )
    k = 1;     % loop counter ( 1 <= k   <= N )
    I = [];    % index set    ( #I <= n )

    % subcall
    if ( N > 1)
        [Discr, w, t] = cdiscr_subcall (n, N, pos, k, I, X, Y, RankY, 0, 0, 0);
    else
      if ( N == 1 )
        Discr = 1;
        w = X (:,1);
        t = 1;
      else
        % should not happen
        Discr = 0;
        w = 0;
        t = 0;
      end
    end
    
    % computing time (in seconds)
    tm_term = clock;
    tm_diff = tm_term - tm_start;
    
    sec_min  = 60;
    sec_hour = 60 * sec_min;
    sec_day  = 24 * sec_hour;
    
    time = tm_diff (3) * sec_day;
    time = time + tm_diff (4) * sec_hour;
    time = time + tm_diff (5) * sec_min;
    time = time + tm_diff (6);

end


% *************************************************************
%
% SUBROUTINE required by function: cdiscr (X)
%
% *************************************************************

function [Discr, wstar, tstar] = cdiscr_subcall (n, N, pos, k, I, X, Y, RankY, Discr, wstar, tstar)

  % set computational accuracy
  independence = 1e-24;
  gamma_tolerance = 1e-14;

  A = Y (:,I);
  
  for (j = k:N)
    
    if (pos == 1)
      % compute discrepancy
      I(pos) = j;
      w = X (:,j);
      t = 1;
      D = emp_measure (X, w, t);

      if (D > Discr && D > 1/N)
        Discr = D;
        wstar = w;
        tstar = t;
      end

    else
      % regularity check
      b = Y (:,j);
      s = A * (A \ b) - b;

      if (s'*s < independence)
        % not linearly independent
        continue;
      else
        % compute discrepancy
        I(pos) = j;
        boolean = 1;
        
        e = ones (size (I,2), 1);
        z = ( Y (:,I)' * Y (:,I) ) \ e;
        gamma = e' * z;
     
        if (gamma < 1 - gamma_tolerance)
          t = sqrt ( (1 - gamma) / gamma );
          w = ((1 + t*t) / t) * (X (:,I) * z);
          w = w / sqrt (w'*w);
        else
          if (z'*X(:,I)'*X(:,I)*z > gamma_tolerance)
             t = sqrt ( (1 - gamma) / gamma );
             w = ((1 + t*t) / t) * (X (:,I) * z);
             w = w / sqrt (w'*w);
          else
            if (size (I,2) == RankY && rank(X(:,I)) < n)
              t = 0;
              XT = X (:,I)';
              w = null (XT);
              if (size (w,2) > 0)
                w = w (:,1);
              else
                w = XT \ ones (RankY,1);
                w = w / sqrt (w'*w);
                t = XT (1,:) * w;
              end
            else
              boolean = 0;
            end
            
          end
          
        end
        
        if (boolean > 0)
          % discrepancy for (w,t), (-w,-t)
          [emp, c_emp] = emp_measure (X, w, t);
          [cap, c_cap] = cap_measure (n, t);

          D = emp - cap;
          if (D < 0) D = -D; end
          if (D > Discr)
              Discr = D;
              wstar = w;
              tstar = t;
          end

          D = c_emp - c_cap;
          if (D < 0) D = -D; end
          if (D > Discr)
              Discr = D;
              wstar = -w;
              tstar = -t;
          end
          
        end
      end
      
    end
    
    if (pos < n && j < N)
      [Discr, wstar, tstar] = cdiscr_subcall (n, N, pos+1, j+1, I, X, Y, RankY, Discr, wstar, tstar);
    end 
  
  end

end

% **************************************************************