Changeset 4037

Timestamp:

Jan 25, 2021, 5:08:28 AM (5 months ago)

Author:

Peter

Message:

skeleton for new Ranking container. refs #710.

Location:

branches/kendall-score

Files:

• test/Makefile.am (modified) (1 diff)
• test/kendall.cc (modified) (3 diffs)
• test/ranking.cc (added)
• yat/statistics/Kendall.cc (modified) (3 diffs)
• yat/utility/Makefile.am (modified) (3 diffs)
• yat/utility/Ranking.cc (added)
• yat/utility/Ranking.h (added)

branches/kendall-score/test/Makefile.am

@@ -101 +101 @@
   test/poisson.test \
   test/queue.test test/queue2.test \
-  test/random_distribution.test \
-  test/range.test test/regression.test test/rnd.test \
+  test/range.test \
+  test/ranking.test \
+  test/regression.test test/rnd.test \
   test/rng-mt.test \
   test/roc.test \

branches/kendall-score/test/kendall.cc

@@ -27 +27 @@
 #include "yat/statistics/Kendall.h"
 #include "yat/utility/Matrix.h"
+#include "yat/utility/utility.h"
 
 #include <gsl/gsl_cdf.h>
@@ -47 +48 @@
 {
   test::Suite suite(argc, argv);
-
-  test_add(suite);
-  test_copy(suite);
-  test_p(suite);
-  test_p_exact(suite);
-  test_p_with_ties(suite);
-  test_score(suite);
-  test_score_with_ties(suite);
+  using std::throw_with_nested;
+  try {
+    try {
+      test_add(suite);
+    }
+    catch (std::exception& e) {
+      throw_with_nested(std::runtime_error("test_add() failed"));
+    }
+    test_copy(suite);
+    test_p(suite);
+    test_p_exact(suite);
+    test_p_with_ties(suite);
+    try {
+      test_score(suite);
+    }
+    catch (std::exception& e) {
+      throw_with_nested(std::runtime_error("test_score() failed"));
+    }
+    try {
+      test_score_with_ties(suite);
+    }
+    catch (std::exception& e) {
+      throw_with_nested(std::runtime_error("test_score_with_ties() failed"));
+    }
+  }
+  catch (std::exception& e) {
+    suite.add(false);
+    suite.err() << "error: exception caught: what(): ";
+    utility::print_what(e, suite.err());
+  }
   return suite.return_value();
 }
@@ -62 +85 @@
 {
   suite.out() << YAT_TEST_PROLOGUE;
-  Kendall kendall;
+  suite.err() << "constructor\n";
+  Kendall kendall;
+  suite.err() << "::add(0,0)\n";
   kendall.add(0,0);
-  kendall.add(1,1);
-  double score = kendall.score();
-  if (!suite.add(suite.equal(score, 1.0)))
-    suite.err() << "score not 1.0\n";
+  suite.err() << "::add(1,1)\n";
+  kendall.add(1,1);
+  suite.err() << "::score()\n";
+  try {
+    double score = kendall.score();
+    if (!suite.add(suite.equal(score, 1.0)))
+      suite.err() << "score not 1.0\n";
+  }
+  catch (std::exception& e) {
+    std::throw_with_nested(std::runtime_error("::score(void) failed"));
+  }
 }
 

branches/kendall-score/yat/statistics/Kendall.cc

@@ -24 +24 @@
 #include "Kendall.h"
 
+#include <yat/utility/Ranking.h>
+#include <yat/utility/stl_utility.h>
+
 #include <gsl/gsl_cdf.h>
+
+#include <boost/scoped_ptr.hpp>
 
 #include <algorithm>
 #include <cassert>
 #include <cmath>
+#include <iterator>
 #include <limits>
 #include <map>
+#include <set>
 #include <vector>
 
