error.cpp 19.6 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
lint  
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
lint  
Amir MOHAMMADI committed
17 18
#include <bob.core/assert.h>
#include <bob.core/cast.h>
19
#include <bob.core/logging.h>
André Anjos's avatar
André Anjos committed
20

21
#include <bob.math/linsolve.h>
Amir MOHAMMADI's avatar
lint  
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
lint  
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
lint  
Amir MOHAMMADI committed
32
    bob::core::array::ccopy(a, b);
33
    bob::core::array::sort<T>(b);
34 35 36
  }
}

Amir MOHAMMADI's avatar
lint  
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 FAR or FRR with threshold NaN.\n";
43 44
    return std::make_pair(1.,1.);
  }
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
45 46 47 48
  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
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
lint  
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
lint  
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
lint  
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
lint  
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
lint  
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
lint  
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
lint  
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
lint  
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
lint  
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
lint  
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
lint  
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
lint  
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
lint  
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
  // sort the array, if necessary
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
123
  blitz::Array<double, 1> neg;
124
  sort(negatives, neg, is_sorted);
125
  
126 127 128 129
  // far == 1 is a corner case
  if (far_value >= 1 - 1e-12)
    return neg(0) - 1e-12;

130 131 132
  // Move towards the beginning of array changing the threshold until we pass
  // the desired FAR value. Start with a threshold that corresponds to FAR ==
  // 0.
133 134 135 136
  int last_pivot = neg.extent(0);
  int pivot = last_pivot / 2; 

  double threshold = neg(last_pivot-1) + 1e-12;
137
  double future_far;
138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
  while (pivot < neg.extent(0)) {
    future_far = blitz::count(neg >= neg(pivot)) / (double)neg.extent(0);
    if (future_far <= far_value){
      // Moving towards the beggining of the list
      last_pivot = pivot;
      pivot = pivot / 2;      
    }
    else{
      //Checking if we are in the boundary
      if(pivot==neg.extent(0)-1){
          threshold = neg(pivot) + 1e-12;
          break;
      }      
      future_far = blitz::count(neg >= neg(pivot+1)) / (double)neg.extent(0);

      // If we are in the boundary we are done....
      if (future_far > far_value)
        //Let's keep searching towards the end of the list
        pivot = (last_pivot + pivot) / 2;
      else{
        // I found you....
        threshold = neg(pivot+1);
        break;
      }
    }
  }  
164
  return threshold;
André Anjos's avatar
André Anjos committed
165 166
}

167
double bob::measure::frrThreshold(const blitz::Array<double, 1> &negatives,
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
168 169
                                  const blitz::Array<double, 1> &positives,
                                  double frr_value, bool is_sorted) {
André Anjos's avatar
André Anjos committed
170 171 172

  // check the parameters are valid
  if (frr_value < 0. || frr_value > 1.) {
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
173 174
    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
175 176 177 178
    m % frr_value;
    throw std::runtime_error(m.str());
  }
  if (positives.size() < 2) {
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
179 180
    throw std::runtime_error(
        "the number of positive scores must be at least 2");
André Anjos's avatar
André Anjos committed
181 182
  }

183
  // sort positive scores descendantly, if necessary
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
184
  blitz::Array<double, 1> pos;
185
  sort(positives, pos, is_sorted);
André Anjos's avatar
André Anjos committed
186

187 188 189 190
  // frr == 1 is a corner case
  if (frr_value >= 1 - 1e-12)
    return pos(pos.extent(0)-1) + 1e-12;

191 192 193 194 195
  // Move towards the end of array changing the threshold until we pass
  // the desired FRR value. Start with a threshold that corresponds to FRR ==
  // 0.
  int index = 0;
  double threshold = pos(index) - 1e-12;
196
  double future_frr;
197 198
  while (index < pos.extent(0)) {
    future_frr = blitz::count(pos < pos(index)) / (double)pos.extent(0);
199 200
    if (future_frr > frr_value)
      break;
201 202
    threshold = pos(index);
    ++index;
André Anjos's avatar
André Anjos committed
203
  }
204
  return threshold;
André Anjos's avatar
André Anjos committed
205 206 207 208 209 210 211 212 213
}

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

  double m_weight; ///< The weighting factor

Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
214 215 216 217 218 219
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
220 221
  }

Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
222 223
  inline double operator()(double far, double frr) const {
    return (m_weight * far) + ((1.0 - m_weight) * frr);
André Anjos's avatar
André Anjos committed
224 225 226
  }
};

Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
227 228 229 230
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;
231 232
  sort(negatives, neg, is_sorted);
  sort(positives, pos, is_sorted);
André Anjos's avatar
André Anjos committed
233
  weighted_error predicate(cost);
234
  return bob::measure::minimizingThreshold(neg, pos, predicate);
