Commit c2697a01 authored by Amir Mohammadi's avatar Amir Mohammadi
Browse files

lint

parent dfd1ea75
Pipeline #8620 passed with stages
in 7 minutes and 43 seconds
...@@ -66,12 +66,13 @@ def min_cllr(negatives, positives): ...@@ -66,12 +66,13 @@ def min_cllr(negatives, positives):
pos = sorted(positives) pos = sorted(positives)
N = len(neg) N = len(neg)
P = len(pos) P = len(pos)
I = N+P I = N + P
# now, iterate through both score sets and add a 0 for negative and 1 for positive scores # now, iterate through both score sets and add a 0 for negative and 1 for
n, p = 0,0 # positive scores
n, p = 0, 0
ideal = numpy.zeros(I) ideal = numpy.zeros(I)
neg_indices = [0]*N neg_indices = [0] * N
pos_indices = [0]*P pos_indices = [0] * P
for i in range(I): for i in range(I):
if p < P and (n == N or neg[n] > pos[p]): if p < P and (n == N or neg[n] > pos[p]):
pos_indices[p] = i pos_indices[p] = i
...@@ -88,12 +89,12 @@ def min_cllr(negatives, positives): ...@@ -88,12 +89,12 @@ def min_cllr(negatives, positives):
# disable runtime warnings for a short time since log(0) will raise a warning # disable runtime warnings for a short time since log(0) will raise a warning
old_warn_setup = numpy.seterr(divide='ignore') old_warn_setup = numpy.seterr(divide='ignore')
# ... compute logs # ... compute logs
posterior_log_odds = numpy.log(popt)-numpy.log(1.-popt); posterior_log_odds = numpy.log(popt) - numpy.log(1. - popt)
log_prior_odds = math.log(float(P)/float(N)); log_prior_odds = math.log(float(P) / float(N))
# ... activate old warnings # ... activate old warnings
numpy.seterr(**old_warn_setup) numpy.seterr(**old_warn_setup)
llrs = posterior_log_odds - log_prior_odds; llrs = posterior_log_odds - log_prior_odds
# some weired addition # some weired addition
# for i in range(I): # for i in range(I):
......
...@@ -7,153 +7,179 @@ ...@@ -7,153 +7,179 @@
* Copyright (C) Idiap Research Institute, Martigny, Switzerland * Copyright (C) Idiap Research Institute, Martigny, Switzerland
*/ */
#include <stdexcept>
#include <algorithm> #include <algorithm>
#include <limits>
#include <boost/format.hpp> #include <boost/format.hpp>
#include <limits>
#include <stdexcept>
#include <bob.core/assert.h>
#include <bob.core/cast.h>
#include <bob.core/array_copy.h> #include <bob.core/array_copy.h>
#include <bob.core/array_sort.h> #include <bob.core/array_sort.h>
#include <bob.core/assert.h>
#include <bob.core/cast.h>
#include <bob.math/pavx.h>
#include <bob.math/linsolve.h> #include <bob.math/linsolve.h>
#include <bob.math/pavx.h>
#include "error.h" #include "error.h"
template <typename T> template <typename T>
static void sort(const blitz::Array<T,1>& a, blitz::Array<T,1>& b, bool is_sorted){ static void sort(const blitz::Array<T, 1> &a, blitz::Array<T, 1> &b,
if (is_sorted){ bool is_sorted) {
if (is_sorted) {
b.reference(a); b.reference(a);
} else { } else {
bob::core::array::ccopy(a,b); bob::core::array::ccopy(a, b);
bob::core::array::sort<T>(b); bob::core::array::sort<T>(b);
} }
} }
std::pair<double, double> bob::measure::farfrr(const blitz::Array<double,1>& negatives, std::pair<double, double>
const blitz::Array<double,1>& positives, double threshold) { bob::measure::farfrr(const blitz::Array<double, 1> &negatives,
if (!negatives.size()) throw std::runtime_error("Cannot compute FAR when no negatives are given"); const blitz::Array<double, 1> &positives,
if (!positives.size()) throw std::runtime_error("Cannot compute FRR when no positives are given"); 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");
blitz::sizeType total_negatives = negatives.extent(blitz::firstDim); blitz::sizeType total_negatives = negatives.extent(blitz::firstDim);
blitz::sizeType total_positives = positives.extent(blitz::firstDim); blitz::sizeType total_positives = positives.extent(blitz::firstDim);
blitz::sizeType false_accepts = blitz::count(negatives >= threshold); blitz::sizeType false_accepts = blitz::count(negatives >= threshold);
blitz::sizeType false_rejects = blitz::count(positives < threshold); blitz::sizeType false_rejects = blitz::count(positives < threshold);
return std::make_pair(false_accepts/(double)total_negatives, return std::make_pair(false_accepts / (double)total_negatives,
false_rejects/(double)total_positives); false_rejects / (double)total_positives);
} }
std::pair<double, double> bob::measure::precision_recall(const blitz::Array<double,1>& negatives, std::pair<double, double>
const blitz::Array<double,1>& positives, double threshold) { bob::measure::precision_recall(const blitz::Array<double, 1> &negatives,
if (!negatives.size() || !positives.size()) throw std::runtime_error("Cannot compute precision or recall when no positives or no negatives are given"); 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");
blitz::sizeType total_positives = positives.extent(blitz::firstDim); blitz::sizeType total_positives = positives.extent(blitz::firstDim);
blitz::sizeType false_positives = blitz::count(negatives >= threshold); blitz::sizeType false_positives = blitz::count(negatives >= threshold);
blitz::sizeType true_positives = blitz::count(positives >= threshold); blitz::sizeType true_positives = blitz::count(positives >= threshold);
blitz::sizeType total_classified_positives = true_positives + false_positives; blitz::sizeType total_classified_positives = true_positives + false_positives;
if (!total_classified_positives) total_classified_positives = 1; //avoids division by zero if (!total_classified_positives)
if (!total_positives) total_positives = 1; //avoids division by zero total_classified_positives = 1; // avoids division by zero
return std::make_pair(true_positives/(double)(total_classified_positives), if (!total_positives)
true_positives/(double)(total_positives)); total_positives = 1; // avoids division by zero
return std::make_pair(true_positives / (double)(total_classified_positives),
true_positives / (double)(total_positives));
} }
double bob::measure::f_score(const blitz::Array<double, 1> &negatives,
double bob::measure::f_score(const blitz::Array<double,1>& negatives, const blitz::Array<double, 1> &positives,
const blitz::Array<double,1>& positives, double threshold, double weight) { double threshold, double weight) {
std::pair<double, double> ratios = std::pair<double, double> ratios =
bob::measure::precision_recall(negatives, positives, threshold); bob::measure::precision_recall(negatives, positives, threshold);
double precision = ratios.first; double precision = ratios.first;
double recall = ratios.second; double recall = ratios.second;
if (weight <= 0) weight = 1; if (weight <= 0)
weight = 1;
if (precision == 0 && recall == 0) if (precision == 0 && recall == 0)
return 0; return 0;
return (1 + weight*weight) * precision * recall / (weight * weight * precision + recall); return (1 + weight * weight) * precision * recall /
(weight * weight * precision + recall);
} }
double eer_predicate(double far, double frr) { double eer_predicate(double far, double frr) { return std::abs(far - frr); }
return std::abs(far - frr);
}
double bob::measure::eerThreshold(const blitz::Array<double,1>& negatives, const blitz::Array<double,1>& positives, bool is_sorted) { double bob::measure::eerThreshold(const blitz::Array<double, 1> &negatives,
blitz::Array<double,1> neg, pos; const blitz::Array<double, 1> &positives,
bool is_sorted) {
blitz::Array<double, 1> neg, pos;
sort(negatives, neg, is_sorted); sort(negatives, neg, is_sorted);
sort(positives, pos, is_sorted); sort(positives, pos, is_sorted);
return bob::measure::minimizingThreshold(neg, pos, eer_predicate); 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)); 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, bool is_sorted) { double bob::measure::farThreshold(const blitz::Array<double, 1> &negatives,
const blitz::Array<double, 1> &,
double far_value, bool is_sorted) {
// check the parameters are valid // check the parameters are valid
if (far_value < 0. || far_value > 1.) { 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.]"); boost::format m("the argument for `far_value' cannot take the value %f - "
"the value must be in the interval [0.,1.]");
m % far_value; m % far_value;
throw std::runtime_error(m.str()); throw std::runtime_error(m.str());
} }
if (negatives.size() < 2) { if (negatives.size() < 2) {
throw std::runtime_error("the number of negative scores must be at least 2"); throw std::runtime_error(
"the number of negative scores must be at least 2");
} }
// sort the array, if necessary // sort the array, if necessary
blitz::Array<double,1> neg; blitz::Array<double, 1> neg;
sort(negatives, neg, is_sorted); sort(negatives, neg, is_sorted);
// compute position of the threshold // compute position of the threshold
double crr = 1.-far_value; // (Correct Rejection Rate; = 1 - FAR) double crr = 1. - far_value; // (Correct Rejection Rate; = 1 - FAR)
double crr_index = crr * neg.extent(0); double crr_index = crr * neg.extent(0);
// compute the index above the current CRR value // compute the index above the current CRR value
int index = std::min((int)std::floor(crr_index), neg.extent(0)-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 // correct index if we have multiple score values at the requested position
while (index && neg(index) == neg(index-1)) --index; while (index && neg(index) == neg(index - 1))
--index;
// we compute a correction term to assure that we are in the middle of two cases // we compute a correction term to assure that we are in the middle of two
// cases
double correction; double correction;
if (index){ if (index) {
// assure that we are in the middle of two cases // assure that we are in the middle of two cases
correction = 0.5 * (neg(index) - neg(index-1)); correction = 0.5 * (neg(index) - neg(index - 1));
} else { } else {
// add an overall correction term // add an overall correction term
correction = 0.5 * (neg(neg.extent(0)-1) - neg(0)) / neg.extent(0); correction = 0.5 * (neg(neg.extent(0) - 1) - neg(0)) / neg.extent(0);
} }
return neg(index) - correction; return neg(index) - correction;
} }
double bob::measure::frrThreshold(const blitz::Array<double,1>&, const blitz::Array<double,1>& positives, double frr_value, bool is_sorted) { double bob::measure::frrThreshold(const blitz::Array<double, 1> &,
const blitz::Array<double, 1> &positives,
double frr_value, bool is_sorted) {
// check the parameters are valid // check the parameters are valid
if (frr_value < 0. || frr_value > 1.) { if (frr_value < 0. || frr_value > 1.) {
boost::format m("the argument for `frr_value' cannot take the value %f - the value must be in the interval [0.,1.]"); boost::format m("the argument for `frr_value' cannot take the value %f - "
"the value must be in the interval [0.,1.]");
m % frr_value; m % frr_value;
throw std::runtime_error(m.str()); throw std::runtime_error(m.str());
} }
if (positives.size() < 2) { if (positives.size() < 2) {
throw std::runtime_error("the number of positive scores must be at least 2"); throw std::runtime_error(
"the number of positive scores must be at least 2");
} }
// sort positive scores descendantly, if necessary // sort positive scores descendantly, if necessary
blitz::Array<double,1> pos; blitz::Array<double, 1> pos;
sort(positives, pos, is_sorted); sort(positives, pos, is_sorted);
// compute position of the threshold // compute position of the threshold
double frr_index = frr_value * pos.extent(0); double frr_index = frr_value * pos.extent(0);
// compute the index above the current CAR value // compute the index above the current CAR value
int index = std::min((int)std::ceil(frr_index), pos.extent(0)-1); 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 // correct index if we have multiple score values at the requested position
while (index < pos.extent(0)-1 && pos(index) == pos(index+1)) ++index; while (index < pos.extent(0) - 1 && 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 to assure that we are in the middle of two
// cases
double correction; double correction;
if (index < pos.extent(0)-1){ if (index < pos.extent(0) - 1) {
// assure that we are in the middle of two cases // assure that we are in the middle of two cases
correction = 0.5 * (pos(index+1) - pos(index)); correction = 0.5 * (pos(index + 1) - pos(index));
} else { } else {
// add an overall correction term // add an overall correction term
correction = 0.5 * (pos(pos.extent(0)-1) - pos(0)) / pos.extent(0); correction = 0.5 * (pos(pos.extent(0) - 1) - pos(0)) / pos.extent(0);
} }
return pos(index) + correction; return pos(index) + correction;
...@@ -166,54 +192,59 @@ class weighted_error { ...@@ -166,54 +192,59 @@ class weighted_error {
double m_weight; ///< The weighting factor double m_weight; ///< The weighting factor
public: //api public: // api
weighted_error(double weight) : m_weight(weight) {
weighted_error(double weight): m_weight(weight) { if (weight > 1.0)
if (weight > 1.0) m_weight = 1.0; m_weight = 1.0;
if (weight < 0.0) m_weight = 0.0; if (weight < 0.0)
m_weight = 0.0;
} }
inline double operator() (double far, double frr) const { inline double operator()(double far, double frr) const {
return (m_weight*far) + ((1.0-m_weight)*frr); return (m_weight * far) + ((1.0 - m_weight) * frr);
} }
}; };
double bob::measure::minWeightedErrorRateThreshold(const blitz::Array<double,1>& negatives, const blitz::Array<double,1>& positives, double cost, bool is_sorted) { double bob::measure::minWeightedErrorRateThreshold(
blitz::Array<double,1> neg, pos; const blitz::Array<double, 1> &negatives,
const blitz::Array<double, 1> &positives, double cost, bool is_sorted) {
blitz::Array<double, 1> neg, pos;
sort(negatives, neg, is_sorted); sort(negatives, neg, is_sorted);
sort(positives, pos, is_sorted); sort(positives, pos, is_sorted);
weighted_error predicate(cost); weighted_error predicate(cost);
return bob::measure::minimizingThreshold(neg, pos, predicate); return bob::measure::minimizingThreshold(neg, pos, predicate);
} }
blitz::Array<double,2> bob::measure::roc(const blitz::Array<double,1>& negatives, blitz::Array<double, 2>
const blitz::Array<double,1>& positives, size_t points) { bob::measure::roc(const blitz::Array<double, 1> &negatives,
const blitz::Array<double, 1> &positives, size_t points) {
double min = std::min(blitz::min(negatives), blitz::min(positives)); double min = std::min(blitz::min(negatives), blitz::min(positives));
double max = std::max(blitz::max(negatives), blitz::max(positives)); double max = std::max(blitz::max(negatives), blitz::max(positives));
double step = (max-min)/((double)points-1.0); double step = (max - min) / ((double)points - 1.0);
blitz::Array<double,2> retval(2, points); blitz::Array<double, 2> retval(2, points);
for (int i=0; i<(int)points; ++i) { for (int i = 0; i < (int)points; ++i) {
std::pair<double, double> ratios = std::pair<double, double> ratios =
bob::measure::farfrr(negatives, positives, min + i*step); bob::measure::farfrr(negatives, positives, min + i * step);
// preserve X x Y ordering (FAR x FRR) // preserve X x Y ordering (FAR x FRR)
retval(0,i) = ratios.first; retval(0, i) = ratios.first;
retval(1,i) = ratios.second; retval(1, i) = ratios.second;
} }
return retval; return retval;
} }
blitz::Array<double,2> bob::measure::precision_recall_curve(const blitz::Array<double,1>& negatives, blitz::Array<double, 2>
const blitz::Array<double,1>& positives, size_t points) { bob::measure::precision_recall_curve(const blitz::Array<double, 1> &negatives,
const blitz::Array<double, 1> &positives,
size_t points) {
double min = std::min(blitz::min(negatives), blitz::min(positives)); double min = std::min(blitz::min(negatives), blitz::min(positives));
double max = std::max(blitz::max(negatives), blitz::max(positives)); double max = std::max(blitz::max(negatives), blitz::max(positives));
double step = (max-min)/((double)points-1.0); double step = (max - min) / ((double)points - 1.0);
blitz::Array<double,2> retval(2, points); blitz::Array<double, 2> retval(2, points);
for (int i=0; i<(int)points; ++i) { for (int i = 0; i < (int)points; ++i) {
std::pair<double, double> ratios = std::pair<double, double> ratios =
bob::measure::precision_recall(negatives, positives, min + i*step); bob::measure::precision_recall(negatives, positives, min + i * step);
retval(0,i) = ratios.first; retval(0, i) = ratios.first;
retval(1,i) = ratios.second; retval(1, i) = ratios.second;
} }
return retval; return retval;
} }
...@@ -221,124 +252,115 @@ blitz::Array<double,2> bob::measure::precision_recall_curve(const blitz::Array<d ...@@ -221,124 +252,115 @@ blitz::Array<double,2> bob::measure::precision_recall_curve(const blitz::Array<d
/** /**
* Structure for getting permutations when sorting an array * Structure for getting permutations when sorting an array
*/ */
struct ComparePairs struct ComparePairs {
{ ComparePairs(const blitz::Array<double, 1> &v) : m_v(v) {}
ComparePairs(const blitz::Array<double,1> &v):
m_v(v)
{
}
bool operator()(size_t a, size_t b) bool operator()(size_t a, size_t b) { return m_v(a) < m_v(b); }
{
return m_v(a) < m_v(b);
}
blitz::Array<double,1> m_v; blitz::Array<double, 1> m_v;
}; };
/** /**
* Sort an array and get the permutations (using stable_sort) * Sort an array and get the permutations (using stable_sort)
*/ */
void sortWithPermutation(const blitz::Array<double,1>& values, std::vector<size_t>& v) void sortWithPermutation(const blitz::Array<double, 1> &values,
{ std::vector<size_t> &v) {
int N = values.extent(0); int N = values.extent(0);
bob::core::array::assertSameDimensionLength(N, v.size()); bob::core::array::assertSameDimensionLength(N, v.size());
for(int i=0; i<N; ++i) for (int i = 0; i < N; ++i)
v[i] = i; v[i] = i;
std::stable_sort(v.begin(), v.end(), ComparePairs(values)); std::stable_sort(v.begin(), v.end(), ComparePairs(values));
} }
blitz::Array<double,2> bob::measure::rocch(const blitz::Array<double,1>& negatives, blitz::Array<double, 2>
const blitz::Array<double,1>& positives) bob::measure::rocch(const blitz::Array<double, 1> &negatives,
{ const blitz::Array<double, 1> &positives) {
// Number of positive and negative scores // Number of positive and negative scores
size_t Nt = positives.extent(0); size_t Nt = positives.extent(0);
size_t Nn = negatives.extent(0); size_t Nn = negatives.extent(0);
size_t N = Nt + Nn; size_t N = Nt + Nn;
// Create a big array with all scores // Create a big array with all scores
blitz::Array<double,1> scores(N); blitz::Array<double, 1> scores(N);
blitz::Range rall = blitz::Range::all(); blitz::Range rall = blitz::Range::all();
scores(blitz::Range(0,Nt-1)) = positives(rall); scores(blitz::Range(0, Nt - 1)) = positives(rall);
scores(blitz::Range(Nt,N-1)) = negatives(rall); scores(blitz::Range(Nt, N - 1)) = negatives(rall);
// It is important here that scores that are the same (i.e. already in order) should NOT be swapped. // It is important here that scores that are the same (i.e. already in order)
// should NOT be swapped.
// std::stable_sort has this property. // std::stable_sort has this property.
std::vector<size_t> perturb(N); std::vector<size_t> perturb(N);
sortWithPermutation(scores, perturb); sortWithPermutation(scores, perturb);
// Apply permutation // Apply permutation
blitz::Array<size_t,1> Pideal(N); blitz::Array<size_t, 1> Pideal(N);
for(size_t i=0; i<N; ++i) for (size_t i = 0; i < N; ++i)
Pideal(i) = (perturb[i] < Nt ? 1 : 0); Pideal(i) = (perturb[i] < Nt ? 1 : 0);
blitz::Array<double,1> Pideal_d = bob::core::array::cast<double>(Pideal); blitz::Array<double, 1> Pideal_d = bob::core::array::cast<double>(Pideal);
// Apply the PAVA algorithm // Apply the PAVA algorithm
blitz::Array<double,1> Popt(N); blitz::Array<double, 1> Popt(N);
blitz::Array<size_t,1> width = bob::math::pavxWidth(Pideal_d, Popt); blitz::Array<size_t, 1> width = bob::math::pavxWidth(Pideal_d, Popt);
// Allocate output // Allocate output
int nbins = width.extent(0); int nbins = width.extent(0);
blitz::Array<double,2> retval(2,nbins+1); // FAR, FRR blitz::Array<double, 2> retval(2, nbins + 1); // FAR, FRR
// Fill in output // Fill in output
size_t left = 0; size_t left = 0;
size_t fa = Nn; size_t fa = Nn;
size_t miss = 0; size_t miss = 0;
for(int i=0; i<nbins; ++i) for (int i = 0; i < nbins; ++i) {
{ retval(0, i) = fa / (double)Nn; // pfa
retval(0,i) = fa / (double)Nn; // pfa