+
+#include <iostream> // debug
 namespace theplu {
 namespace yat {
@@ -66 +75 @@
   class Kendall::Pimpl
   {
+    class Count
+    {
+    public:
+      Count(const std::multiset<std::pair<double, double> >& data);
+      long int count(void) const;
+      double score(void) const;
+      double variance(void) const;
+      class Ties
+      {
+      public:
+        Ties(void);
+        void add(size_t n);
+        bool have_ties(void) const { return n_pairs_; }
+        // \return \sum x * (x-1)
+        unsigned long int n_pairs(void) const { return n_pairs_;}
+        // \return \sum x * (x-1) * (x-2)
+        unsigned long int n_triples(void) const { return n_triples_; }
+        // \return \sum x * (x-1) * (2*x+5)
+        unsigned long int v_correction(void) const { return v_correction_; }
+      private:
+        unsigned long int n_pairs_;
+        unsigned long int n_triples_;
+        unsigned long int v_correction_;
+      };
+
+    private:
+      Ties x_ties_;
+      Ties y_ties_;
+
+      // # pairs such that (x_i-x_j)(y_i-y_j) > 0
+      long int concordant_;
+      // # pairs such that (x_i-x_j)(y_i-y_j) < 0
+      long int disconcordant_;
+      // # pairs such that x_i!=x_j && y_i==y_j
+      long int extraX_;
+      // # pairs such that x_i==x_j && y_i!=y_j
+      long int extraY_;
+      // # pairs such that x_i==x_j && y_i==y_j
+      long int spare_;
+
+      template<typename Iterator>
+      void calculate_ties(Iterator first, Iterator last, Ties& ties)
+      {
+        while (first != last) {
+          Iterator upper = first;
+          size_t n = 1;
+          ++upper;
+          while (upper!=last && *upper==*first) {
+            ++n;
+            ++upper;
+          }
+          ties.add(n);
+          first = upper;
+        }
+      }
+    };
+
   public:
-    Pimpl(void) : updated_(false) {}
-
-    void add(double x, double y)
-    {
-      x_.push_back(x);
-      std::map<double, size_t>::iterator it = x_count_.lower_bound(x);
-      if (it==x_count_.end() || it->first != x)
-        x_count_.insert(it, std::make_pair(x, 1));
-      else
-        ++it->second;
-      y_.push_back(y);
-      it = y_count_.lower_bound(y);
-      if (it==y_count_.end() || it->first != y)
-        y_count_.insert(it, std::make_pair(y, 1));
-      else
-        ++it->second;
-      updated_ = false;
+    Pimpl(void);
+    Pimpl(const Pimpl& other);
+    Pimpl& operator=(const Pimpl& rhs);
+    void add(double x, double y);
+    size_t n(void) const;
+    /// \return one-sided p-value
+    double p_approx(bool right) const;
+    double p_exact(bool right, bool left) const;
+    void reset(void);
+    double score(void) const;
+  private:
+    // return # concordant pairs minus # discordant pairs
+    long int count(void) const;
+    // data always sort wrt first and then second (if first are equal)
+    std::multiset<std::pair<double, double> > data_;
+    // calculated in score(void)
+    boost::scoped_ptr<Count> count_;
+  };
+
+
+  // Kendall class
+  Kendall::Kendall(void)
+    : pimpl_(new Pimpl)
+  {
+  }
+
+
+  Kendall::Kendall(const Kendall& rhs)
+    : pimpl_(new Pimpl(*rhs.pimpl_))
+  {
+  }
+
+
+  Kendall::~Kendall(void)
+  {
+    delete pimpl_;
+  }
+
+
+  void Kendall::add(double x, double y)
+  {
+    pimpl_->add(x, y);
+  }
+
+
+  size_t Kendall::n(void) const
+  {
+    return pimpl_->n();
+  }
+
+
+  double Kendall::score(void) const
+  {
+    return pimpl_->score();
+  }
+
+
+  double Kendall::p_left(bool exact) const
+  {
+    if (!exact)
+      return pimpl_->p_approx(false);
+    return pimpl_->p_exact(false, true);
+  }
+
+
+  double Kendall::p_right(bool exact) const
+  {
+    if (!exact)
+      return pimpl_->p_approx(true);
+    return pimpl_->p_exact(true, false);
+  }
+
+
+  double Kendall::p_value(bool exact) const
+  {
+    if (exact)
+      return pimpl_->p_exact(true, true);
+    if (score()>0.0)
+      return 2*p_right(false);
+    return 2*p_left(false);
+  }
+
+
+  void Kendall::reset(void)
+  {
+    pimpl_->reset();
+  }
+
+
+  Kendall& Kendall::operator=(const Kendall& rhs)
+  {
+    if (&rhs == this)
+      return *this;
+
+    assert(pimpl_);
+    assert(rhs.pimpl_);
+    *pimpl_ = *rhs.pimpl_;
+    return *this;
+  }
+
+
+  Kendall::Pimpl::Count::Count(const std::multiset<std::pair<double,double>>& data)
+  {
+    // We follow 3 Algorithm SDTau for some-duplicate datasets in
+    // 'Fast Algorithms For The Calculation Of Kendall's Tau'
+    // by David Christen (Computational Statistics, March 2005)
+
+    // data is sorted w.r.t. ::first
+    calculate_ties(utility::pair_first_iterator(data.begin()),
+                   utility::pair_first_iterator(data.end()),
+                   x_ties_);
+    unsigned int d = 0;
+    unsigned int e = 0;
+    utility::Ranking<double> Y;
+
+    for (auto it=data.cbegin(); it!=data.end(); ++it) {
+      if (it != data.begin()) {
+        auto previous = std::prev(it);
+        if (it->first != previous->first) { // x not equal
+          d = 0;
+          e = 1;
+        }
+        else if (it->second == previous->second) // y also equal
+          ++e;
+        else { // x equal, y not equal
+          d += e;
+          e = 1;
+        }
+      }
+
+      auto lower = Y.lower_bound(it->second);
+      // number of element in Y smaller than it->second
+      int n_smaller = Y.ranking(lower);
+      // number of element in Y equal to it->second
+      int n_equal = 1;
+      auto upper = std::next(lower);
+      while (upper!=Y.cend() && *upper==*lower) {
+        ++upper;
+        ++n_equal;
+      }
+      Y.insert(lower, it->second);
+      size_t i = Y.size();
+
+      long int a = n_smaller - d;
+      long int b = n_equal - e;
+      long int c = i - (a+b+d+e-1);
+
+      extraY_ += d;
+      extraX_ += b;
+      concordant_ += a;
+      disconcordant_ += c;
     }
-
-
-    size_t n(void) const
-    {
-      assert(y_.size()==x_.size());
-      return x_.size();
+    assert(0);
+  }
+
+
+  long int Kendall::Pimpl::Count::count(void) const
+  {
+    return concordant_ - disconcordant_;
+  }
+
+
+  double Kendall::Pimpl::Count::score(void) const
+  {
+    double numerator = count();
+    double denominator = concordant_ + disconcordant_;
+    if (extraX_ || extraY_) {
+      denominator =
+        std::sqrt((denominator + extraX_)*(denominator + extraY_));
     }
-
-
-    /// \return one-sided p-value
-    double p_approx(bool right) const
-    {
-      double k = count();
-      if (!right)
-        k = -k;
-      // avoid overflow
-      double n = this->n();
-      /*
-        According to wikipedia,
-        z = k / sqrt(v)
-        is approximately standard normal
-        v = (v0 - vt - vu)/18 + v1 + v2
-        v0 = n(n-1)(2n+5)
-        vt = \sum t(t-1)(2t+5)
-        vu = \sum u(u-1)(2u+5)
-        v1 = \sum t(t-1)) * \sum u(u-1) / (2n(n-1))
-        v2 = sum t(t-1)(t-2) \sum u(u-1)(u-2) / (9n(n-1)(n-2))
-
-        where t is number of equal values in group i
-        and similarly u for y.
-       */
-      double v0 = n*(n-1)*(2*n+5);
-      double vt = 0;
-      double vu = 0;
-      double v1 = 0;
-      double v2 = 0;
-      // all correction terms above are zero in absence of ties
-      bool x_have_ties = x_.size() != x_count_.size();
-      bool y_have_ties = y_.size() != y_count_.size();
-      if (x_have_ties || y_have_ties) {
-        if (x_have_ties)
-          vt = v_correction(x_count_);
-        if (y_have_ties)
-          vu = v_correction(y_count_);
-        if (x_have_ties && y_have_ties) {
-          v1 = count_pair(x_count_) * (count_pair(y_count_) / (2*n*(n-1)));
-          v2 = count_triple(x_count_);
+    return numerator / denominator;
+  }
+
+
+  double Kendall::Pimpl::Count::variance(void) const
+  {
+    /*
+      According to wikipedia,
+      z = k / sqrt(v)
+      is approximately standard normal
+      v = (v0 - vt - vu)/18 + v1 + v2
+      v0 = n(n-1)(2n+5)
+      vt = \sum t(t-1)(2t+5)
+      vu = \sum u(u-1)(2u+5)
+      v1 = \sum t(t-1)) * \sum u(u-1) / (2n(n-1))
+      v2 = sum t(t-1)(t-2) \sum u(u-1)(u-2) / (9n(n-1)(n-2))
+
+      where t is number of equal values in group i and similarly u for
+      y.
+    */
+    double n = score();
+    double v0 = n*(n-1)*(2*n+5);
+    double vt = 0;
+    double vu = 0;
+    double v1 = 0;
+    double v2 = 0;
+    // all correction terms above are zero in absence of ties
+    bool x_have_ties = x_ties_.have_ties();
+    bool y_have_ties = y_ties_.have_ties();
+    if (x_have_ties || y_have_ties) {
+      if (x_have_ties)
+        vt = x_ties_.v_correction();
+      if (y_have_ties) {
+        vu = y_ties_.v_correction();
+        if (x_have_ties) {
+          v1 = x_ties_.n_pairs() * (y_ties_.n_pairs() / (2*n*(n-1)));
+          v2 = x_ties_.n_triples();
           if (v2)
-            v2 *= count_triple(y_count_) / (9*n*(n-1)*(n-2));
+            v2 *= y_ties_.n_triples() / (9*n*(n-1)*(n-2));
         }
       }
-      double v = (v0 - vt - vu)/18 + v1 + v2;
-      return gsl_cdf_gaussian_Q(k, std::sqrt(v));
     }
-
-
-    double p_exact(bool right, bool left) const
-    {
+    return (v0 - vt - vu)/18 + v1 + v2;
+  }
+
+
+  Kendall::Pimpl::Count::Ties::Ties(void)
+    : n_pairs_(0), n_triples_(0), v_correction_(0)
+  {}
+
+
+  void Kendall::Pimpl::Count::Ties::add(size_t n)
+  {
+    unsigned long int factor = n * (n-1);
+    n_pairs_ += factor;
+    n_triples_ += factor * (n-2);
+    v_correction_ += factor * (2*n+5);
+  }
+
+
+  Kendall::Pimpl::Pimpl(void)
+  {}
+
+
+  Kendall::Pimpl::Pimpl(const Pimpl& other)
+    : data_(other.data_)
+  {}
+
+
+  Kendall::Pimpl& Kendall::Pimpl::operator=(const Pimpl& rhs)
+  {
+    data_ = rhs.data_;
+    count_.reset();
+    return *this;
+  }
+
+
+  void Kendall::Pimpl::add(double x, double y)
+  {
+    data_.insert(std::make_pair(x, y));
+    count_.reset();
+  }
+
+
+  size_t Kendall::Pimpl::n(void) const
+  {
+    return data_.size();
+  }
+
+
+  double Kendall::Pimpl::p_approx(bool right) const
+  {
+    double k = count();
+    if (!right)
+      k = -k;
+    return gsl_cdf_gaussian_Q(k, std::sqrt(count_->variance()));
+  }
+
+
+  double Kendall::Pimpl::p_exact(bool right, bool left) const
+  {
+    assert(0);
+    return 0.0;
+    /*
       long int upper = 0;
       long int lower = 0;
-      if (right && left) {
-        upper = std::max(count(), -count());
-        lower = std::min(count(), -count());
-      }
-      else if (right && !left) {
-        upper = count();
-        lower = std::numeric_limits<long int>::min();
-      }
-      else if (!right && left) {
-        upper = std::numeric_limits<long int>::max();
-        lower = count();
+      if (right) {
+      if (left) {
+      upper = std::max(count(), -count());
+      lower = std::min(count(), -count());
       }
       else {
-        assert(false && "should not end up here");
-      }
-
+      upper = count();
+      lower = std::numeric_limits<long int>::min();
+      }
+      }
+      else {
+      assert(left && "left or right must be true");
+      upper = std::numeric_limits<long int>::max();
+      lower = count();
+      }
+    */
+    /*
       std::vector<double> x(x_);
       std::sort(x.begin(), x.end());
@@ -166 +431 @@
       unsigned int total = 0;
       do {
-        long int k = statistics::count(x.begin(), x.end(), y_.begin());
-        if (k>=upper || k<=lower)
-          ++n;
-        ++total;
+      long int k = statistics::count(x.begin(), x.end(), y_.begin());
+      if (k>=upper || k<=lower)
+      ++n;
+      ++total;
       }
       while (std::next_permutation(x.begin(), x.end()));
       return static_cast<double>(n)/static_cast<double>(total);
-    }
-
-
-    void reset(void)
-    {
-      Pimpl tmp;
-      *this = tmp;
-    }
-
-
-    double score(void) const
-    {
-      double denominator = 0.5*n()*(n()-1.0);
-      if (have_ties()) {
-        double pairs = denominator;
-        denominator = std::sqrt((pairs - 0.5*count_pair(x_count_))*
-                                (pairs - 0.5*count_pair(y_count_)));
-      }
-      return count() / denominator;
-    }
-
-  private:
-    long int count(void) const
-    {
-      if (!updated_) {
-        count_ = statistics::count(x_.begin(), x_.end(), y_.begin());
-        updated_ = true;
-      }
-      return count_;
-    }
-
-    bool have_ties(void) const
-    {
-      return (x_.size() != x_count_.size()) || (y_.size() != y_count_.size());
-    }
-
-
-    /*
-      \return \sum x * (x-1) * (2*x+5) where x is second in map \a count
-     */
-    unsigned long int v_correction(const std::map<double, size_t>& count) const
-    {
-      unsigned long int result = 0;
-      std::map<double, size_t>::const_iterator end=count.end();
-      for (std::map<double, size_t>::const_iterator i=count.begin();i!=end;++i)
-        result += i->second * (i->second - 1) * (2*i->second+5);
-      return result;
-    }
-
-
-    /*
-      \return \sum x * (x-1) where x is second in map \a count
-     */
-    unsigned long int count_pair(const std::map<double, size_t>& count) const
-    {
-      unsigned long int result = 0;
-      typedef std::map<double, size_t>::const_iterator iterator;
-      iterator end=count.end();
-      for (iterator i=count.begin(); i!=end;++i)
-        result += i->second * (i->second - 1);
-      return result;
-    }
-
-    /*
-      \return \sum x * (x-1) * (x-2) where x is second in map \a count
-     */
-    unsigned long int count_triple(const std::map<double, size_t>& count) const
-    {
-      unsigned long int result = 0;
-      std::map<double, size_t>::const_iterator end=count.end();
-      for (std::map<double, size_t>::const_iterator i=count.begin();i!=end;++i)
-        result += i->second * (i->second - 1) * (i->second-2);
-      return result;
-    }
-
-    std::vector<double> x_;
-    std::vector<double> y_;
-    std::map<double, size_t> x_count_;
-    std::map<double, size_t> y_count_;
-    mutable long int count_;
-    mutable bool updated_;
-  };
-
-
-  // Kendall class
-  Kendall::Kendall(void)
-    : pimpl_(new Pimpl)
-  {
-  }
-
-
-  Kendall::Kendall(const Kendall& rhs)
-    : pimpl_(new Pimpl(*rhs.pimpl_))
-  {
-  }
-
-
-  Kendall::~Kendall(void)
-  {
-    delete pimpl_;
-  }
-
-
-  void Kendall::add(double x, double y)
-  {
-    pimpl_->add(x, y);
-  }
-
-
-  size_t Kendall::n(void) const
-  {
-    return pimpl_->n();
-  }
-
-
-  double Kendall::score(void) const
-  {
-    return pimpl_->score();
-  }
-
-
-  double Kendall::p_left(bool exact) const
-  {
-    if (!exact)
-      return pimpl_->p_approx(false);
-    return pimpl_->p_exact(false, true);
-  }
-
-
-  double Kendall::p_right(bool exact) const
-  {
-    if (!exact)
-      return pimpl_->p_approx(true);
-    return pimpl_->p_exact(true, false);
-  }
-
-
-  double Kendall::p_value(bool exact) const
-  {
-    if (exact)
-      return pimpl_->p_exact(true, true);
-    if (score()>0.0)
-      return 2*p_right(false);
-    return 2*p_left(false);
-  }
-
-
-  void Kendall::reset(void)
-  {
-    pimpl_->reset();
-  }
-
-
-  Kendall& Kendall::operator=(const Kendall& rhs)
-  {
-    if (&rhs == this)
-      return *this;
-
-    assert(pimpl_);
-    assert(rhs.pimpl_);
-    *pimpl_ = *rhs.pimpl_;
-    return *this;
+    */
+    return 0;
+  }
+
+
+  void Kendall::Pimpl::reset(void)
+  {
+    Pimpl tmp;
+    *this = tmp;
+  }
+
+
+  double Kendall::Pimpl::score(void) const
+  {
+    count();
+    assert(count_.get());
+    return count_->score();
+  }
+
+
+  long int Kendall::Pimpl::count(void) const
+  {
+    if (!count_)
+      // const_cast to allow lazy eval is more restrictive than
+      // making count_ mutable.
+      const_cast<Pimpl*>(this)->count_.reset(new Count(data_));
+    return count_->count();
   }
 

branches/kendall-score/yat/utility/Makefile.am

@@ -2 +2 @@
 
 # Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
-# Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2020 Peter Johansson
+# Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017 Peter Johansson
 #
 # This file is part of the yat library, http://dev.thep.lu.se/yat
@@ -44 +44 @@
   yat/utility/OptionSwitch.cc \
   yat/utility/PCA.cc \
+  yat/utility/Ranking.cc \
   yat/utility/split.cc \
   yat/utility/stl_utility.cc \
@@ -121 +122 @@
   $(srcdir)/yat/utility/StreamRedirect.h \
   $(srcdir)/yat/utility/Range.h \
+  $(srcdir)/yat/utility/Ranking.h \
   $(srcdir)/yat/utility/stl_utility.h \
   $(srcdir)/yat/utility/StrideIterator.h \