André Anjos's avatar
André Anjos committed
235 236
}

Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
237 238 239
blitz::Array<double, 2>
bob::measure::roc(const blitz::Array<double, 1> &negatives,
                  const blitz::Array<double, 1> &positives, size_t points) {
André Anjos's avatar
André Anjos committed
240 241
  double min = std::min(blitz::min(negatives), blitz::min(positives));
  double max = std::max(blitz::max(negatives), blitz::max(positives));
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
242 243 244
  double step = (max - min) / ((double)points - 1.0);
  blitz::Array<double, 2> retval(2, points);
  for (int i = 0; i < (int)points; ++i) {
André Anjos's avatar
André Anjos committed
245
    std::pair<double, double> ratios =
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
246
        bob::measure::farfrr(negatives, positives, min + i * step);
André Anjos's avatar
André Anjos committed
247
    // preserve X x Y ordering (FAR x FRR)
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
248 249
    retval(0, i) = ratios.first;
    retval(1, i) = ratios.second;
André Anjos's avatar
André Anjos committed
250 251 252 253
  }
  return retval;
}

Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
254 255 256 257
blitz::Array<double, 2>
bob::measure::precision_recall_curve(const blitz::Array<double, 1> &negatives,
                                     const blitz::Array<double, 1> &positives,
                                     size_t points) {
André Anjos's avatar
André Anjos committed
258 259
  double min = std::min(blitz::min(negatives), blitz::min(positives));
  double max = std::max(blitz::max(negatives), blitz::max(positives));
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
260 261 262
  double step = (max - min) / ((double)points - 1.0);
  blitz::Array<double, 2> retval(2, points);
  for (int i = 0; i < (int)points; ++i) {
André Anjos's avatar
André Anjos committed
263
    std::pair<double, double> ratios =
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
264 265 266
        bob::measure::precision_recall(negatives, positives, min + i * step);
    retval(0, i) = ratios.first;
    retval(1, i) = ratios.second;
André Anjos's avatar
André Anjos committed
267 268 269 270 271 272 273
  }
  return retval;
}

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

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

Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
279
  blitz::Array<double, 1> m_v;
André Anjos's avatar
André Anjos committed
280 281 282 283 284
};

/**
  * Sort an array and get the permutations (using stable_sort)
  */
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
285 286
void sortWithPermutation(const blitz::Array<double, 1> &values,
                         std::vector<size_t> &v) {
André Anjos's avatar
André Anjos committed
287 288
  int N = values.extent(0);
  bob::core::array::assertSameDimensionLength(N, v.size());
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
289
  for (int i = 0; i < N; ++i)
André Anjos's avatar
André Anjos committed
290 291
    v[i] = i;

292
  std::stable_sort(v.begin(), v.end(), ComparePairs(values));
André Anjos's avatar
André Anjos committed
293 294
}

Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
295 296 297
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
298 299 300 301 302 303
  // 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
lint  
Amir MOHAMMADI committed
304
  blitz::Array<double, 1> scores(N);
André Anjos's avatar
André Anjos committed
305
  blitz::Range rall = blitz::Range::all();
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
306 307
  scores(blitz::Range(0, Nt - 1)) = positives(rall);
  scores(blitz::Range(Nt, N - 1)) = negatives(rall);
André Anjos's avatar
André Anjos committed
308

Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
309 310
  // 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
311 312 313 314 315
  // std::stable_sort has this property.
  std::vector<size_t> perturb(N);
  sortWithPermutation(scores, perturb);

  // Apply permutation
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
316 317
  blitz::Array<size_t, 1> Pideal(N);
  for (size_t i = 0; i < N; ++i)
André Anjos's avatar
André Anjos committed
318
    Pideal(i) = (perturb[i] < Nt ? 1 : 0);
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
319
  blitz::Array<double, 1> Pideal_d = bob::core::array::cast<double>(Pideal);
André Anjos's avatar
André Anjos committed
320 321

  // Apply the PAVA algorithm
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
322 323
  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
324 325 326

  // Allocate output
  int nbins = width.extent(0);
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
327
  blitz::Array<double, 2> retval(2, nbins + 1); // FAR, FRR
André Anjos's avatar
André Anjos committed
328 329 330 331 332

  // Fill in output
  size_t left = 0;
  size_t fa = Nn;
  size_t miss = 0;
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
333 334 335
  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
336
    left += width(i);
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
337 338
    if (left >= 1)
      miss = blitz::sum(Pideal(blitz::Range(0, left - 1)));
André Anjos's avatar
André Anjos committed
339 340
    else
      miss = 0;
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
341 342 343
    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
344 345 346
    else
      fa = 0;
  }
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
347 348
  retval(0, nbins) = fa / (double)Nn;   // pfa
  retval(1, nbins) = miss / (double)Nt; // pmiss
