error.cpp 20.1 KB
Newer Older
André Anjos's avatar
André Anjos committed
1 2 3 4 5 6 7 8 9 10 11
/**
 * @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 <algorithm>
#include <boost/format.hpp>
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
12 13
#include <limits>
#include <stdexcept>
André Anjos's avatar
André Anjos committed
14

15 16
#include <bob.core/array_copy.h>
#include <bob.core/array_sort.h>
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
17 18
#include <bob.core/assert.h>
#include <bob.core/cast.h>
19
#include <bob.core/logging.h>
20

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

24
#include "error.h"
25

26
template <typename T>
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
27 28 29
static void sort(const blitz::Array<T, 1> &a, blitz::Array<T, 1> &b,
                 bool is_sorted) {
  if (is_sorted) {
30 31
    b.reference(a);
  } else {
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
32
    bob::core::array::ccopy(a, b);
33
    bob::core::array::sort<T>(b);
34 35 36
  }
}

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
37 38 39 40
std::pair<double, double>
bob::measure::farfrr(const blitz::Array<double, 1> &negatives,
                     const blitz::Array<double, 1> &positives,
                     double threshold) {
41
  if (std::isnan(threshold)){
42
    bob::core::error << "Cannot compute FPR (FAR) or FNR (FRR) with threshold NaN.\n";
43 44
    return std::make_pair(1.,1.);
  }
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
45
  if (!negatives.size())
46
    throw std::runtime_error("Cannot compute FPR (FAR) when no negatives are given");
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
47
  if (!positives.size())
48
    throw std::runtime_error("Cannot compute FNR (FRR) when no positives are given");
André Anjos's avatar
André Anjos committed
49 50 51 52
  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);
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
53 54
  return std::make_pair(false_accepts / (double)total_negatives,
                        false_rejects / (double)total_positives);
André Anjos's avatar
André Anjos committed
55 56
}

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
57 58 59 60 61 62 63
std::pair<double, double>
bob::measure::precision_recall(const blitz::Array<double, 1> &negatives,
                               const blitz::Array<double, 1> &positives,
                               double threshold) {
  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
64 65 66 67
  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;
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
68 69 70 71 72 73
  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));
André Anjos's avatar
André Anjos committed
74 75
}

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
76 77 78
double bob::measure::f_score(const blitz::Array<double, 1> &negatives,
                             const blitz::Array<double, 1> &positives,
                             double threshold, double weight) {
André Anjos's avatar
André Anjos committed
79 80 81 82
  std::pair<double, double> ratios =
      bob::measure::precision_recall(negatives, positives, threshold);
  double precision = ratios.first;
  double recall = ratios.second;
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
83 84
  if (weight <= 0)
    weight = 1;
André Anjos's avatar
André Anjos committed
85 86
  if (precision == 0 && recall == 0)
    return 0;
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
87 88
  return (1 + weight * weight) * precision * recall /
         (weight * weight * precision + recall);
André Anjos's avatar
André Anjos committed
89 90
}

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
91
double eer_predicate(double far, double frr) { return std::abs(far - frr); }
André Anjos's avatar
André Anjos committed
92

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
93 94 95 96
double bob::measure::eerThreshold(const blitz::Array<double, 1> &negatives,
                                  const blitz::Array<double, 1> &positives,
                                  bool is_sorted) {
  blitz::Array<double, 1> neg, pos;
97 98
  sort(negatives, neg, is_sorted);
  sort(positives, pos, is_sorted);
99
  return bob::measure::minimizingThreshold(neg, pos, eer_predicate);
André Anjos's avatar
André Anjos committed
100 101
}

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
102 103
double bob::measure::eerRocch(const blitz::Array<double, 1> &negatives,
                              const blitz::Array<double, 1> &positives) {
André Anjos's avatar
André Anjos committed
104 105 106
  return bob::measure::rocch2eer(bob::measure::rocch(negatives, positives));
}

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
107
double bob::measure::farThreshold(const blitz::Array<double, 1> &negatives,
108
                                  const blitz::Array<double, 1> &positives,
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
109
                                  double far_value, bool is_sorted) {
André Anjos's avatar
André Anjos committed
110 111
  // check the parameters are valid
  if (far_value < 0. || far_value > 1.) {
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
112 113
    boost::format m("the argument for `far_value' cannot take the value %f - "
                    "the value must be in the interval [0.,1.]");
André Anjos's avatar
André Anjos committed
114 115 116 117
    m % far_value;
    throw std::runtime_error(m.str());
  }
  if (negatives.size() < 2) {
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
118 119
    throw std::runtime_error(
        "the number of negative scores must be at least 2");
André Anjos's avatar
André Anjos committed
120 121
  }

122 123 124
  // sort the negatives array, if necessary, and keep it in the scores variable
  blitz::Array<double, 1> scores;
  sort(negatives, scores, is_sorted);
André Anjos's avatar
André Anjos committed
125

126
  double epsilon = std::numeric_limits<double>::epsilon();
127
  // handle special case of far == 1 without any iterating
128
  if (far_value >= 1 - epsilon)
129
    return std::nextafter(scores(0), scores(0)-1);
130 131 132 133 134 135 136 137 138 139 140 141

  // Reverse negatives so the end is the start. This way the code below will be
  // very similar to the implementation in the frrThreshold function. The
  // implementations are not exactly the same though.
  scores.reverseSelf(0);
  // Move towards the end of array changing the threshold until we pass the
  // desired FAR value. Starting with a threshold that corresponds to FAR == 0.
  int total_count = scores.extent(0);
  int current_position = 0;
  // since the comparison is `if score >= threshold then accept as genuine`, we
  // can choose the largest score value + eps as the threshold so that we can
  // get for 0% FAR.
142
  double valid_threshold = std::nextafter(scores(current_position), scores(current_position)+1);
143
  double current_threshold;
144
  double future_far;
145 146 147 148 149 150 151 152
  while (current_position < total_count) {
    current_threshold = scores(current_position);
    // keep iterating if values are repeated
    while (current_position < total_count-1 && scores(current_position+1) == current_threshold)
      current_position++;
    // All the scores up to the current position and including the current
    // position will be accepted falsely.
    future_far = (double)(current_position+1) / (double)total_count;
153 154
    if (future_far > far_value)
      break;
155 156
    valid_threshold = current_threshold;
    current_position++;
André Anjos's avatar
André Anjos committed
157
  }
158
  return valid_threshold;
André Anjos's avatar
André Anjos committed
159 160
}

161
double bob::measure::frrThreshold(const blitz::Array<double, 1> &negatives,
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
162 163
                                  const blitz::Array<double, 1> &positives,
                                  double frr_value, bool is_sorted) {
André Anjos's avatar
André Anjos committed
164 165 166

  // check the parameters are valid
  if (frr_value < 0. || frr_value > 1.) {
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
167 168
    boost::format m("the argument for `frr_value' cannot take the value %f - "
                    "the value must be in the interval [0.,1.]");
André Anjos's avatar
André Anjos committed
169 170 171 172
    m % frr_value;
    throw std::runtime_error(m.str());
  }
  if (positives.size() < 2) {
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
173 174
    throw std::runtime_error(
        "the number of positive scores must be at least 2");
André Anjos's avatar
André Anjos committed
175 176
  }

177 178 179
  // sort the positives array, if necessary, and keep it in the scores variable
  blitz::Array<double, 1> scores;
  sort(positives, scores, is_sorted);
André Anjos's avatar
André Anjos committed
180

181
  double epsilon = std::numeric_limits<double>::epsilon();
182
  // handle special case of frr == 1 without any iterating
183
  if (frr_value >= 1 - epsilon)
184
    return std::nextafter(scores(scores.extent(0)-1), scores(scores.extent(0)-1)+1);
185 186 187 188 189 190 191 192 193

  // Move towards the end of array changing the threshold until we pass the
  // desired FRR value. Starting with a threshold that corresponds to FRR == 0.
  int total_count = scores.extent(0);
  int current_position = 0;
  // since the comparison is `if score >= threshold then accept as genuine`, we
  // can use the smallest positive score as the threshold for 0% FRR.
  double valid_threshold = scores(current_position);
  double current_threshold;
194
  double future_frr;
195 196 197 198 199 200 201 202
  while (current_position < total_count) {
    current_threshold = scores(current_position);
    // keep iterating if values are repeated
    while (current_position < total_count-1 && scores(current_position+1) == current_threshold)
      current_position++;
    // All the scores up to the current_position but not including
    // current_position will be rejected falsely.
    future_frr = (double)current_position / (double)total_count;
203 204
    if (future_frr > frr_value)
      break;
205 206
    valid_threshold = current_threshold;
    current_position++;
André Anjos's avatar
André Anjos committed
207
  }
208
  return valid_threshold;
André Anjos's avatar
André Anjos committed
209 210 211 212 213 214 215 216 217
}

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

  double m_weight; ///< The weighting factor

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
218 219 220 221 222 223
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;
André Anjos's avatar
André Anjos committed
224 225
  }

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
226 227
  inline double operator()(double far, double frr) const {
    return (m_weight * far) + ((1.0 - m_weight) * frr);
André Anjos's avatar
André Anjos committed
228 229 230
  }
};

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
231 232 233 234
double bob::measure::minWeightedErrorRateThreshold(
    const blitz::Array<double, 1> &negatives,
    const blitz::Array<double, 1> &positives, double cost, bool is_sorted) {
  blitz::Array<double, 1> neg, pos;
235 236
  sort(negatives, neg, is_sorted);
  sort(positives, pos, is_sorted);
André Anjos's avatar
André Anjos committed
237
  weighted_error predicate(cost);
238
  return bob::measure::minimizingThreshold(neg, pos, predicate);
André Anjos's avatar
André Anjos committed
239 240
}

241 242 243 244 245 246 247 248 249 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
blitz::Array<double, 1>
bob::measure::log_values(size_t points_, int min_power) {
  int points = (int)points_;
  blitz::Array<double, 1> retval(points);
  double counts_per_step = points / (-min_power) ;
  for (int i = 1-points; i <= 0; ++i) {
    retval(i+points-1) = std::pow(10., (double)i/counts_per_step);
  }
  return retval;
}


blitz::Array<double, 1>
bob::measure::meaningfulThresholds(
    const blitz::Array<double, 1> &negatives,
    const blitz::Array<double, 1> &positives,
    size_t points_, int min_far, bool is_sorted) {

  int points = (int)points_;
  int half_points = points / 2;
  blitz::Array<double, 1> thresholds(points);
  blitz::Array<double, 1> neg, pos;

  // sort negatives and positives
  sort(negatives, neg, is_sorted);
  sort(positives, pos, is_sorted);

  // Create an far_list and frr_list
  auto frr_list = bob::measure::log_values(half_points, min_far);
  auto far_list = bob::measure::log_values(points - half_points, min_far);

  // Compute thresholds for far_list and frr_list
  for (int i = 0; i < points; ++i) {
    if (i < half_points)
      thresholds(i) = bob::measure::frrThreshold(neg, pos, frr_list(i), true);
    else
      thresholds(i) = bob::measure::farThreshold(neg, pos, far_list(i-half_points), true);
  }

  // Sort the thresholds
  bob::core::array::sort(thresholds);

  return thresholds;
}


Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
287 288
blitz::Array<double, 2>
bob::measure::roc(const blitz::Array<double, 1> &negatives,
289 290 291 292 293 294 295 296 297 298 299 300 301
                  const blitz::Array<double, 1> &positives,
                  size_t points_, int min_far) {
  int points = (int)points_;
  blitz::Array<double, 2> retval(2, points);

  auto thresholds = bob::measure::meaningfulThresholds(
    negatives, positives, points_, min_far);

  // compute far and frr based on these thresholds
  for (int i = 0; i < points; ++i) {
    auto temp = bob::measure::farfrr(negatives, positives, thresholds(i));
    retval(0, i) = temp.first;
    retval(1, i) = temp.second;
André Anjos's avatar
André Anjos committed
302
  }
303
  return retval;
André Anjos's avatar
André Anjos committed
304 305
}

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
306 307 308 309 310
blitz::Array<double, 2>
bob::measure::precision_recall_curve(const blitz::Array<double, 1> &negatives,
                                     const blitz::Array<double, 1> &positives,
                                     size_t points) {
  blitz::Array<double, 2> retval(2, points);
311 312
  auto thresholds = bob::measure::meaningfulThresholds(
    negatives, positives, points);
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
313
  for (int i = 0; i < (int)points; ++i) {
314
    auto ratios = bob::measure::precision_recall(negatives, positives, thresholds(i));
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
315 316
    retval(0, i) = ratios.first;
    retval(1, i) = ratios.second;
André Anjos's avatar
André Anjos committed
317 318 319 320 321 322 323
  }
  return retval;
}

/**
  * Structure for getting permutations when sorting an array
  */
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
324 325
struct ComparePairs {
  ComparePairs(const blitz::Array<double, 1> &v) : m_v(v) {}
André Anjos's avatar
André Anjos committed
326

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
327
  bool operator()(size_t a, size_t b) { return m_v(a) < m_v(b); }
André Anjos's avatar
André Anjos committed
328

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
329
  blitz::Array<double, 1> m_v;
André Anjos's avatar
André Anjos committed
330 331 332 333 334
};

