Merge branch 'new_error'

Conflicts:
	bob/measure/cpp/error.h

bob/measure/__init__.py

@@ -13,7 +13,9 @@ from . import openbr
 import numpy
 
 def mse (estimation, target):
-  """Calculates the mean square error between a set of outputs and target
+  """mse(estimation, target) -> error
+
+  Calculates the mean square error between a set of outputs and target
   values using the following formula:
 
   .. math::
@@ -28,7 +30,9 @@ def mse (estimation, target):
   return numpy.mean((estimation - target)**2, 0)
 
 def rmse (estimation, target):
-  """Calculates the root mean square error between a set of outputs and target
+  """rmse(estimation, target) -> error
+
+  Calculates the root mean square error between a set of outputs and target
   values using the following formula:
 
   .. math::
@@ -43,12 +47,14 @@ def rmse (estimation, target):
   return numpy.sqrt(mse(estimation, target))
 
 def relevance (input, machine):
-  """Calculates the relevance of every input feature to the estimation process
+  """relevance (input, machine) -> relevances
+
+  Calculates the relevance of every input feature to the estimation process
   using the following definition from:
 
-  Neural Triggering System Operating on High Resolution Calorimetry
-  Information, Anjos et al, April 2006, Nuclear Instruments and Methods in
-  Physics Research, volume 559, pages 134-138
+    Neural Triggering System Operating on High Resolution Calorimetry
+    Information, Anjos et al, April 2006, Nuclear Instruments and Methods in
+    Physics Research, volume 559, pages 134-138
 
   .. math::
 
@@ -73,10 +79,12 @@ def relevance (input, machine):
   return retval
 
 def recognition_rate(cmc_scores):
-  """Calculates the recognition rate from the given input, which is identical
+  """recognition_rate(cmc_scores) -> RR
+
+  Calculates the recognition rate from the given input, which is identical
   to the rank 1 (C)MC value.
 
-  The input has a specific format, which is a list of two-element tuples.  Each
+  The input has a specific format, which is a list of two-element tuples. Each
   of the tuples contains the negative and the positive scores for one test
   item.  To read the lists from score files in 4 or 5 column format, please use
   the :py:func:`bob.measure.load.cmc_four_column` or
@@ -86,10 +94,20 @@ def recognition_rate(cmc_scores):
   positive score is greater than or equal to all negative scores, divided by
   the number of all test items.  If several positive scores for one test item
   exist, the **highest** score is taken.
+
+  **Parameters:**
+
+  ``cmc_scores`` : [(array_like(1D, float), array_like(1D, float))]
+    A list of tuples, where each tuple contains the ``negative`` and ``positive`` scores for one probe of the database
+
+  **Returns:**
+
+  ``RR`` : float
+    The rank 1 recognition rate, i.e., the relative number of correctly identified identities
   """
   # If no scores are given, the recognition rate is exactly 0.
   if not cmc_scores:
-    return 0
+    return 0.
 
   correct = 0.
   for neg, pos in cmc_scores:
@@ -98,15 +116,17 @@ def recognition_rate(cmc_scores):
     max_pos = numpy.max(pos)
     # check if the positive score is smaller than all negative scores
     if (neg < max_pos).all():
-      correct += 1
+      correct += 1.
 
   # return relative number of correctly matched scores
   return correct / float(len(cmc_scores))
 
 def cmc(cmc_scores):
-  """Calculates the cumulative match characteristic (CMC) from the given input.
+  """cmc(cmc_scores) -> curve
+
+  Calculates the cumulative match characteristic (CMC) from the given input.
 
-  The input has a specific format, which is a list of two-element tuples.  Each
+  The input has a specific format, which is a list of two-element tuples. Each
   of the tuples contains the negative and the positive scores for one test
   item.  To read the lists from score files in 4 or 5 column format, please use
   the :py:func:`bob.measure.load.cmc_four_column` or
@@ -117,6 +137,16 @@ def cmc(cmc_scores):
   higher than the positive score.  If several positive scores for one test item
   exist, the **highest** positive score is taken. The CMC finally computes how
   many test items have rank r or higher.
+
+  **Parameters:**
+
+  ``cmc_scores`` : [(array_like(1D, float), array_like(1D, float))]
+    A list of tuples, where each tuple contains the ``negative`` and ``positive`` scores for one probe of the database
+
+  **Returns:**
+
+  ``curve`` : array_like(2D, float)
+    The CMC curve, with the Rank in the first column and the number of correctly classified clients (in this rank) in the second column.
   """
 
   # If no scores are given, we cannot plot anything
