error.cpp 18 KB
Newer Older
André Anjos's avatar
André Anjos committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14
/**
 * @date Wed Apr 20 08:02:30 2011 +0200
 * @author Andre Anjos <andre.anjos@idiap.ch>
 *
 * @brief Implements the error evaluation routines
 *
 * Copyright (C) Idiap Research Institute, Martigny, Switzerland
 */

#include <stdexcept>
#include <algorithm>
#include <limits>
#include <boost/format.hpp>

André Anjos's avatar
André Anjos committed
15
#include <bob.core/assert.h>
André Anjos's avatar
André Anjos committed
16
#include <bob.core/cast.h>
17 18
#include <bob.core/array_copy.h>
#include <bob.core/array_sort.h>
André Anjos's avatar
André Anjos committed
19

20 21
#include <bob.math/pavx.h>
#include <bob.math/linsolve.h>
André Anjos's avatar
André Anjos committed
22

23
#include "error.h"
24

25
template <typename T>
26 27 28 29 30
static void sort(const blitz::Array<T,1>& a, blitz::Array<T,1>& b, bool isSorted){
  if (isSorted){
    b.reference(a);
  } else {
    bob::core::array::ccopy(a,b);
31
    bob::core::array::sort<T>(b);
32 33 34
  }
}

André Anjos's avatar
André Anjos committed
35 36
std::pair<double, double> bob::measure::farfrr(const blitz::Array<double,1>& negatives,
    const blitz::Array<double,1>& positives, double threshold) {
37 38
  if (!negatives.size()) throw std::runtime_error("Cannot compute FAR when no negatives are given");
  if (!positives.size()) throw std::runtime_error("Cannot compute FRR when no positives are given");
André Anjos's avatar
André Anjos committed
39 40 41 42 43 44 45 46 47 48
  blitz::sizeType total_negatives = negatives.extent(blitz::firstDim);
  blitz::sizeType total_positives = positives.extent(blitz::firstDim);
  blitz::sizeType false_accepts = blitz::count(negatives >= threshold);
  blitz::sizeType false_rejects = blitz::count(positives < threshold);
  return std::make_pair(false_accepts/(double)total_negatives,
      false_rejects/(double)total_positives);
}

std::pair<double, double> bob::measure::precision_recall(const blitz::Array<double,1>& negatives,
    const blitz::Array<double,1>& positives, double threshold) {
49
  if (!negatives.size() || !positives.size()) throw std::runtime_error("Cannot compute precision or recall when no positives or no negatives are given");
André Anjos's avatar
André Anjos committed
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76
  blitz::sizeType total_positives = positives.extent(blitz::firstDim);
  blitz::sizeType false_positives = blitz::count(negatives >= threshold);
  blitz::sizeType true_positives = blitz::count(positives >= threshold);
  blitz::sizeType total_classified_positives = true_positives + false_positives;
  if (!total_classified_positives) total_classified_positives = 1; //avoids division by zero
  if (!total_positives) total_positives = 1; //avoids division by zero
  return std::make_pair(true_positives/(double)(total_classified_positives),
      true_positives/(double)(total_positives));
}


double bob::measure::f_score(const blitz::Array<double,1>& negatives,
    const blitz::Array<double,1>& positives, double threshold, double weight) {
  std::pair<double, double> ratios =
      bob::measure::precision_recall(negatives, positives, threshold);
  double precision = ratios.first;
  double recall = ratios.second;
  if (weight <= 0) weight = 1;
  if (precision == 0 && recall == 0)
    return 0;
  return (1 + weight*weight) * precision * recall / (weight * weight * precision + recall);
}

double eer_predicate(double far, double frr) {
  return std::abs(far - frr);
}

77 78 79 80 81
double bob::measure::eerThreshold(const blitz::Array<double,1>& negatives, const blitz::Array<double,1>& positives, bool isSorted) {
  blitz::Array<double,1> neg, pos;
  sort(negatives, neg, isSorted);
  sort(positives, pos, isSorted);
  return bob::measure::minimizingThreshold(neg, pos, eer_predicate);
André Anjos's avatar
André Anjos committed
82 83
}

