function w = svmpauc(X, y, alpha_val, beta_val, C, eps)
    % SVMpAUC - optimize pAUC in the FPR interval [alpha_val, beta_val]
    %         - learns a linear score function parameterized by 'w'
    % INPUT:
    %  <X, y> - training set <features, label>,
    %        where X is of size 'num_of_training_examples x dimension'
    %        and y is of size 'num_of_training_examples x 1'
    %  C > 0 - regularization paramter
    %  eps > 0 - tolerance parameter 'epsilon' (in the cutting plane method)
    %
    % OUTPUT:
    %  w - model learnt; w is of size 'dimension x 1'
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%This code uses Andrea Vedaldi's MATLAB API for implementing structural SVMs (http://www.robots.ox.ac.uk/~vedaldi/code/svm-struct-matlab.html);
	%this in turn is a wrapper over Thorsten Joachims' C API (http://www.cs.cornell.edu/people/tj/svm_light/svm_struct.html).
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		

      % Global variables visible to other functions    
      global m n neg_alpha_ceil neg_alpha_ceil_next salpha alpha neg_beta_floor neg_beta_ceil sbeta beta;

      %--Error checking--
      % Check the number of inputs
      if nargin ~= 6
          error('myApp:argChk', 'Wrong number of input arguments to svmpauc_special: Please give five input arguments')
      end  

      % Check if 'X' and 'y' have same number of points
      if(size(X, 1) ~= size(y, 1))
          error('Mismatch between sizes of X and y');
      end   

      % Check if 'y' is a column vector
      if(size(y, 2) ~= 1)
          error('y must be a column vector');
      end   

      % Check if 'y' is binary && 'y' has only -1 and +1
      if(size(unique(y), 1) > 2 || min(y) ~= - 1 || max(y) ~= 1)
          error('y (set of training labels) can contain only two values: -1 or 1; Also, the training set should contain both positive and negative examples');
      end   

      % Check if 0 <= alpha_val <= 1
      if(alpha_val < 0 || alpha_val > 1)
          error('alpha_val should be a value in [0, 1]');
      end  

      % Check if 0 <= beta_val <= 1
      if(beta_val < 0 || beta_val > 1)
          error('beta_val should be a value in [0, 1]');
      end   

      % Check if alpha_val < beta_val
      if(alpha_val >= beta_val)
          error('alpha_val should be lesser than beta_val');
      end  

      % Check if C > 0
      if(C <= 0)
          error('Regularization parameter C must be greater than 0');
      end  

      % Check if epsilon > 0
      if(eps <= 0)
          error('Error tolerance parameter eps must be greater than 0');
      end    
      %------------------  

      d = size(X,2); % number of features

      pos = find(y==1); % indices of positive objects
      neg = find(y==-1); % indices of positive objects

      m = size(pos,1); % number of positive objects
      n = size(neg,1); % number of negative objects

      alpha = alpha_val;
      neg_alpha_ceil = ceil(n * alpha); %number of negatives of interest ceil(n x alpha)
      neg_alpha_ceil_next =  neg_alpha_ceil + 1; %number of negatives of interest ceil(n x alpha) + 1
      salpha = neg_alpha_ceil - n*alpha; %additional factor for left out portion of the curve on the left

      beta = beta_val;
      neg_beta_floor = floor(n * beta); %number of negatives of interest floor(n x beta)
      neg_beta_ceil =  ceil(n * beta); %number of negatives of interest ceil(n x beta)
      sbeta = n*beta - neg_beta_floor; %additional factor for left out portion of the curve on the right

      % SVMStruct API requires two cell structures patterns and labels
      % containing the training set
      patterns = {}; 
      labels = {};

      % Separating the positive and negative training set examples
      Xnew=[X(pos,:);X(neg,:)];  

      patterns{1} = Xnew;  %input data 
      labels{1}  = [n*ones(m,1); m*ones(n, 1)]; %ideal output: all 1's 1^(mxn)

      % ------------------------------------------------------------------
      %                                                    Run SVM struct
      % ------------------------------------------------------------------

      parm.patterns = patterns; %Input data
      parm.labels = labels; %Output data
      parm.lossFn = @lossCB; %Loss function
      parm.constraintFn  = @constraintCB; %Find Most Violated Constraint
      parm.featureFn = @featureCB; %Joint feature map
      parm.dimension = d; %number of features

      model = svm_struct_learn([' -c ' num2str(C) ' -o 2 -v 0 -w 3 -e ' num2str(eps)], parm);
      % Call API for StructSVM with margin rescaling (-o 2), no printed output (-v 0), 
      % 1-slack algorithm - dual (-w 3), with epsilon set to eps,
      % and regularization constant set to C

      w = model.w; %Output learnt model 'w'