/**
  * Sort an array and get the permutations (using stable_sort)
  */
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
335 336
void sortWithPermutation(const blitz::Array<double, 1> &values,
                         std::vector<size_t> &v) {
André Anjos's avatar
André Anjos committed
337 338
  int N = values.extent(0);
  bob::core::array::assertSameDimensionLength(N, v.size());
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
339
  for (int i = 0; i < N; ++i)
André Anjos's avatar
André Anjos committed
340 341
    v[i] = i;

342
  std::stable_sort(v.begin(), v.end(), ComparePairs(values));
André Anjos's avatar
André Anjos committed
343 344
}

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
345 346 347
blitz::Array<double, 2>
bob::measure::rocch(const blitz::Array<double, 1> &negatives,
                    const blitz::Array<double, 1> &positives) {
André Anjos's avatar
André Anjos committed
348 349 350 351 352 353
  // 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
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
354
  blitz::Array<double, 1> scores(N);
André Anjos's avatar
André Anjos committed
355
  blitz::Range rall = blitz::Range::all();
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
356 357
  scores(blitz::Range(0, Nt - 1)) = positives(rall);
  scores(blitz::Range(Nt, N - 1)) = negatives(rall);
André Anjos's avatar
André Anjos committed
358

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
359 360
  // It is important here that scores that are the same (i.e. already in order)
  // should NOT be swapped.
André Anjos's avatar
André Anjos committed
361 362 363 364 365
  // std::stable_sort has this property.
  std::vector<size_t> perturb(N);
  sortWithPermutation(scores, perturb);

  // Apply permutation
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
366 367
  blitz::Array<size_t, 1> Pideal(N);
  for (size_t i = 0; i < N; ++i)
André Anjos's avatar
André Anjos committed
368
    Pideal(i) = (perturb[i] < Nt ? 1 : 0);
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
369
  blitz::Array<double, 1> Pideal_d = bob::core::array::cast<double>(Pideal);
André Anjos's avatar
André Anjos committed
370 371

  // Apply the PAVA algorithm
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
372 373
  blitz::Array<double, 1> Popt(N);
  blitz::Array<size_t, 1> width = bob::math::pavxWidth(Pideal_d, Popt);
André Anjos's avatar
André Anjos committed
374 375 376

  // Allocate output
  int nbins = width.extent(0);
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
377
  blitz::Array<double, 2> retval(2, nbins + 1); // FAR, FRR
André Anjos's avatar
André Anjos committed
378 379 380 381 382

  // Fill in output
  size_t left = 0;
  size_t fa = Nn;
  size_t miss = 0;
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
383 384 385
  for (int i = 0; i < nbins; ++i) {
    retval(0, i) = fa / (double)Nn;   // pfa
    retval(1, i) = miss / (double)Nt; // pmiss
André Anjos's avatar
André Anjos committed
386
    left += width(i);
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
387 388
    if (left >= 1)
      miss = blitz::sum(Pideal(blitz::Range(0, left - 1)));
André Anjos's avatar
André Anjos committed
389 390
    else
      miss = 0;
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
391 392 393
    if (Pideal.extent(0) - 1 >= (int)left)
      fa = N - left -
           blitz::sum(Pideal(blitz::Range(left, Pideal.extent(0) - 1)));
André Anjos's avatar
André Anjos committed
394 395 396
    else
      fa = 0;
  }
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
397 398
  retval(0, nbins) = fa / (double)Nn;   // pfa
  retval(1, nbins) = miss / (double)Nt; // pmiss
André Anjos's avatar
André Anjos committed
399 400 401 402

  return retval;
}

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
403
double bob::measure::rocch2eer(const blitz::Array<double, 2> &pfa_pmiss) {
André Anjos's avatar
André Anjos committed
404 405 406 407
  bob::core::array::assertSameDimensionLength(2, pfa_pmiss.extent(0));
  const int N = pfa_pmiss.extent(1);

  double eer = 0.;
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
408 409
  blitz::Array<double, 2> XY(2, 2);
  blitz::Array<double, 1> one(2);
André Anjos's avatar
André Anjos committed
410
  one = 1.;
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
411 412 413 414 415
  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);