84
double bob::measure::eerRocch(const blitz::Array<double,1>& negatives, const blitz::Array<double,1>& positives) {
André Anjos's avatar
André Anjos committed
85 86 87
  return bob::measure::rocch2eer(bob::measure::rocch(negatives, positives));
}

88
double bob::measure::farThreshold(const blitz::Array<double,1>& negatives, const blitz::Array<double,1>&, double far_value, bool isSorted) {
André Anjos's avatar
André Anjos committed
89 90 91 92 93 94 95 96 97 98
  // check the parameters are valid
  if (far_value < 0. || far_value > 1.) {
    boost::format m("the argument for `far_value' cannot take the value %f - the value must be in the interval [0.,1.]");
    m % far_value;
    throw std::runtime_error(m.str());
  }
  if (negatives.size() < 2) {
    throw std::runtime_error("the number of negative scores must be at least 2");
  }

99 100 101
  // sort the array, if necessary
  blitz::Array<double,1> neg;
  sort(negatives, neg, isSorted);
André Anjos's avatar
André Anjos committed
102 103 104

  // compute position of the threshold
  double crr = 1.-far_value; // (Correct Rejection Rate; = 1 - FAR)
105
  double crr_index = crr * neg.extent(0);
André Anjos's avatar
André Anjos committed
106
  // compute the index above the current CRR value
107
  int index = std::min((int)std::floor(crr_index), neg.extent(0)-1);
André Anjos's avatar
André Anjos committed
108 109

  // correct index if we have multiple score values at the requested position
110
  while (index && neg(index) == neg(index-1)) --index;
André Anjos's avatar
André Anjos committed
111

112
  // we compute a correction term to assure that we are in the middle of two cases
André Anjos's avatar
André Anjos committed
113 114 115
  double correction;
  if (index){
    // assure that we are in the middle of two cases
116
    correction = 0.5 * (neg(index) - neg(index-1));
André Anjos's avatar
André Anjos committed
117 118
  } else {
    // add an overall correction term
119
    correction = 0.5 * (neg(neg.extent(0)-1) - neg(0)) / neg.extent(0);
André Anjos's avatar
André Anjos committed
120 121
  }

122
  return neg(index) - correction;
André Anjos's avatar
André Anjos committed
123 124
}

125
double bob::measure::frrThreshold(const blitz::Array<double,1>&, const blitz::Array<double,1>& positives, double frr_value, bool isSorted) {
André Anjos's avatar
André Anjos committed
126 127 128 129 130 131 132 133 134 135 136

  // check the parameters are valid
  if (frr_value < 0. || frr_value > 1.) {
    boost::format m("the argument for `frr_value' cannot take the value %f - the value must be in the interval [0.,1.]");
    m % frr_value;
    throw std::runtime_error(m.str());
  }
  if (positives.size() < 2) {
    throw std::runtime_error("the number of positive scores must be at least 2");
  }

137 138
  // sort positive scores descendantly, if necessary
  blitz::Array<double,1> pos;
139
  sort(positives, pos, isSorted);
André Anjos's avatar
André Anjos committed
140 141

  // compute position of the threshold
142 143 144
  double frr_index = frr_value * pos.extent(0);
  // compute the index above the current CAR value
  int index = std::min((int)std::ceil(frr_index), pos.extent(0)-1);
André Anjos's avatar
André Anjos committed
145 146

  // correct index if we have multiple score values at the requested position
147
  while (index < pos.extent(0)-1 && pos(index) == pos(index+1)) ++index;
André Anjos's avatar
André Anjos committed
148 149 150

  // we compute a correction term to assure that we are in the middle of two cases
  double correction;
151
  if (index < pos.extent(0)-1){
André Anjos's avatar
André Anjos committed
152
    // assure that we are in the middle of two cases
153
    correction = 0.5 * (pos(index+1) - pos(index));
André Anjos's avatar
André Anjos committed
154 155
  } else {
    // add an overall correction term
156
    correction = 0.5 * (pos(pos.extent(0)-1) - pos(0)) / pos.extent(0);
André Anjos's avatar
André Anjos committed
157 158
  }

159
  return pos(index) + correction;
André Anjos's avatar
André Anjos committed
160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181
}

