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>
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 42 43 44
  if (std::isnan(threshold)){
    bob::core::error << "Cannot compute FAR or FRR with threshold NaN";
    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 108 109
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
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);
André Anjos's avatar
André Anjos committed
125 126

  // compute position of the threshold
Amir Mohammadi's avatar
lint  
Amir Mohammadi committed
127
  double crr = 1. - far_value; // (Correct Rejection Rate; = 1 - FAR)
128
  double crr_index = crr * neg.extent(0) - 1.;
André Anjos's avatar
André Anjos committed
129
  // compute the index above the current CRR value
130
  int index = (int)std::ceil(crr_index);
André Anjos's avatar
André Anjos committed
131

132
  // increase the threshold when we have several negatives with the same score
133 134
  while (index < neg.extent(0)-1 && neg(index) == neg(index+1))
    ++index;
André Anjos's avatar
André Anjos committed
135

136
  if (index < neg.extent(0)-1){
137 138
    // return the threshold that is just above the desired FAR
    return neg(index);
André Anjos's avatar
André Anjos committed
139
  } else {
140 141
    // 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
142 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;
163
  sort(positives, pos, is_sorted);
André Anjos's avatar
André Anjos committed
164 165

  // compute position of the threshold
166
  double frr_index = frr_value * pos.extent(0) - 1.;
167 168
  // compute the index below the current FAR value
  int index = (int)std::ceil(frr_index);
André Anjos's avatar
André Anjos committed
169

170
  // lower the threshold when several positives have the same score
171 172
  while (index && pos(index) == pos(index-1))
    --index;
André Anjos's avatar
André Anjos committed
173

174
  if (index){
175 176 177
    // return the FRR threshold that is just above the desired FRR
    // We have to add a little noise to since the FRR calculation excludes the threshold
    return pos(index) + 1e-8 * pos(index);
André Anjos's avatar
André Anjos committed
178
  } else {
179 180
    // 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
181 182 183 184 185 186 187 188 189 190
  }
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  return retval;
}

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

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

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

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

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

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

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

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

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

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

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

  return retval;
}

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

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

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

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

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

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