André Anjos's avatar
André Anjos committed
416 417

  double eerseg = 0.;
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
418
  for (int i = 0; i < N - 1; ++i) {
André Anjos's avatar
André Anjos committed
419
    // Define XY matrix
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
420 421 422 423
    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
André Anjos's avatar
André Anjos committed
424 425 426 427 428 429
    // 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);
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
430
    if (std::min(abs_dd0, abs_dd1) < std::numeric_limits<double>::epsilon())
André Anjos's avatar
André Anjos committed
431
      eerseg = 0.;
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
432
    else {
André Anjos's avatar
André Anjos committed
433
      // Find line coefficients seg s.t. XY.seg = 1,
434
      bob::math::linsolve_(XY, one, seg);
André Anjos's avatar
André Anjos committed
435 436 437 438 439 440 441 442 443 444 445 446
      // 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
447
 * negative values at the given FPR (FAR) positions.
André Anjos's avatar
André Anjos committed
448 449 450
 *
 * @param negatives Impostor scores
 * @param positives Client scores
451
 * @param far_list  The list of FPR (FAR) values where the FNR (FRR) should be calculated
André Anjos's avatar
André Anjos committed
452
 *
453
 * @return The ROC curve with the FPR in the first row and the FNR in the
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
454
 * second.
André Anjos's avatar
André Anjos committed
455
 */
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
456 457 458 459 460
blitz::Array<double, 2>
bob::measure::roc_for_far(const blitz::Array<double, 1> &negatives,
                          const blitz::Array<double, 1> &positives,
                          const blitz::Array<double, 1> &far_list,
                          bool is_sorted) {
André Anjos's avatar
André Anjos committed
461 462
  int n_points = far_list.extent(0);

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
463 464 465 466
  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.");
André Anjos's avatar
André Anjos committed
467

468
  // sort negative and positive scores ascendantly
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
469
  blitz::Array<double, 1> neg, pos;
470 471
  sort(negatives, neg, is_sorted);
  sort(positives, pos, is_sorted);
André Anjos's avatar
André Anjos committed
472

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
473
  blitz::Array<double, 2> retval(2, n_points);
André Anjos's avatar
André Anjos committed
474

475
  // index into the FAR list
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
476
  int far_index = n_points - 1;
André Anjos's avatar
André Anjos committed
477

478 479 480 481 482 483 484 485 486 487
  // Get the threshold for the requested far values and calculate far and frr
  // values based on the threshold.
  while(far_index >= 0) {
    // calculate the threshold for the requested far
    auto threshold = bob::measure::farThreshold(neg, pos, far_list(far_index), true);
    // calculate the frr and re-calculate the far
    auto farfrr = bob::measure::farfrr(neg, pos, threshold);
    retval(0, far_index) = farfrr.first;
    retval(1, far_index) = farfrr.second;
    far_index--;
André Anjos's avatar
André Anjos committed
488 489 490 491 492
  }

  return retval;
}

493

André Anjos's avatar
André Anjos committed
494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509
/**
 * 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) {
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
510 511
  // some constants we need for the calculation.
  // these come from the NIST implementation...
André Anjos's avatar
André Anjos committed
512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530
  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;

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
531 532 533 534
  if (p >= 1.0)
    p = 1 - eps;
  if (p <= 0.0)
    p = eps;
André Anjos's avatar
André Anjos committed
535 536 537 538 539 540

  double q = p - 0.5;

  if (std::abs(q) <= SPLIT) {
    double r = q * q;
    retval = q * (((A3 * r + A2) * r + A1) * r + A0) /
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
541 542 543 544 545 546 547
             ((((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));
André Anjos's avatar
André Anjos committed
548
    retval = (((C3 * r + C2) * r + C1) * r + C0) / ((D2 * r + D1) * r + 1.0);
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
549 550
    if (q < 0)
      retval *= -1.0;
André Anjos's avatar
André Anjos committed
551 552 553 554 555 556
  }

  return retval;
}

namespace blitz {
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
557
BZ_DECLARE_FUNCTION(_ppndf) ///< A blitz::Array binding
André Anjos's avatar
André Anjos committed
558 559
}

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
560
double bob::measure::ppndf(double value) { return _ppndf(value); }
André Anjos's avatar
André Anjos committed
561

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
562 563
blitz::Array<double, 2>
bob::measure::det(const blitz::Array<double, 1> &negatives,
564
                  const blitz::Array<double, 1> &positives, size_t points, int min_far) {
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
565
  blitz::Array<double, 2> retval(2, points);
566
  retval = blitz::_ppndf(bob::measure::roc(negatives, positives, points, min_far));
André Anjos's avatar
André Anjos committed
567 568 569
  return retval;
}

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
570 571 572 573 574 575
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,
                  const blitz::Array<double, 1> &test_positives, size_t points,
                  bool is_sorted, bool thresholds) {
576

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
577
  blitz::Array<double, 1> dev_neg, dev_pos;
578 579
  sort(dev_negatives, dev_neg, is_sorted);
  sort(dev_positives, dev_pos, is_sorted);
580

Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
581
  double step = 1.0 / ((double)points - 1.0);
582
  auto retval_shape0 = (thresholds) ? 3 : 2;
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
583 584 585 586 587 588
  blitz::Array<double, 2> retval(retval_shape0, points);
  for (int i = 0; i < (int)points; ++i) {
    double alpha = (double)i * step;
    retval(0, i) = alpha;
    double threshold = bob::measure::minWeightedErrorRateThreshold(
        dev_neg, dev_pos, alpha, true);
André Anjos's avatar
André Anjos committed
589
    std::pair<double, double> ratios =
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
590 591
        bob::measure::farfrr(test_negatives, test_positives, threshold);
    retval(1, i) = (ratios.first + ratios.second) / 2;
592
    if (thresholds) {
Amir MOHAMMADI's avatar
Amir MOHAMMADI committed
593
      retval(2, i) = threshold;
594
    }
André Anjos's avatar
André Anjos committed
595 596 597
  }
  return retval;
}