/**
 * Provides a functor predicate for weighted error calculation
 */
class weighted_error {

  double m_weight; ///< The weighting factor

  public: //api

  weighted_error(double weight): m_weight(weight) {
    if (weight > 1.0) m_weight = 1.0;
    if (weight < 0.0) m_weight = 0.0;
  }

  inline double operator() (double far, double frr) const {
    return (m_weight*far) + ((1.0-m_weight)*frr);
  }

};

182 183 184 185
double bob::measure::minWeightedErrorRateThreshold(const blitz::Array<double,1>& negatives, const blitz::Array<double,1>& positives, double cost, bool isSorted) {
  blitz::Array<double,1> neg, pos;
  sort(negatives, neg, isSorted);
  sort(positives, pos, isSorted);
André Anjos's avatar
André Anjos committed
186
  weighted_error predicate(cost);
187
  return bob::measure::minimizingThreshold(neg, pos, predicate);
André Anjos's avatar
André Anjos committed
188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248
}

blitz::Array<double,2> bob::measure::roc(const blitz::Array<double,1>& negatives,
 const blitz::Array<double,1>& positives, size_t points) {
  double min = std::min(blitz::min(negatives), blitz::min(positives));
  double max = std::max(blitz::max(negatives), blitz::max(positives));
  double step = (max-min)/((double)points-1.0);
  blitz::Array<double,2> retval(2, points);
  for (int i=0; i<(int)points; ++i) {
    std::pair<double, double> ratios =
      bob::measure::farfrr(negatives, positives, min + i*step);
    // preserve X x Y ordering (FAR x FRR)
    retval(0,i) = ratios.first;
    retval(1,i) = ratios.second;
  }
  return retval;
}

blitz::Array<double,2> bob::measure::precision_recall_curve(const blitz::Array<double,1>& negatives,
 const blitz::Array<double,1>& positives, size_t points) {
  double min = std::min(blitz::min(negatives), blitz::min(positives));
  double max = std::max(blitz::max(negatives), blitz::max(positives));
  double step = (max-min)/((double)points-1.0);
  blitz::Array<double,2> retval(2, points);
  for (int i=0; i<(int)points; ++i) {
    std::pair<double, double> ratios =
      bob::measure::precision_recall(negatives, positives, min + i*step);
    retval(0,i) = ratios.first;
    retval(1,i) = ratios.second;
  }
  return retval;
}

/**
  * Structure for getting permutations when sorting an array
  */
struct ComparePairs
{
  ComparePairs(const blitz::Array<double,1> &v):
    m_v(v)
  {
  }

  bool operator()(size_t a, size_t b)
  {
    return m_v(a) < m_v(b);
  }

  blitz::Array<double,1> m_v;
};

/**
  * Sort an array and get the permutations (using stable_sort)
  */
void sortWithPermutation(const blitz::Array<double,1>& values, std::vector<size_t>& v)
{
  int N = values.extent(0);
  bob::core::array::assertSameDimensionLength(N, v.size());
  for(int i=0; i<N; ++i)
    v[i] = i;

249
  std::stable_sort(v.begin(), v.end(), ComparePairs(values));
André Anjos's avatar
André Anjos committed
250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364
}