end

% ------------------------------------------------------------------
%                                               SVM struct callbacks
% ------------------------------------------------------------------

function delta = lossCB(param, a, abar)
    %LOSS FUNCTION: delta(a, abar)
    % a - ideal output Vs. abar - given output

      % Global variables visible to other functions    
      global m n neg_alpha_ceil neg_alpha_ceil_next salpha alpha neg_beta_floor neg_beta_ceil sbeta beta;

      aminus = abar(m+1: m+n); %  abar = [aplus; aminus]

      snpos = sort(m - aminus, 'descend');  % Sort negatives based on no. of positives below them

      if(neg_alpha_ceil == 0)
            delta = sum(snpos(1:neg_beta_floor)) + sbeta*snpos(neg_beta_ceil); 
      else
            delta = salpha*snpos(neg_alpha_ceil) + sum(snpos(neg_alpha_ceil_next:neg_beta_floor)) + sbeta*snpos(neg_beta_ceil); 
      end
      % 1 - pAUC(0, beta), where  
      %     pAUC(0, beta) = 
      %             (fraction of area between neg_alpha_ceil and neg_alpha_ceil_next) 
      %                         + (AUC between neg_alpha_ceil_next and neg_beta_floor) + (fraction of area between neg_beta_floor + neg_beta_ceil)
      delta = delta / m / (n * (beta - alpha)); %Normalize
end

