error.cpp 19.2 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>
André Anjos's avatar
André Anjos committed
19

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

23
#include "error.h"
24

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

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
36
37
38
39
40
41
42
43
std::pair<double, double>
bob::measure::farfrr(const blitz::Array<double, 1> &negatives,
                     const blitz::Array<double, 1> &positives,
                     double threshold) {
  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
44
45
46
47
  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
48
49
  return std::make_pair(false_accepts / (double)total_negatives,
                        false_rejects / (double)total_positives);
André Anjos's avatar
André Anjos committed
50
51
}

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
52
53
54
55
56
57
58
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
59
60
61
62
  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
63
64
65
66
67
68
  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
69
70
}

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
71
72
73
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
74
75
76
77
  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
78
79
  if (weight <= 0)
    weight = 1;
André Anjos's avatar
André Anjos committed
80
81
  if (precision == 0 && recall == 0)
    return 0;
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
82
83
  return (1 + weight * weight) * precision * recall /
         (weight * weight * precision + recall);
André Anjos's avatar
André Anjos committed
84
85
}

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
86
double eer_predicate(double far, double frr) { return std::abs(far - frr); }
André Anjos's avatar
André Anjos committed
87

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
88
89
90
91
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;
Amir Mohammadi's avatar
Amir Mohammadi committed
92
93
  sort(negatives, neg, is_sorted);
  sort(positives, pos, is_sorted);
94
  return bob::measure::minimizingThreshold(neg, pos, eer_predicate);
André Anjos's avatar
André Anjos committed
95
96
}

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
97
98
double bob::measure::eerRocch(const blitz::Array<double, 1> &negatives,
                              const blitz::Array<double, 1> &positives) {
André Anjos's avatar
André Anjos committed
99
100
101
  return bob::measure::rocch2eer(bob::measure::rocch(negatives, positives));
}

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
102
103
104
double bob::measure::farThreshold(const blitz::Array<double, 1> &negatives,
                                  const blitz::Array<double, 1> &,
                                  double far_value, bool is_sorted) {
André Anjos's avatar
André Anjos committed
105
106
  // check the parameters are valid
  if (far_value < 0. || far_value > 1.) {
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
107
108
    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
109
110
111
112
    m % far_value;
    throw std::runtime_error(m.str());
  }
  if (negatives.size() < 2) {
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
113
114
    throw std::runtime_error(
        "the number of negative scores must be at least 2");
André Anjos's avatar
André Anjos committed
115
116
  }

117
  // sort the array, if necessary
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
118
  blitz::Array<double, 1> neg;
Amir Mohammadi's avatar
Amir Mohammadi committed
119
  sort(negatives, neg, is_sorted);
André Anjos's avatar
André Anjos committed
120
121

  // compute position of the threshold
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
122
  double crr = 1. - far_value; // (Correct Rejection Rate; = 1 - FAR)
123
  double crr_index = crr * neg.extent(0);
André Anjos's avatar
André Anjos committed
124
  // compute the index above the current CRR value
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
125
  int index = std::min((int)std::floor(crr_index), neg.extent(0) - 1);
André Anjos's avatar
André Anjos committed
126
127

  // correct index if we have multiple score values at the requested position
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
128
129
  while (index && neg(index) == neg(index - 1))
    --index;
André Anjos's avatar
André Anjos committed
130

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
131
132
  // we compute a correction term to assure that we are in the middle of two
  // cases
André Anjos's avatar
André Anjos committed
133
  double correction;
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
134
  if (index) {
André Anjos's avatar
André Anjos committed
135
    // assure that we are in the middle of two cases
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
136
    correction = 0.5 * (neg(index) - neg(index - 1));
André Anjos's avatar
André Anjos committed
137
138
  } else {
    // add an overall correction term
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
139
    correction = 0.5 * (neg(neg.extent(0) - 1) - neg(0)) / neg.extent(0);
André Anjos's avatar
André Anjos committed
140
141
  }

142
  return neg(index) - correction;
André Anjos's avatar
André Anjos committed
143
144
}

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
145
146
147
double bob::measure::frrThreshold(const blitz::Array<double, 1> &,
                                  const blitz::Array<double, 1> &positives,
                                  double frr_value, bool is_sorted) {
André Anjos's avatar
André Anjos committed
148
149
150

  // check the parameters are valid
  if (frr_value < 0. || frr_value > 1.) {
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
151
152
    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
153
154
155
156
    m % frr_value;
    throw std::runtime_error(m.str());
  }
  if (positives.size() < 2) {
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
157
158
    throw std::runtime_error(
        "the number of positive scores must be at least 2");
André Anjos's avatar
André Anjos committed
159
160
  }

161
  // sort positive scores descendantly, if necessary
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
162
  blitz::Array<double, 1> pos;
Amir Mohammadi's avatar
Amir Mohammadi committed
163
  sort(positives, pos, is_sorted);
André Anjos's avatar
André Anjos committed
164
165

  // compute position of the threshold
166
167
  double frr_index = frr_value * pos.extent(0);
  // compute the index above the current CAR value
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
168
  int index = std::min((int)std::ceil(frr_index), pos.extent(0) - 1);
André Anjos's avatar
André Anjos committed
169
170

  // correct index if we have multiple score values at the requested position
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
171
172
  while (index < pos.extent(0) - 1 && pos(index) == pos(index + 1))
    ++index;
André Anjos's avatar
André Anjos committed
173

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
174
175
  // we compute a correction term to assure that we are in the middle of two
  // cases
André Anjos's avatar
André Anjos committed
176
  double correction;
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
177
  if (index < pos.extent(0) - 1) {
André Anjos's avatar
André Anjos committed
178
    // assure that we are in the middle of two cases
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
179
    correction = 0.5 * (pos(index + 1) - pos(index));
André Anjos's avatar
André Anjos committed
180
181
  } else {
    // add an overall correction term
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
182
    correction = 0.5 * (pos(pos.extent(0) - 1) - pos(0)) / pos.extent(0);
André Anjos's avatar
André Anjos committed
183
184
  }

185
  return pos(index) + correction;
André Anjos's avatar
André Anjos committed
186
187
188
189
190
191
192
193
194
}

