error.cpp 19.3 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)
Manuel Günther's avatar
Manuel Günther committed
123
  double crr_index = std::max(crr * neg.extent(0) - 1., 0.);
André Anjos's avatar
André Anjos committed
124
  // compute the index above the current CRR value
125
  int index = std::min((int)std::ceil(crr_index), neg.extent(0)-1);
André Anjos's avatar
André Anjos committed
126

127
128
  // increase the threshold when we have several negatives with the same score
  while (index < neg.extent(0)-1 && neg(index) == neg(index+1)) ++index;
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
129
    --index;
André Anjos's avatar
André Anjos committed
130

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
131
  // we compute a correction term to assure that we are in the middle of two
132
  if (index < neg.extent(0)-1){
André Anjos's avatar
André Anjos committed
133
    // assure that we are in the middle of two cases
134
135
    double correction = 0.5 * (neg(index+1) - neg(index));
    return neg(index) + correction;
André Anjos's avatar
André Anjos committed
136
  } else {
137
138
    // We cannot reach the desired threshold, as we have too many identical lowest scores, or the number of scores is too low
    return std::numeric_limits<double>::quiet_NaN();
André Anjos's avatar
André Anjos committed
139
140
141
  }
}

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
142
143
144
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
145
146
147

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

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

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

167
168
  // lower the threshold when several positives have the same score
  while (index && pos(index) == pos(index-1)) --index;
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
169
    ++index;
André Anjos's avatar
André Anjos committed
170

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
171
  // we compute a correction term to assure that we are in the middle of two
172
  if (index){
André Anjos's avatar
André Anjos committed
173
    // assure that we are in the middle of two cases
174
175
    double correction = 0.5 * (pos(index) - pos(index-1));
    return pos(index) - correction;
André Anjos's avatar
André Anjos committed
176
  } else {
177
178
    // We cannot reach the desired threshold, as we have too many identical highest scores
    return std::numeric_limits<double>::quiet_NaN();
André Anjos's avatar
André Anjos committed
179
180
181
182
183
184
185
186
187
188
  }
}

/**
 * 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
189
190
191
192
193
194
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
195
196
  }

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
197
198
  inline double operator()(double far, double frr) const {
    return (m_weight * far) + ((1.0 - m_weight) * frr);
André Anjos's avatar
André Anjos committed
199
200
201
  }
};

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
202
203
204
205
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
206
207
  sort(negatives, neg, is_sorted);
  sort(positives, pos, is_sorted);
André Anjos's avatar
André Anjos committed
208
  weighted_error predicate(cost);
209
  return bob::measure::minimizingThreshold(neg, pos, predicate);
André Anjos's avatar
André Anjos committed
210
211
}

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

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
229
230
231
232
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
233
234
  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
235
236
237
  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
238
    std::pair<double, double> ratios =
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
239
240
241
        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
242
243
244
245
246
247
248
  }
  return retval;
}

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

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

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
254
  blitz::Array<double, 1> m_v;
André Anjos's avatar
André Anjos committed
255
256
257
258
259
};

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

267
  std::stable_sort(v.begin(), v.end(), ComparePairs(values));
André Anjos's avatar
André Anjos committed
268
269
}

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
270
271
272
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
273
274
275
276
277
278
  // 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
279
  blitz::Array<double, 1> scores(N);
André Anjos's avatar
André Anjos committed
280
  blitz::Range rall = blitz::Range::all();
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
281
282
  scores(blitz::Range(0, Nt - 1)) = positives(rall);
  scores(blitz::Range(Nt, N - 1)) = negatives(rall);
André Anjos's avatar
André Anjos committed
283

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

  // Apply the PAVA algorithm
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
297
298
  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
299
300
301

  // Allocate output
  int nbins = width.extent(0);
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
302
  blitz::Array<double, 2> retval(2, nbins + 1); // FAR, FRR
André Anjos's avatar
André Anjos committed
303
304
305
306
307

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

  return retval;
}

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

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

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
388
389
390
391
  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
392

393
  // sort negative and positive scores ascendantly
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
394
  blitz::Array<double, 1> neg, pos;
Amir Mohammadi's avatar
Amir Mohammadi committed
395
396
  sort(negatives, neg, is_sorted);
  sort(positives, pos, is_sorted);
André Anjos's avatar
André Anjos committed
397
398

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

  // index into the FAR and FRR list
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
402
  int far_index = n_points - 1;
André Anjos's avatar
André Anjos committed
403
  int pos_index = 0, neg_index = 0;
404
  int n_pos = pos.extent(0), n_neg = neg.extent(0);
André Anjos's avatar
André Anjos committed
405
406

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

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

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

442
  // check if all FRR values have been set
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
443
444
445
446
447
448
449
  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
450
451
452
    // fill in the remaining elements of the CAR list
    do {
      // copy the FAR value
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
453
454
455
456
      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)) {
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
      } else {
Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
460
461
462
        // 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
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
      }
    } 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
486
487
  // some constants we need for the calculation.
  // these come from the NIST implementation...
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
  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
507
508
509
510
  if (p >= 1.0)
    p = 1 - eps;
  if (p <= 0.0)
    p = eps;
André Anjos's avatar
André Anjos committed
511
512
513
514
515
516

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

  return retval;
}

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

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

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
538
539
540
541
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
542
543
544
545
  retval = blitz::_ppndf(bob::measure::roc(negatives, positives, points));
  return retval;
}

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
546
547
548
549
550
551
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) {
552

Amir Mohammadi's avatar
lint    
Amir Mohammadi committed
553
  blitz::Array<double, 1> dev_neg, dev_pos;
Amir Mohammadi's avatar
Amir Mohammadi committed
554
555
  sort(dev_negatives, dev_neg, is_sorted);
  sort(dev_positives, dev_pos, is_sorted);
556

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