blitz::Array<double,2> bob::measure::rocch(const blitz::Array<double,1>& negatives,
 const blitz::Array<double,1>& positives)
{
  // Number of positive and negative scores
  size_t Nt = positives.extent(0);
  size_t Nn = negatives.extent(0);
  size_t N = Nt + Nn;

  // Create a big array with all scores
  blitz::Array<double,1> scores(N);
  blitz::Range rall = blitz::Range::all();
  scores(blitz::Range(0,Nt-1)) = positives(rall);
  scores(blitz::Range(Nt,N-1)) = negatives(rall);

  // It is important here that scores that are the same (i.e. already in order) should NOT be swapped.
  // std::stable_sort has this property.
  std::vector<size_t> perturb(N);
  sortWithPermutation(scores, perturb);

  // Apply permutation
  blitz::Array<size_t,1> Pideal(N);
  for(size_t i=0; i<N; ++i)
    Pideal(i) = (perturb[i] < Nt ? 1 : 0);
  blitz::Array<double,1> Pideal_d = bob::core::array::cast<double>(Pideal);

  // Apply the PAVA algorithm
  blitz::Array<double,1> Popt(N);
  blitz::Array<size_t,1> width = bob::math::pavxWidth(Pideal_d, Popt);

  // Allocate output
  int nbins = width.extent(0);
  blitz::Array<double,2> retval(2,nbins+1); // FAR, FRR

  // Fill in output
  size_t left = 0;
  size_t fa = Nn;
  size_t miss = 0;
  for(int i=0; i<nbins; ++i)
  {
    retval(0,i) = fa / (double)Nn; // pfa
    retval(1,i) = miss / (double)Nt; // pmiss
    left += width(i);
    if(left >= 1)
      miss = blitz::sum(Pideal(blitz::Range(0,left-1)));
    else
      miss = 0;
    if(Pideal.extent(0)-1 >= (int)left)
      fa = N - left - blitz::sum(Pideal(blitz::Range(left,Pideal.extent(0)-1)));
    else
      fa = 0;
  }
  retval(0,nbins) = fa / (double)Nn; // pfa
  retval(1,nbins) = miss / (double)Nt; // pmiss

  return retval;
}

double bob::measure::rocch2eer(const blitz::Array<double,2>& pfa_pmiss)
{
  bob::core::array::assertSameDimensionLength(2, pfa_pmiss.extent(0));
  const int N = pfa_pmiss.extent(1);

  double eer = 0.;
  blitz::Array<double,2> XY(2,2);
  blitz::Array<double,1> one(2);
  one = 1.;
  blitz::Array<double,1> seg(2);
  double& XY00 = XY(0,0);
  double& XY01 = XY(0,1);
  double& XY10 = XY(1,0);
  double& XY11 = XY(1,1);

  double eerseg = 0.;
  for(int i=0; i<N-1; ++i)
  {
    // Define XY matrix
    XY00 = pfa_pmiss(0,i); // pfa
    XY10 = pfa_pmiss(0,i+1); // pfa
    XY01 = pfa_pmiss(1,i); // pmiss
    XY11 = pfa_pmiss(1,i+1); // pmiss
    // xx and yy should be sorted:
    assert(XY10 <= XY00 && XY01 <= XY11);

    // Commpute "dd"
    double abs_dd0 = std::fabs(XY00 - XY10);
    double abs_dd1 = std::fabs(XY01 - XY11);
    if(std::min(abs_dd0,abs_dd1) < std::numeric_limits<double>::epsilon())
      eerseg = 0.;
    else
    {
      // Find line coefficients seg s.t. XY.seg = 1,
      bob::math::linsolve_(XY, seg, one);
      // Candidate for the EER (to be compared to current value)
      eerseg = 1. / blitz::sum(seg);
    }

    eer = std::max(eer, eerseg);
  }

  return eer;
}

/**
 * This function computes the ROC coordinates for the given positive and
 * negative values at the given FAR positions.
 *
 * @param negatives Impostor scores
 * @param positives Client scores
 * @param far_list  The list of FAR values where the FRR should be calculated
 *
 * @return The ROC curve with the FAR in the first row and the FRR in the second.
 */