/**
 * 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
195
196
197
198
199
200
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
201
202
  }

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
203
204
  inline double operator()(double far, double frr) const {
    return (m_weight * far) + ((1.0 - m_weight) * frr);
André Anjos's avatar
André Anjos committed
205
206
207
  }
};

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
208
209
210
211
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;
Amir Mohammadi's avatar
Amir Mohammadi committed
212
213
  sort(negatives, neg, is_sorted);
  sort(positives, pos, is_sorted);
André Anjos's avatar
André Anjos committed
214
  weighted_error predicate(cost);
215
  return bob::measure::minimizingThreshold(neg, pos, predicate);
André Anjos's avatar
André Anjos committed
216
217
}

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
218
219
220
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
221
222
  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
223
224
225
  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
226
    std::pair<double, double> ratios =
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
227
        bob::measure::farfrr(negatives, positives, min + i * step);
André Anjos's avatar
André Anjos committed
228
    // preserve X x Y ordering (FAR x FRR)
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
229
230
    retval(0, i) = ratios.first;
    retval(1, i) = ratios.second;
André Anjos's avatar
André Anjos committed
231
232
233
234
  }
  return retval;
}

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
235
236
237
238
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
239
240
  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
241
242
243
  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
244
    std::pair<double, double> ratios =
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
245
246
247
        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
248
249
250
251
252
253
254
  }
  return retval;
}

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

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
258
  bool operator()(size_t a, size_t b) { return m_v(a) < m_v(b); }
André Anjos's avatar
André Anjos committed
259

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
260
  blitz::Array<double, 1> m_v;
André Anjos's avatar
André Anjos committed
261
262
263
264
265
};

/**
  * Sort an array and get the permutations (using stable_sort)
  */
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
266
267
void sortWithPermutation(const blitz::Array<double, 1> &values,
                         std::vector<size_t> &v) {
André Anjos's avatar
André Anjos committed
268
269
  int N = values.extent(0);
  bob::core::array::assertSameDimensionLength(N, v.size());
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
270
  for (int i = 0; i < N; ++i)
André Anjos's avatar
André Anjos committed
271
272
    v[i] = i;

273
  std::stable_sort(v.begin(), v.end(), ComparePairs(values));
André Anjos's avatar
André Anjos committed
274
275
}

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
276
277
278
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
279
280
281
282
283
284
  // 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