@@ -147,6 +177,7 @@ def cmc(cmc_scores):
 def get_config():
   """Returns a string containing the configuration information.
   """
+  import bob.extension
   return bob.extension.get_config(__name__, version.externals)
 
 

bob/measure/calibration.py

@@ -11,7 +11,20 @@ import math
 import numpy
 
 def cllr(negatives, positives):
-  """Computes the 'cost of log likelihood ratio' measure as given in the bosaris toolkit"""
+  """cllr(negatives, positives) -> cllr
+
+  Computes the 'cost of log likelihood ratio' (:math:`C_{llr}`) measure as given in the Bosaris toolkit
+
+  **Parameters:**
+
+  ``negatives, positives`` : array_like(1D, float)
+    The scores computed by comparing elements from different classes and the same class, respectively.
+
+  **Returns**
+
+  ``cllr`` : float
+    The computed :math:`C_{llr}` value.
+  """
   sum_pos, sum_neg = 0., 0.
   for pos in positives:
     sum_pos += math.log(1. + math.exp(-pos), 2.)
@@ -21,7 +34,20 @@ def cllr(negatives, positives):
 
 
 def min_cllr(negatives, positives):
-  """Computes the 'minimum cost of log likelihood ratio' measure as given in the bosaris toolkit"""
+  """cllr(negatives, positives) -> cllr
+
+  Computes the 'minimum cost of log likelihood ratio' (:math:`C_{llr}^{min}`) measure as given in the bosaris toolkit
+
+  **Parameters:**
+
+  ``negatives, positives`` : array_like(1D, float)
+    The scores computed by comparing elements from different classes and the same class, respectively.
+
+  **Returns**
+
+  ``min_cllr`` : float
+    The computed :math:`C_{llr}^{min}` value.
+  """
 
   from bob.math import pavx
 

bob/measure/error.cpp → bob/measure/cpp/error.cpp

@@ -7,8 +7,6 @@
  * Copyright (C) Idiap Research Institute, Martigny, Switzerland
  */
 
-#include "error.h"
-
 #include <stdexcept>
 #include <algorithm>
 #include <limits>
@@ -16,10 +14,24 @@
 
 #include <bob.core/assert.h>
 #include <bob.core/cast.h>
+#include <bob.core/array_copy.h>
+#include <bob.core/array_sort.h>
 
 #include <bob.math/pavx.h>
 #include <bob.math/linsolve.h>
 