blitz::Array<double,2> bob::measure::roc_for_far(const blitz::Array<double,1>& negatives,
365
 const blitz::Array<double,1>& positives, const blitz::Array<double,1>& far_list, bool isSorted) {
André Anjos's avatar
André Anjos committed
366 367 368 369 370
  int n_points = far_list.extent(0);

  if (negatives.extent(0) == 0) throw std::runtime_error("The given set of negatives is empty.");
  if (positives.extent(0) == 0) throw std::runtime_error("The given set of positives is empty.");

371 372 373 374
  // sort negative and positive scores ascendantly
  blitz::Array<double,1> neg, pos;
  sort(negatives, neg, isSorted);
  sort(positives, pos, isSorted);
André Anjos's avatar
André Anjos committed
375 376 377 378 379 380 381

  // do some magic to compute the FRR list
  blitz::Array<double,2> retval(2, n_points);

  // index into the FAR and FRR list
  int far_index = n_points-1;
  int pos_index = 0, neg_index = 0;
382
  int n_pos = pos.extent(0), n_neg = neg.extent(0);
André Anjos's avatar
André Anjos committed
383 384

  // iterators into the result lists
385
  auto pos_it = pos.begin(), neg_it = neg.begin();
André Anjos's avatar
André Anjos committed
386 387 388
  // do some fast magic to compute the FRR values ;-)
  do{
    // check whether the current positive value is less than the current negative one
389
    if (*pos_it < *neg_it){
André Anjos's avatar
André Anjos committed
390 391 392 393 394 395 396 397 398 399 400 401
      // increase the positive count
      ++pos_index;
      // go to the next positive value
      ++pos_it;
    }else{
      // increase the negative count
      ++neg_index;
      // go to the next negative value
      ++neg_it;
    }
    // check, if we have reached a new FAR limit,
    // i.e. if the relative number of negative similarities is greater than 1-FAR (which is the CRR)
402

403 404

    if (((double)neg_index / (double)n_neg > 1. - far_list(far_index)) &&
405
       !(bob::core::isClose  ((double)neg_index / (double)n_neg, 1. - far_list(far_index), 1e-9, 1e-9)))   {
André Anjos's avatar
André Anjos committed
406 407
      // copy the far value
      retval(0,far_index) = far_list(far_index);
408 409
      // calculate the FRR for the current FAR
      retval(1,far_index) = (double)pos_index / (double)n_pos;
André Anjos's avatar
André Anjos committed
410 411 412 413 414
      // go to the next FAR value
      --far_index;
    }

  // do this, as long as there are elements in both lists left and not all FRR elements where calculated yet
415
  } while (pos_it != pos.end() && neg_it != neg.end() && far_index >= 0);
André Anjos's avatar
André Anjos committed
416

417
  // check if all FRR values have been set
André Anjos's avatar
André Anjos committed
418 419
  if (far_index >= 0){
    // walk to the end of both lists; at least one of both lists should already have reached its limit.
420 421
    while (pos_it++ != pos.end()) ++pos_index;
    while (neg_it++ != neg.end()) ++neg_index;
André Anjos's avatar
André Anjos committed
422 423 424 425 426 427
    // fill in the remaining elements of the CAR list
    do {
      // copy the FAR value
      retval(0,far_index) = far_list(far_index);
      // check if the criterion is fulfilled (should be, as long as the lowest far is not below 0)
      if ((double)neg_index / (double)n_neg > 1. - far_list(far_index)){
428 429
        // calculate the FRR for the current FAR
        retval(1,far_index) = (double)pos_index / (double)n_pos;
André Anjos's avatar
André Anjos committed
430
      } else {
431 432
        // set FRR to 1 (this should never happen, but might be due to numerical issues)
        retval(1,far_index) = 1.;
André Anjos's avatar
André Anjos committed
433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516
      }
    } while (far_index--);
  }

  return retval;
}