285
  blitz::Array<double, 1> scores(N);
André Anjos's avatar
André Anjos committed
286
  blitz::Range rall = blitz::Range::all();
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
287
288
  scores(blitz::Range(0, Nt - 1)) = positives(rall);
  scores(blitz::Range(Nt, N - 1)) = negatives(rall);
André Anjos's avatar
André Anjos committed
289

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
290
291
  // 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
292
293
294
295
296
  // 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
297
298
  blitz::Array<size_t, 1> Pideal(N);
  for (size_t i = 0; i < N; ++i)
André Anjos's avatar
André Anjos committed
299
    Pideal(i) = (perturb[i] < Nt ? 1 : 0);
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
300
  blitz::Array<double, 1> Pideal_d = bob::core::array::cast<double>(Pideal);
André Anjos's avatar
André Anjos committed
301
302

  // Apply the PAVA algorithm
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
303
304
  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
305
306
307

  // Allocate output
  int nbins = width.extent(0);
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
308
  blitz::Array<double, 2> retval(2, nbins + 1); // FAR, FRR
André Anjos's avatar
André Anjos committed
309
310
311
312
313

  // Fill in output
  size_t left = 0;
  size_t fa = Nn;
  size_t miss = 0;
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
314
315
316
  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
317
    left += width(i);
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
318
319
    if (left >= 1)
      miss = blitz::sum(Pideal(blitz::Range(0, left - 1)));
André Anjos's avatar
André Anjos committed
320
321
    else
      miss = 0;
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
322
323
324
    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
325
326
327
    else
      fa = 0;
  }
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
328
329
  retval(0, nbins) = fa / (double)Nn;   // pfa
  retval(1, nbins) = miss / (double)Nt; // pmiss