+#include "error.h"
+
+template <typename T>
+static void sort(const blitz::Array<T,1>& a, blitz::Array<T,1>& b, bool isSorted){
+  if (isSorted){
+    b.reference(a);
+  } else {
+    bob::core::array::ccopy(a,b);
+    bob::core::array::sort<T>(b);
+  }
+}
+
 std::pair<double, double> bob::measure::farfrr(const blitz::Array<double,1>& negatives,
     const blitz::Array<double,1>& positives, double threshold) {
   blitz::sizeType total_negatives = negatives.extent(blitz::firstDim);
@@ -61,18 +73,18 @@ double eer_predicate(double far, double frr) {
   return std::abs(far - frr);
 }
 
-double bob::measure::eerThreshold(const blitz::Array<double,1>& negatives,
-    const blitz::Array<double,1>& positives) {
-  return bob::measure::minimizingThreshold(negatives, positives, eer_predicate);
+double bob::measure::eerThreshold(const blitz::Array<double,1>& negatives, const blitz::Array<double,1>& positives, bool isSorted) {
+  blitz::Array<double,1> neg, pos;
+  sort(negatives, neg, isSorted);
+  sort(positives, pos, isSorted);
+  return bob::measure::minimizingThreshold(neg, pos, eer_predicate);
 }
 
-double bob::measure::eerRocch(const blitz::Array<double,1>& negatives,
-    const blitz::Array<double,1>& positives) {
+double bob::measure::eerRocch(const blitz::Array<double,1>& negatives, const blitz::Array<double,1>& positives) {
   return bob::measure::rocch2eer(bob::measure::rocch(negatives, positives));
 }
 
-double bob::measure::farThreshold(const blitz::Array<double,1>& negatives,
-  const blitz::Array<double,1>&, double far_value) {
+double bob::measure::farThreshold(const blitz::Array<double,1>& negatives, const blitz::Array<double,1>&, double far_value, bool isSorted) {
   // check the parameters are valid
   if (far_value < 0. || far_value > 1.) {
     boost::format m("the argument for `far_value' cannot take the value %f - the value must be in the interval [0.,1.]");
@@ -83,35 +95,33 @@ double bob::measure::farThreshold(const blitz::Array<double,1>& negatives,
     throw std::runtime_error("the number of negative scores must be at least 2");
   }
 
-  // sort negative scores ascendingly
-  std::vector<double> negatives_(negatives.shape()[0]);
-  std::copy(negatives.begin(), negatives.end(), negatives_.begin());
-  std::sort(negatives_.begin(), negatives_.end(), std::less<double>());
+  // sort the array, if necessary
+  blitz::Array<double,1> neg;
+  sort(negatives, neg, isSorted);
 
   // compute position of the threshold
   double crr = 1.-far_value; // (Correct Rejection Rate; = 1 - FAR)
-  double crr_index = crr * negatives_.size();
+  double crr_index = crr * neg.extent(0);
   // compute the index above the current CRR value
-  int index = std::min((int)std::floor(crr_index), (int)negatives_.size()-1);
+  int index = std::min((int)std::floor(crr_index), neg.extent(0)-1);
 
   // correct index if we have multiple score values at the requested position
-  while (index && negatives_[index] == negatives_[index-1]) --index;
+  while (index && neg(index) == neg(index-1)) --index;
 
-  // we compute a correction term
+  // we compute a correction term to assure that we are in the middle of two cases
   double correction;
   if (index){
     // assure that we are in the middle of two cases
-    correction = 0.5 * (negatives_[index] - negatives_[index-1]);
+    correction = 0.5 * (neg(index) - neg(index-1));
   } else {
     // add an overall correction term
-    correction = 0.5 * (negatives_.back() - negatives_.front()) / negatives_.size();
+    correction = 0.5 * (neg(neg.extent(0)-1) - neg(0)) / neg.extent(0);
   }
 
-  return negatives_[index] - correction;
+  return neg(index) - correction;
 }
 
-double bob::measure::frrThreshold(const blitz::Array<double,1>&,
-  const blitz::Array<double,1>& positives, double frr_value) {
+double bob::measure::frrThreshold(const blitz::Array<double,1>&, const blitz::Array<double,1>& positives, double frr_value, bool isSorted) {
 
   // check the parameters are valid
   if (frr_value < 0. || frr_value > 1.) {
@@ -123,32 +133,29 @@ double bob::measure::frrThreshold(const blitz::Array<double,1>&,
     throw std::runtime_error("the number of positive scores must be at least 2");
   }
 
-  // sort positive scores descendingly
-  std::vector<double> positives_(positives.shape()[0]);
-  std::copy(positives.begin(), positives.end(), positives_.begin());
-  std::sort(positives_.begin(), positives_.end(), std::greater<double>());
+  // sort positive scores descendantly, if necessary
+  blitz::Array<double,1> pos;
+  sort(positives, pos, isSorted);
 
   // compute position of the threshold
-  double car = 1.-frr_value; // (Correct Acceptance Rate; = 1 - FRR)
-  double car_index = car * positives_.size();
-  // compute the index above the current CRR value
-  int index = std::min((int)std::floor(car_index), (int)positives_.size()-1);
+  double frr_index = frr_value * pos.extent(0);
+  // compute the index above the current CAR value
+  int index = std::min((int)std::ceil(frr_index), pos.extent(0)-1);
 
   // correct index if we have multiple score values at the requested position
-  while (index && positives_[index] == positives_[index-1]) --index;
+  while (index < pos.extent(0) && pos(index) == pos(index+1)) ++index;
 
   // we compute a correction term to assure that we are in the middle of two cases
-  // we compute a correction term
   double correction;
-  if (index){
+  if (index < pos.extent(0)-1){
     // assure that we are in the middle of two cases
-    correction = 0.5 * (positives_[index-1] - positives_[index]);
+    correction = 0.5 * (pos(index+1) - pos(index));
   } else {
     // add an overall correction term
-    correction = 0.5 * (positives_.front() - positives_.back()) / positives_.size();
+    correction = 0.5 * (pos(pos.extent(0)-1) - pos(0)) / pos.extent(0);
   }
 
-  return positives_[index] + correction;
+  return pos(index) + correction;
 }
 
 /**
@@ -171,11 +178,12 @@ class weighted_error {
 
 };
 
-double bob::measure::minWeightedErrorRateThreshold
-(const blitz::Array<double,1>& negatives,
- const blitz::Array<double,1>& positives, double cost) {
+double bob::measure::minWeightedErrorRateThreshold(const blitz::Array<double,1>& negatives, const blitz::Array<double,1>& positives, double cost, bool isSorted) {
+  blitz::Array<double,1> neg, pos;
+  sort(negatives, neg, isSorted);
+  sort(positives, pos, isSorted);
   weighted_error predicate(cost);
-  return bob::measure::minimizingThreshold(negatives, positives, predicate);
+  return bob::measure::minimizingThreshold(neg, pos, predicate);
 }
 
 blitz::Array<double,2> bob::measure::roc(const blitz::Array<double,1>& negatives,
@@ -227,11 +235,6 @@ struct ComparePairs
   blitz::Array<double,1> m_v;
 };
 
-ComparePairs CreateComparePairs(const blitz::Array<double,1>& v)
-{
-  return ComparePairs(v);
-}
-
 /**
   * Sort an array and get the permutations (using stable_sort)
   */
@@ -242,7 +245,7 @@ void sortWithPermutation(const blitz::Array<double,1>& values, std::vector<size_
   for(int i=0; i<N; ++i)
     v[i] = i;
 
-  std::stable_sort(v.begin(), v.end(), CreateComparePairs(values));
+  std::stable_sort(v.begin(), v.end(), ComparePairs(values));
 }
 
 blitz::Array<double,2> bob::measure::rocch(const blitz::Array<double,1>& negatives,
@@ -358,20 +361,16 @@ double bob::measure::rocch2eer(const blitz::Array<double,2>& pfa_pmiss)
  * @return The ROC curve with the FAR in the first row and the FRR in the second.
  */
 blitz::Array<double,2> bob::measure::roc_for_far(const blitz::Array<double,1>& negatives,
- const blitz::Array<double,1>& positives, const blitz::Array<double,1>& far_list) {
+ const blitz::Array<double,1>& positives, const blitz::Array<double,1>& far_list, bool isSorted) {
   int n_points = far_list.extent(0);
 
   if (negatives.extent(0) == 0) throw std::runtime_error("The given set of negatives is empty.");
   if (positives.extent(0) == 0) throw std::runtime_error("The given set of positives is empty.");
 
-  // sort negative scores ascendingly
-  std::vector<double> negatives_(negatives.extent(0));;
-  std::copy(negatives.begin(), negatives.end(), negatives_.begin());
-  std::sort(negatives_.begin(), negatives_.end());
-  // sort positive scores ascendingly
-  std::vector<double> positives_(positives.extent(0));;
-  std::copy(positives.begin(), positives.end(), positives_.begin());
-  std::sort(positives_.begin(), positives_.end());
+  // sort negative and positive scores ascendantly
+  blitz::Array<double,1> neg, pos;
+  sort(negatives, neg, isSorted);
+  sort(positives, pos, isSorted);
 
   // do some magic to compute the FRR list
   blitz::Array<double,2> retval(2, n_points);
@@ -379,10 +378,10 @@ blitz::Array<double,2> bob::measure::roc_for_far(const blitz::Array<double,1>& n
   // index into the FAR and FRR list
   int far_index = n_points-1;
   int pos_index = 0, neg_index = 0;
-  int n_pos = positives_.size(), n_neg = negatives_.size();
+  int n_pos = pos.extent(0), n_neg = neg.extent(0);
 
   // iterators into the result lists
-  std::vector<double>::const_iterator pos_it = positives_.begin(), neg_it = negatives_.begin();
+  auto pos_it = pos.begin(), neg_it = neg.begin();
   // do some fast magic to compute the FRR values ;-)
   do{
     // check whether the current positive value is less than the current negative one
@@ -409,13 +408,13 @@ blitz::Array<double,2> bob::measure::roc_for_far(const blitz::Array<double,1>& n
     }
 
   // do this, as long as there are elements in both lists left and not all FRR elements where calculated yet
-  } while (pos_it != positives_.end() && neg_it != negatives_.end() && far_index >= 0);
+  } while (pos_it != pos.end() && neg_it != neg.end() && far_index >= 0);
 
   // check if all FRR values have been set
   if (far_index >= 0){
     // walk to the end of both lists; at least one of both lists should already have reached its limit.
-    pos_index += positives_.end() - pos_it;
-    neg_index += negatives_.end() - neg_it;
+    while (pos_it++ != pos.end()) ++pos_index;
+    while (neg_it++ != neg.end()) ++neg_index;
     // fill in the remaining elements of the CAR list
     do {
       // copy the FAR value
@@ -511,14 +510,19 @@ 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) {
+ const blitz::Array<double,1>& test_positives, size_t points, bool isSorted) {
+
+  blitz::Array<double,1> dev_neg, dev_pos;
+  sort(dev_negatives, dev_neg, isSorted);
+  sort(dev_positives, dev_pos, isSorted);
+
   double step = 1.0/((double)points-1.0);
   blitz::Array<double,2> retval(2, points);
   for (int i=0; i<(int)points; ++i) {
     double alpha = (double)i*step;
     retval(0,i) = alpha;
-    double threshold = bob::measure::minWeightedErrorRateThreshold(dev_negatives,
-        dev_positives, alpha);
+    double threshold = bob::measure::minWeightedErrorRateThreshold(dev_neg,
+        dev_pos, alpha, true);
     std::pair<double, double> ratios =
       bob::measure::farfrr(test_negatives, test_positives, threshold);
     retval(1,i) = (ratios.first + ratios.second) / 2;

bob/measure/error.h → bob/measure/cpp/error.h

@@ -62,7 +62,7 @@ namespace bob { namespace measure {
    * or "client"). 'negatives' holds the score information for samples that are
    * labeled *not* to belong to the class (a.k.a., "noise" or "impostor").
    *
-   * For more precise details about how the method considers error rates, please refer to the documentation of the method bob.measure.farfrr. 
+   * For more precise details about how the method considers error rates, please refer to the documentation of the method bob.measure.farfrr.
    *
    * It is possible that scores are inverted in the negative/positive sense. In
    * some setups the designer may have setup the system so 'positive' samples
@@ -107,66 +107,12 @@ namespace bob { namespace measure {
       return blitz::Array<bool,1>(negatives < threshold);
     }
 
-  /**
-   * Recursively minimizes w.r.t. to the given predicate method. Please refer
-   * to minimizingThreshold() for a full explanation. This method is only
-   * supposed to be used through that method.
-   */
-  template <typename T>
-  static double recursive_minimization(const blitz::Array<double,1>& negatives,
-      const blitz::Array<double,1>& positives, T& predicate,
-      double min, double max, size_t steps) {
-    static const double QUIT_THRESHOLD = 1e-10;
-    const double diff = max - min;
-    const double too_small = std::abs(diff/max);
-
-    //if the difference between max and min is too small, we quit.
-    if ( too_small < QUIT_THRESHOLD ) return min; //or max, does not matter...
-
-    double step_size = diff/(double)steps;
-    double min_value = predicate(1.0, 0.0); ///< to the left of the range
-
-    //the accumulator holds the thresholds that given the minimum value for the
-    //input predicate.
-    std::vector<double> accumulator;
-    accumulator.reserve(steps);
-
-    for (size_t i=0; i<steps; ++i) {
-      double threshold = ((double)i * step_size) + min;
-
-      std::pair<double, double> ratios =
-        farfrr(negatives, positives, threshold);
-
-      double current_cost = predicate(ratios.first, ratios.second);
-
-      if (current_cost < min_value) {
-        min_value = current_cost;
-        accumulator.clear(); ///< clean-up, we got a better minimum
-        accumulator.push_back(threshold); ///< remember this threshold
-      }
-      else if (std::abs(current_cost - min_value) < 1e-16) {
-        //accumulate to later decide...
-        accumulator.push_back(threshold);
-      }
-    }
-
-    //we stop when it doesn't matter anymore to threshold.
-    if (accumulator.size() != steps) {
-      //still needs some refinement: pick-up the middle of the range and go
-      return recursive_minimization(negatives, positives, predicate,
-          accumulator[accumulator.size()/2]-step_size,
-          accumulator[accumulator.size()/2]+step_size,
-          steps);
-    }
-
-    return accumulator[accumulator.size()/2];
-  }
-
   /**
    * This method can calculate a threshold based on a set of scores (positives
    * and negatives) given a certain minimization criteria, input as a
    * functional predicate. For a discussion on 'positive' and 'negative' see
    * bob::measure::farfrr().
+   * Here, it is expected that the positives and the negatives are sorted ascendantly.
    *
    * The predicate method gives back the current minimum given false-acceptance
    * (FA) and false-rejection (FR) ratios for the input data. As a predicate,
@@ -178,46 +124,100 @@ namespace bob { namespace measure {
    * Please note that this method will only work with single-minimum smooth
    * predicates.
    *
-   * The minimization is carried out in a recursive manner. First, we identify
-   * the threshold that minimizes the predicate given a set of N (N=100)
-   * thresholds between the min(negatives, positives) and the max(negatives,
-   * positives). If the minimum lies in a range of values, the center value is
-   * picked up.
+   * The minimization is carried out in a data-driven way.
+   * Starting from the lowest score (might be a positive or a negative), it
+   * increases the threshold based on the distance between the current score
+   * and the following higher score (also keeping track of duplicate scores)
+   * and computes the predicate for each possible threshold.
    *
-   * In a second round of minimization new minimum and maximum bounds are
-   * defined based on the center value plus/minus the step (max-min/N) and a
-   * new minimization round is restarted for N samples within the new bounds.
-   *
-   * The procedure continues until all calculated predicates in a given round
-   * give the same minimum. At this point, the center threshold is picked up and
-   * returned.
+   * Finally, that threshold is returned, for which the predicate returned the
+   * lowest value.
    */
-  template <typename T> double
-    minimizingThreshold(const blitz::Array<double,1>& negatives,
-        const blitz::Array<double,1>& positives, T& predicate) {
-      const size_t N = 100; ///< number of steps in each iteration
-      double min = std::min(blitz::min(negatives), blitz::min(positives));
-      double max = std::max(blitz::max(negatives), blitz::max(positives));
-      return recursive_minimization(negatives, positives, predicate, min,
-          max, N);
+  template <typename T>
+  double minimizingThreshold(const blitz::Array<double,1>& negatives, const blitz::Array<double,1>& positives, T& predicate){
+    // iterate over the whole set of points
+    auto pos_it = positives.begin(), neg_it = negatives.begin();
+
+    // iterate over all possible far and frr points and compute the predicate for each possible threshold...
+    double min_predicate = 1e8;
+    double min_threshold = 1e8;
+    double current_predicate = 1e8;
+    // we start with the extreme values for far and frr
+    double far = 1., frr = 0.;
+    // the decrease/increase for far/frr when moving one negative/positive
+    double far_decrease = 1./negatives.extent(0), frr_increase = 1./positives.extent(0);
+    // we start with the threshold based on the minimum score
+    double current_threshold = std::min(*pos_it, *neg_it);
+    // now, iterate over both lists, in a sorted order
+    while (pos_it != positives.end() && neg_it != negatives.end()){
+      // compute predicate
+      current_predicate = predicate(far, frr);
+      if (current_predicate <= min_predicate){
+        min_predicate = current_predicate;
+        min_threshold = current_threshold;
+      }
+      if (*pos_it >= *neg_it){
+        // compute current threshold
+        current_threshold = *neg_it;
+        // go to the next negative value
+        ++neg_it;
+        far -= far_decrease;
+      } else {
+        // compute current threshold
+        current_threshold = *pos_it;
+        // go to the next positive
+        ++pos_it;
+        frr += frr_increase;
+      }
+      // increase positive and negative as long as they contain the same value
+      while (neg_it != negatives.end() && *neg_it == current_threshold) {
+        // go to next negative
+        ++neg_it;
+        far -= far_decrease;
+      }
+      while (pos_it != positives.end() && *pos_it == current_threshold) {
+        // go to next positive
+        ++pos_it;
+        frr += frr_increase;
+      }
+      // compute a new threshold based on the center between last and current score, if we are not already at the end of the score lists
+      if (neg_it != negatives.end() || pos_it != positives.end()){
+        if (neg_it != negatives.end() && pos_it != positives.end())
+          current_threshold += std::min(*pos_it, *neg_it);
+        else if (neg_it != negatives.end())
+          current_threshold += *neg_it;
+        else
+          current_threshold += *pos_it;
+        current_threshold /= 2;
+      }
+    } // while
+
+    // now, we have reached the end of one list (usually the negatives)
+    // so, finally compute predicate for the last time
+    current_predicate = predicate(far, frr);
+    if (current_predicate < min_predicate){
+      min_predicate = current_predicate;
+      min_threshold = current_threshold;
     }