error.cpp 17.9 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
519
520
521
522
 const blitz::Array<double,1>& test_positives, size_t points, bool isSorted) {

  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
523
524
525
526
527
  double step = 1.0/((double)points-1.0);
  blitz::Array<double,2> retval(2, points);
  for (int i=0; i<(int)points; ++i) {
    double alpha = (double)i*step;
    retval(0,i) = alpha;
528
529
    double threshold = bob::measure::minWeightedErrorRateThreshold(dev_neg,
        dev_pos, alpha, true);
André Anjos's avatar
André Anjos committed
530
531
532
533
534
535
    std::pair<double, double> ratios =
      bob::measure::farfrr(test_negatives, test_positives, threshold);
    retval(1,i) = (ratios.first + ratios.second) / 2;
  }
  return retval;
}