André Anjos's avatar
André Anjos committed
349 350 351 352

  return retval;
}

Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
353
double bob::measure::rocch2eer(const blitz::Array<double, 2> &pfa_pmiss) {
André Anjos's avatar
André Anjos committed
354 355 356 357
  bob::core::array::assertSameDimensionLength(2, pfa_pmiss.extent(0));
  const int N = pfa_pmiss.extent(1);

  double eer = 0.;
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
358 359
  blitz::Array<double, 2> XY(2, 2);
  blitz::Array<double, 1> one(2);
André Anjos's avatar
André Anjos committed
360
  one = 1.;
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
361 362 363 364 365
  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
366 367

  double eerseg = 0.;
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
368
  for (int i = 0; i < N - 1; ++i) {
André Anjos's avatar
André Anjos committed
369
    // Define XY matrix
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
370 371 372 373
    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
374 375 376 377 378 379
    // 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
lint  
Amir MOHAMMADI committed
380
    if (std::min(abs_dd0, abs_dd1) < std::numeric_limits<double>::epsilon())
André Anjos's avatar
André Anjos committed
381
      eerseg = 0.;
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
382
    else {
André Anjos's avatar
André Anjos committed
383
      // Find line coefficients seg s.t. XY.seg = 1,
384
      bob::math::linsolve_(XY, one, seg);
André Anjos's avatar
André Anjos committed
385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402
      // 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
 *
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
403 404
 * @return The ROC curve with the FAR in the first row and the FRR in the
 * second.
André Anjos's avatar
André Anjos committed
405
 */
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
406 407 408 409 410
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
411 412
  int n_points = far_list.extent(0);

Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
413 414 415 416
  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
417

418
  // sort negative and positive scores ascendantly
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
419
  blitz::Array<double, 1> neg, pos;
420 421
  sort(negatives, neg, is_sorted);
  sort(positives, pos, is_sorted);
André Anjos's avatar
André Anjos committed
422 423

  // do some magic to compute the FRR list
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
424
  blitz::Array<double, 2> retval(2, n_points);
André Anjos's avatar
André Anjos committed
425 426

  // index into the FAR and FRR list
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
427
  int far_index = n_points - 1;
André Anjos's avatar
André Anjos committed
428
  int pos_index = 0, neg_index = 0;
429
  int n_pos = pos.extent(0), n_neg = neg.extent(0);
André Anjos's avatar
André Anjos committed
430 431

  // iterators into the result lists
432
  auto pos_it = pos.begin(), neg_it = neg.begin();
André Anjos's avatar
André Anjos committed
433
  // do some fast magic to compute the FRR values ;-)
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
434 435 436 437
  do {
    // check whether the current positive value is less than the current
    // negative one
    if (*pos_it < *neg_it) {
André Anjos's avatar
André Anjos committed
438 439 440 441
      // increase the positive count
      ++pos_index;
      // go to the next positive value
      ++pos_it;
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
442
    } else {
André Anjos's avatar
André Anjos committed
443 444 445 446 447 448
      // increase the negative count
      ++neg_index;
      // go to the next negative value
      ++neg_it;
    }
    // check, if we have reached a new FAR limit,
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
449 450
    // i.e. if the relative number of negative similarities is greater than
    // 1-FAR (which is the CRR)
451 452

    if (((double)neg_index / (double)n_neg > 1. - far_list(far_index)) &&
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
453 454
        !(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
455
      // copy the far value
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
456
      retval(0, far_index) = far_list(far_index);
457
      // calculate the FRR for the current FAR
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
458
      retval(1, far_index) = (double)pos_index / (double)n_pos;
André Anjos's avatar
André Anjos committed
459 460 461 462
      // go to the next FAR value
      --far_index;
    }

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

467
  // check if all FRR values have been set
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
468 469 470 471 472 473 474
  if (far_index >= 0) {
    // walk to the end of both lists; at least one of both lists should already
    // have reached its limit.
    while (pos_it++ != pos.end())
      ++pos_index;
    while (neg_it++ != neg.end())
      ++neg_index;
André Anjos's avatar
André Anjos committed
475 476 477
    // fill in the remaining elements of the CAR list
    do {
      // copy the FAR value
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
478 479 480 481
      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)) {
482
        // calculate the FRR for the current FAR
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
483
        retval(1, far_index) = (double)pos_index / (double)n_pos;
André Anjos's avatar
André Anjos committed
484
      } else {
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
485 486 487
        // 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
488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510
      }
    } 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) {
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
511 512
  // some constants we need for the calculation.
  // these come from the NIST implementation...
André Anjos's avatar
André Anjos committed
513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531
  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
lint  
Amir MOHAMMADI committed
532 533 534 535
  if (p >= 1.0)
    p = 1 - eps;
  if (p <= 0.0)
    p = eps;
André Anjos's avatar
André Anjos committed
536 537 538 539 540 541

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

  return retval;
}

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

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

Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
563 564 565 566
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);
André Anjos's avatar
André Anjos committed
567 568 569 570
  retval = blitz::_ppndf(bob::measure::roc(negatives, positives, points));
  return retval;
}

Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
571 572 573 574 575 576
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) {
577

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

Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
582
  double step = 1.0 / ((double)points - 1.0);
583
  auto retval_shape0 = (thresholds) ? 3 : 2;
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
584 585 586 587 588 589
  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
590
    std::pair<double, double> ratios =
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
591 592
        bob::measure::farfrr(test_negatives, test_positives, threshold);
    retval(1, i) = (ratios.first + ratios.second) / 2;
593
    if (thresholds) {
Amir MOHAMMADI's avatar
lint  
Amir MOHAMMADI committed
594
      retval(2, i) = threshold;
595
    }
André Anjos's avatar
André Anjos committed
596 597 598
  }
  return retval;
}