André Anjos's avatar
André Anjos committed
330
331
332
333

  return retval;
}

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
334
double bob::measure::rocch2eer(const blitz::Array<double, 2> &pfa_pmiss) {
André Anjos's avatar
André Anjos committed
335
336
337
338
  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
339
340
  blitz::Array<double, 2> XY(2, 2);
  blitz::Array<double, 1> one(2);
André Anjos's avatar
André Anjos committed
341
  one = 1.;
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
342
343
344
345
346
  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
347
348

  double eerseg = 0.;
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
349
  for (int i = 0; i < N - 1; ++i) {
André Anjos's avatar
André Anjos committed
350
    // Define XY matrix
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
351
352
353
354
    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
355
356
357
358
359
360
    // 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
361
    if (std::min(abs_dd0, abs_dd1) < std::numeric_limits<double>::epsilon())
André Anjos's avatar
André Anjos committed
362
      eerseg = 0.;
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
363
    else {
André Anjos's avatar
André Anjos committed
364
      // Find line coefficients seg s.t. XY.seg = 1,
365
      bob::math::linsolve(XY, one, seg);
André Anjos's avatar
André Anjos committed
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
      // 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
384
385
 * @return The ROC curve with the FAR in the first row and the FRR in the
 * second.
André Anjos's avatar
André Anjos committed
386
 */
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
387
388
389
390
391
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
392
393
  int n_points = far_list.extent(0);

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
394
395
396
397
  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
398

399
  // sort negative and positive scores ascendantly
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
400
  blitz::Array<double, 1> neg, pos;
Amir Mohammadi's avatar
Amir Mohammadi committed
401
402
  sort(negatives, neg, is_sorted);
  sort(positives, pos, is_sorted);
André Anjos's avatar
André Anjos committed
403
404

  // do some magic to compute the FRR list
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
405
  blitz::Array<double, 2> retval(2, n_points);
André Anjos's avatar
André Anjos committed
406
407

  // index into the FAR and FRR list
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
408
  int far_index = n_points - 1;
André Anjos's avatar
André Anjos committed
409
  int pos_index = 0, neg_index = 0;
410
  int n_pos = pos.extent(0), n_neg = neg.extent(0);
André Anjos's avatar
André Anjos committed
411
412

  // iterators into the result lists
413
  auto pos_it = pos.begin(), neg_it = neg.begin();
André Anjos's avatar
André Anjos committed
414
  // do some fast magic to compute the FRR values ;-)
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
415
416
417
418
  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
419
420
421
422
      // increase the positive count
      ++pos_index;
      // go to the next positive value
      ++pos_it;
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
423
    } else {
André Anjos's avatar
André Anjos committed
424
425
426
427
428
429
      // 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
430
431
    // i.e. if the relative number of negative similarities is greater than
    // 1-FAR (which is the CRR)
432
433

    if (((double)neg_index / (double)n_neg > 1. - far_list(far_index)) &&
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
434
435
        !(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
436
      // copy the far value
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
437
      retval(0, far_index) = far_list(far_index);
438
      // calculate the FRR for the current FAR
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
439
      retval(1, far_index) = (double)pos_index / (double)n_pos;
André Anjos's avatar
André Anjos committed
440
441
442
443
      // go to the next FAR value
      --far_index;
    }

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

448
  // check if all FRR values have been set
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
449
450
451
452
453
454
455
  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
456
457
458
    // fill in the remaining elements of the CAR list
    do {
      // copy the FAR value
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
459
460
461
462
      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)) {
463
        // calculate the FRR for the current FAR
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
464
        retval(1, far_index) = (double)pos_index / (double)n_pos;
André Anjos's avatar
André Anjos committed
465
      } else {
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
466
467
468
        // 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
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
      }
    } 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
492
493
  // some constants we need for the calculation.
  // these come from the NIST implementation...
André Anjos's avatar
André Anjos committed
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
  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
513
514
515
516
  if (p >= 1.0)
    p = 1 - eps;
  if (p <= 0.0)
    p = eps;
André Anjos's avatar
André Anjos committed
517
518
519
520
521
522

  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
523
524
525
526
527
528
529
             ((((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
530
    retval = (((C3 * r + C2) * r + C1) * r + C0) / ((D2 * r + D1) * r + 1.0);
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
531
532
    if (q < 0)
      retval *= -1.0;
André Anjos's avatar
André Anjos committed
533
534
535
536
537
538
  }

  return retval;
}

namespace blitz {
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
539
BZ_DECLARE_FUNCTION(_ppndf) ///< A blitz::Array binding
André Anjos's avatar
André Anjos committed
540
541
}

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
542
double bob::measure::ppndf(double value) { return _ppndf(value); }
André Anjos's avatar
André Anjos committed
543

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
544
545
546
547
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
548
549
550
551
  retval = blitz::_ppndf(bob::measure::roc(negatives, positives, points));
  return retval;
}

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
552
553
554
555
556
557
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) {
558

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
559
  blitz::Array<double, 1> dev_neg, dev_pos;
Amir Mohammadi's avatar
Amir Mohammadi committed
560
561
  sort(dev_negatives, dev_neg, is_sorted);
  sort(dev_positives, dev_pos, is_sorted);
562

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
563
  double step = 1.0 / ((double)points - 1.0);
564
  auto retval_shape0 = (thresholds) ? 3 : 2;
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
565
566
567
568
569
570
  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
571
    std::pair<double, double> ratios =
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
572
573
        bob::measure::farfrr(test_negatives, test_positives, threshold);
    retval(1, i) = (ratios.first + ratios.second) / 2;
574
    if (thresholds) {
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
575
      retval(2, i) = threshold;
576
    }
André Anjos's avatar
André Anjos committed
577
578
579
  }
  return retval;
}