function psi = featureCB(param, x, a)
    %JOINT FEATURE MAP: psi(x, a)
    % x - input, y - output

      % Global variables visible to other functions    
      global m n alpha beta;

      Xplus = x(1:m, :); % Positive instances
      Xminus = x(m+1:m+n, :); % Negative instances

      aplus = a(1: m); %  a = [aplus; aminus]
      aminus = a(m+1: m+n); 

      % Joint feature representation: psi = \sum_{i = 1}^m aplus(i)*Xplus(i) -
      %                                      \sum_{j = 1}^n aminus(j)*Xminus(j)
      psi = (aplus-n/2)'*Xplus - (aminus-m/2)'*Xminus; 
      %'-n/2' and '-m/2' ensure that results dont change with previous representation

      psi = psi / m / (n * (beta - alpha));
      psi = sparse(psi'); % Sparse representation of joint feature map (required by API)
end

function ahat = constraintCB(param, model, x, a)
    % FIND MOST VIOLATED CONSTRAINT
    %   margin rescaling: argmax_a delta(ai, a) + <psi(x,a), w>

      % Global variables visible to other functions
      global neg_alpha_ceil;
    
      if(neg_alpha_ceil == 0)
          ahat = constraintCB_special(model, x); %Call special case for [0, beta]
      else
          ahat = constraintCB_general(model, x); %Call general case for [alpha, beta]
      end
end

function ahat = constraintCB_special(model, x)
    % FIND MOST VIOLATED CONSTRAINT
    %   margin rescaling: argmax_a delta(ai, a) + <psi(x,a), w>
    %   for special case [0, beta]

      % Global variables visible to other functions    
      global m n neg_beta_floor neg_beta_ceil sbeta;

      w = model.w; %Given 'w'
      wx = x*w; %Score function w.x

      [wpos, pos_ind] = sort(wx(1:m),'descend'); % Sort positives based on w.x
      [wneg, neg_ind] = sort(wx(m+1:m+n),'descend'); % Sort negatives based on w.x

      % Score function for positive instances: fplus(i) = w^T x^+_i
      fplus = wx(1:m);

      % Score function for negative instances: 
      %         fminus(j) = w^T x^-_j + 1        if (j)_a \in {1, ..., j_\beta}
      %                   = w^T x^-_j + s_\beta  if (j)_a = j_\beta + 1
      %                   = w^T x^-_j             o/w  
      fminus = wx(m+1: m+n);
      fminus(neg_ind(1:neg_beta_floor)) = fminus(neg_ind(1:neg_beta_floor)) + 1;
      if(neg_beta_floor < neg_beta_ceil) %if(floor <> ceil)
          fminus(neg_ind(neg_beta_ceil)) = fminus(neg_ind(neg_beta_ceil)) + sbeta;
      end

      %Sort instances according to f = [fminus; fplus]; %negatives first, to ensure ties are penalized
      [fsort, fsort_indices] = sort([fminus; fplus], 'descend');

      %Labels of sorted instances
      labels = [-ones(n, 1); ones(m, 1)]; 
      slabels = labels(fsort_indices);

      %Convert to output representation
      ahat = convert2output(slabels, pos_ind, neg_ind);
end

function ahat = constraintCB_general(model, x)
    % FIND MOST VIOLATED CONSTRAINT
    %   margin rescaling: argmax_a delta(ai, a) + <psi(x,a), w>
    %   for general case [alpha, beta] 

      % Global variables visible to other functions      
      global m n neg_alpha_ceil neg_alpha_ceil_next neg_beta_floor neg_beta_ceil sbeta;

      w = model.w; %Given 'w'
      wx = x*w; %Score function w.x

      [wpos, pos_ind] = sort(wx(1:m),'descend'); % Sort positives based on w.x
      [wneg, neg_ind] = sort(wx(m+1:m+n),'descend'); % Sort negatives based on w.x

      % Score function for positive instances: fplus(i) = w^T x^+_i
      fplus = wx(1:m);

      %%%%%%%%%Optimize over r_i \in {1, ..., j_\alpha - 1}%%%%%%%%%

      if(neg_alpha_ceil > 0)
          % Score function for negative instances in range [(1)_a, ..., (j_\alpha - 1)_a];
          %         fminus(j) = w^T x^-_j
          fminus1 = wx(m+1: m+n);
	  fminus1 = fminus1(neg_ind(1:neg_alpha_ceil-1));

          %Sort instances according to f = [fminus; fplus]; %negatives first, to ensure ties are penalized
          [fsort1, fsort_indices1] = sort([fminus1; fplus], 'descend');

          %Labels corresponding to ordering zbar^1
          labels1 = [-ones(neg_alpha_ceil - 1, 1); ones(m, 1)];
          slabels1 = [labels1(fsort_indices1); -ones(n - neg_alpha_ceil + 1, 1)];

          %Compute objective values
          Q1 = computeObjective(slabels1, wpos, wneg);
      else
          Q1 = -ones(m, 1); %Make objective values -1 (low value) if j_\alpha = 0
      end
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      %%%%%%%%%Optimize over r_i \in {j_\alpha}%%%%%%%%%
      % Labels corresponding to ordering zbar^2
      slabels2 = [-ones(neg_alpha_ceil, 1); ones(m, 1); -ones(n - neg_alpha_ceil, 1)];

      % Compute objective values
      Q2 = computeObjective(slabels2, wpos, wneg);  
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      %%%%%%%%%Optimize over r_i \in {j_\beta + 2, ..., n}%%%%%%%%%%

      % Score function for negative instances in range [(j_\alpha + 2)_a, ..., (n)_a];
      %         fminus(j) = w^T x^-_j + 1        if (j)_a \in {j_\alpha + 2, ..., j_\beta}
      %                   = w^T x^-_j + s_\beta  if (j)_a = j_\beta + 1
      %                   = w^T x^-_j            if (j)_a \in {j_\beta + 2, ..., n} 
      fminus3 = wx(m+1: m+n);
      fminus3(neg_ind(neg_alpha_ceil_next+1:neg_beta_floor)) = fminus3(neg_ind(neg_alpha_ceil_next+1:neg_beta_floor)) + 1;
      if(neg_beta_floor < neg_beta_ceil) %if(floor <> ceil)
          fminus3(neg_ind(neg_beta_ceil)) = fminus3(neg_ind(neg_beta_ceil)) + sbeta;
      end
      fminus3 = fminus3(neg_ind(neg_alpha_ceil_next+1:n));
      % Sort instances according to f = [fminus3; fplus]; %negatives first, to ensure ties are penalized
      [fsort3, fsort_indices3] = sort([fminus3; fplus], 'descend');

      % Labels corresponding to ordering zbar^3
      labels3 = [-ones(n - neg_alpha_ceil_next, 1); ones(m, 1)];
      slabels3 = [-ones(neg_alpha_ceil_next, 1); labels3(fsort_indices3)];

      % Compute objective values
      Q3 = computeObjective(slabels3, wpos, wneg);   
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      % Construct optimal ordering from Q1, Q2 and Q3
      [max_obj, max_ind] = max([Q1, Q2, Q3], [], 2);  

      num_pos_in_part1 = sum(max_ind == 1);
      num_pos_in_part2 = sum(max_ind == 2);
      num_pos_in_part3 = sum(max_ind == 3);

      [fsort, fsort_indices3] = sort([fminus3; wpos(num_pos_in_part1+num_pos_in_part2+1:m)], 'descend');
      labels3 = [-ones(n - neg_alpha_ceil_next, 1); ones(num_pos_in_part3, 1)];

      % If alpha is non-zero
      if(neg_alpha_ceil > 0)
          [fsort, fsort_indices1] = sort([fminus1; wpos(1:num_pos_in_part1)], 'descend');
          labels1 = [-ones(neg_alpha_ceil - 1, 1); ones(num_pos_in_part1, 1)];

          slabels = [labels1(fsort_indices1); -1; ones(num_pos_in_part2, 1); -1; labels3(fsort_indices3)];
      else
          slabels = [ones(num_pos_in_part2, 1); -1; labels3(fsort_indices3)];
      end

      % Convert ordering to output representation
      ahat = convert2output(slabels, pos_ind, neg_ind);    
end

function Q = computeObjective(labels, wpos, wneg)
    % COMPUTE OBJECTIVE VALUE FOR GIVEN OUTPUT

        % Global variables visible to other functions
        global m neg_alpha_ceil neg_alpha_ceil_next salpha neg_beta_floor neg_beta_ceil sbeta;
    
        % Maintain number of positive and negative examples seen so far
        i = 0; % Number of positive and negative examples seen so far
        j = 0; % Number of positive and negative examples seen so far

        neg1 = 0; % No. of negatives encountered till now
        neg2 = 0; % No. of negatives encountered till now in the range [(j_alpha + 1): (j_beta)]
        flag_alpha = 0; % Has (j_alpha)^th negative been encountered?
        flag_beta = 0; % Has (j_beta + 1)^th negative?
        s = 0; % Sum of scores of negatives encountered till now

        b = zeros(m, 1); % b[i] = no. of negatives in the range [(j_alpha + 1): (j_beta)] above the (i)^th positive
        c = zeros(m, 1); % b[i] = no. of negatives above the (i)^th positive
        d = zeros(m, 1); % c[i] = sum of scores of negatives above the (i)^th positive
        rho_alpha = zeros(m, 1); % rho_alpha(i) = is (j_alpha)^th negative above (i)^th positive?
        rho_beta = zeros(m, 1); % rho_beta(i) = is (j_beta + 1)^th negative above (i)^th positive?

        for l = labels'
            if(l == 1) %Positive label
                i = i + 1;
                b((i)) = neg1;
                c((i)) = neg2;
                d((i)) = s;
                rho_alpha((i)) = flag_alpha;            
                rho_beta((i)) = flag_beta;
            else %Negative label
                j = j + 1;
                if(neg_alpha_ceil > 0 && j == neg_alpha_ceil)
                    flag_alpha = 1;
                elseif(sum(j == neg_alpha_ceil_next:neg_beta_floor))
                    neg1 = neg1 + 1;
                elseif(j == neg_beta_ceil)
                    flag_beta = 1;
                end           
                neg2 = neg2 + 1;
                s = s + wneg(j); 
            end
        end

        %Compute objective
        Q = (salpha*rho_alpha + b + sbeta*rho_beta) - c.*wpos + d;
end

function a = convert2output(labels, pos_ind, neg_ind)    
    % CONVERT GIVEN OUTPUT TO A REPRESENTATION OF SIZE mxn
    
      % Global variables visible to other functions
      global m n;

      % Maintain number of positive and negative examples seen so far
      i = 0; % Number of positive and negative examples seen so far
      j = 0; % Number of positive and negative examples seen so far

      % a = [aplus; aminus]
      aplus = zeros(m, 1);
      aminus = zeros(n, 1);

      pos = 0; % No. of positives encountered
      neg = 0; % No. of negatives encountered

      for l = labels'
          if(l == 1) %Positive label
              i = i + 1;
              pos = pos + 1;
              aplus(pos_ind(i)) = n - neg;
          else %Negative label
              j = j + 1;
              neg = neg + 1;
              aminus(neg_ind(j)) = pos;
          end
      end

      a = [aplus; aminus];
end
