AUC

#define ARMA_64BIT_WORD
#define ARMA_DONT_USE_CXX11
#define ARMA_USE_MKL_ALLOC

#include <armadillo>
#include "mex.h"

using namespace arma;

typedef struct {
  rowvec TPR;
  rowvec FPR;
  double AUC;
} aucResult, *Ptr_aucResult;

struct PRCalculate {
  urowvec& classes;
  vec& scores;
  vec& thresholds;
  uvec& posLocVec;
  uvec& negLocVec;
  rowvec& TPR;
  rowvec& FPR;
  uword P;
  uword N;
  PRCalculate(urowvec& classes, vec& scores, vec& thresholds,
    uvec& posLocVec, uvec& negLocVec, rowvec& TPR,
    rowvec& FPR, uword P, uword N): classes(classes), scores(scores),
    thresholds(thresholds), posLocVec(posLocVec),
    negLocVec(negLocVec), TPR(TPR), FPR(FPR), P(P), N(N) {}

  void operator()(const blocked_range<uword>& range) const {
    for (uword i = range.begin(); i != range.end(); i++)
    {
      urowvec tmp = classes == (scores.t() >= thresholds(i));
      TPR(i) = sum(tmp.elem(posLocVec)) / (double) P;
      FPR(i) = 1 - sum(tmp.elem(negLocVec)) / (double) N;
    }
  }
};

void aucfun(urowvec& classes, vec& scores, Ptr_aucResult ptr_aucResult, bool hyperThread)
{
  uword m = classes.n_elem, P = sum(classes), N = m - P;
  vec thresholds(m+1);
  thresholds.tail(m) = sort(scores, "descend");
  thresholds(0) = thresholds(1);
  ptr_aucResult->TPR = zeros<rowvec>(m+1);
  ptr_aucResult->FPR = zeros<rowvec>(m+1);
  uvec posLocVec = find(classes), negLocVec = find(1-classes);

  int num_threads_now = hyperThread ? 2 : 1;
  mkl_set_num_threads(num_threads_now);
  PRCalculate PR_calculate(classes, scores, thresholds, posLocVec,
    negLocVec, ptr_aucResult->TPR, ptr_aucResult->FPR, P, N);
  parallel_for(blocked_range<uword>(0, m + 1), PR_calculate);
  mkl_set_num_threads(mkl_get_max_threads());
  
  ptr_aucResult->AUC = 1 - 
    sum((ptr_aucResult->TPR.tail(m) - ptr_aucResult->TPR.head(m)) % 
    (ptr_aucResult->FPR.tail(m) + ptr_aucResult->FPR.head(m))) / 2.0;
}