/**
 * The input to this function is a cumulative probability.  The output from
 * this function is the Normal deviate that corresponds to that probability.
 * For example:
 *
 *  INPUT | OUTPUT
 * -------+--------
 *  0.001 | -3.090
 *  0.01  | -2.326
 *  0.1   | -1.282
 *  0.5   |  0.0
 *  0.9   |  1.282
 *  0.99  |  2.326
 *  0.999 |  3.090
 */
static double _ppndf(double p) {
  //some constants we need for the calculation.
  //these come from the NIST implementation...
  static const double SPLIT = 0.42;
  static const double A0 = 2.5066282388;
  static const double A1 = -18.6150006252;
  static const double A2 = 41.3911977353;
  static const double A3 = -25.4410604963;
  static const double B1 = -8.4735109309;
  static const double B2 = 23.0833674374;
  static const double B3 = -21.0622410182;
  static const double B4 = 3.1308290983;
  static const double C0 = -2.7871893113;
  static const double C1 = -2.2979647913;
  static const double C2 = 4.8501412713;
  static const double C3 = 2.3212127685;
  static const double D1 = 3.5438892476;
  static const double D2 = 1.6370678189;
  static const double eps = 2.2204e-16;

  double retval;

  if (p >= 1.0) p = 1 - eps;
  if (p <= 0.0) p = eps;

  double q = p - 0.5;

  if (std::abs(q) <= SPLIT) {
    double r = q * q;
    retval = q * (((A3 * r + A2) * r + A1) * r + A0) /
      ((((B4 * r + B3) * r + B2) * r + B1) * r + 1.0);
  }
  else {
    //r = sqrt (log (0.5 - abs(q)));
    double r = (q > 0.0  ?  1.0 - p : p);
    if (r <= 0.0) throw std::runtime_error("measure::ppndf(): r <= 0.0!");
    r = sqrt ((-1.0) * log (r));
    retval = (((C3 * r + C2) * r + C1) * r + C0) / ((D2 * r + D1) * r + 1.0);
    if (q < 0) retval *= -1.0;
  }

  return retval;
}

namespace blitz {
  BZ_DECLARE_FUNCTION(_ppndf) ///< A blitz::Array binding
}

double bob::measure::ppndf (double value) { return _ppndf(value); }

blitz::Array<double,2> bob::measure::det(const blitz::Array<double,1>& negatives,
    const blitz::Array<double,1>& positives, size_t points) {
  blitz::Array<double,2> retval(2, points);
  retval = blitz::_ppndf(bob::measure::roc(negatives, positives, points));
  return retval;
}

blitz::Array<double,2> bob::measure::epc
(const blitz::Array<double,1>& dev_negatives,
 const blitz::Array<double,1>& dev_positives,
 const blitz::Array<double,1>& test_negatives,
517 518
 const blitz::Array<double,1>& test_positives, size_t points, bool isSorted,
 bool thresholds) {
519 520 521 522 523

  blitz::Array<double,1> dev_neg, dev_pos;
  sort(dev_negatives, dev_neg, isSorted);
  sort(dev_positives, dev_pos, isSorted);

André Anjos's avatar
André Anjos committed
524
  double step = 1.0/((double)points-1.0);
525 526
  auto retval_shape0 = (thresholds) ? 3 : 2;
  blitz::Array<double,2> retval(retval_shape0, points);
André Anjos's avatar
André Anjos committed
527 528 529
  for (int i=0; i<(int)points; ++i) {
    double alpha = (double)i*step;
    retval(0,i) = alpha;
530 531
    double threshold = bob::measure::minWeightedErrorRateThreshold(dev_neg,
        dev_pos, alpha, true);
André Anjos's avatar
André Anjos committed
532 533 534
    std::pair<double, double> ratios =
      bob::measure::farfrr(test_negatives, test_positives, threshold);
    retval(1,i) = (ratios.first + ratios.second) / 2;
535 536 537
    if (thresholds) {
      retval(2,i) = threshold;
    }
André Anjos's avatar
André Anjos committed
538 539 540
  }
  return retval;
}