Percentiler.h

#ifndef _theplu_yat_statistics_percentiler_ 
#define _theplu_yat_statistics_percentiler_ 

// $Id: Percentiler.h 3550 2017-01-03 05:41:02Z peter $

/*
  Copyright (C) 2008, 2009 Jari Häkkinen, Peter Johansson
  Copyright (C) 2010, 2011, 2014, 2016 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include "yat/utility/concept_check.h"
#include "yat/utility/DataWeight.h"
#include "yat/utility/iterator_traits.h"
#include "yat/utility/yat_assert.h"
#include "yat/utility/WeightIterator.h"

#include <boost/concept_check.hpp>
#include <boost/iterator/iterator_concepts.hpp>

#include <algorithm>
#include <cmath>
#include <limits>
#include <numeric>
#include <stdexcept>
#include <vector>

namespace theplu {
namespace yat {
namespace statistics {

  class Percentiler
  {
  public:
    Percentiler(double perc=50, bool sorted=false);

    template<typename RandomAccessIterator>
    double operator()(RandomAccessIterator first,
                      RandomAccessIterator last) const
    {
      BOOST_CONCEPT_ASSERT((utility::DataIteratorConcept<RandomAccessIterator>));
      BOOST_CONCEPT_ASSERT((boost_concepts::RandomAccessTraversal<RandomAccessIterator>));
      if (first==last)
        return std::numeric_limits<double>::quiet_NaN();
      // just to avoid long line
      using utility::weighted_iterator_traits;
      typedef typename weighted_iterator_traits<RandomAccessIterator>::type tag;
      return calculate(first, last, sorted_, tag());
    }
  private:
    double perc_;
    bool sorted_;

    // unweighted version
    template<typename RandomAccessIterator>
    double calculate(RandomAccessIterator first, RandomAccessIterator last,
                     bool sorted, utility::unweighted_iterator_tag tag) const;

    // weighted version
    template<typename RandomAccessIterator>
    double calculate(RandomAccessIterator first, RandomAccessIterator last,
                     bool sorted, utility::weighted_iterator_tag tag) const;

    // using compiler generated copy
    //Percentiler(const Percentiler&);
    //Percentiler& operator=(const Percentiler&);

  };

  // template implementations
  // /////////////////////////

  // unweighted version
  template<typename RandomAccessIterator>
  double
  Percentiler::calculate(RandomAccessIterator first,
                         RandomAccessIterator last,
                         bool sorted,
                         utility::unweighted_iterator_tag tag) const
  {
    size_t n = last - first;
    // range is one value only is a special case
    if (n == 1)
      *first;

    double j = n * perc_ / 100.0;

    // Some special cases
    if (j > n-1) {
      if (sorted)
        return *(--last);
      return *std::max_element(first, last);
    }
    if (j < 1.0) {
      if (sorted)
        return *first;
      return *std::min_element(first, last);
    }

    size_t i = static_cast<size_t>(j);
    // for border case (j integer) we take the average of adjacent bins
    size_t k = (i==j) ? i-1 : i;

    if (sorted) {
      if (i==k)
        return first[i];
      return (first[i]+first[k])/2;
    }

    std::vector<double> vec(first, last);
    // find the ith element
    std::nth_element(vec.begin(), vec.begin()+i, vec.end());
    // if i==k simply return vec[i]
    if (i==k)
      return vec[i];
    YAT_ASSERT(k==i-1);
    // otherwise return average of vec[i] and vec[k].

    // We need to find the kth element. Since we have called
    // nth_element above, we know that all elements in (begin+i, end)
    // are >= vec[i] and all elements in [begin, begin+i) are <=
    // vec[i]. Consequently, kth (k=i-1) element is the max element in
    // range [begin, begin+i).
    return (vec[i] + *std::max_element(vec.begin(), vec.begin()+i))/2;
  }


  // weighted version
  template<typename RandomAccessIterator>
  double Percentiler::calculate(RandomAccessIterator first,
                                RandomAccessIterator last,
                                bool sorted,
                                utility::weighted_iterator_tag tag) const
  {
    if (sorted) {
      utility::iterator_traits<RandomAccessIterator> trait;
      std::vector<double> accum_w;
      accum_w.reserve(last-first);
      std::partial_sum(weight_iterator(first),
                       weight_iterator(last),
                       std::back_inserter(accum_w));

      double w_bound=perc_/100.0*accum_w.back();
      std::vector<double>::const_iterator upper(accum_w.begin());
      double margin=1e-10;
      while (upper!=accum_w.end() && *upper <= w_bound+margin)
        ++upper;
      while (upper!=accum_w.begin() &&
             (upper==accum_w.end() ||
              trait.weight(first+(upper-accum_w.begin()))==0.0))
        --upper;
      std::vector<double>::const_iterator lower(upper);
      while ( *(lower-1)>=w_bound-margin && lower>accum_w.begin())
        --lower;

      return (trait.data(first+(upper-accum_w.begin()))+
              trait.data(first+(lower-accum_w.begin())))/2;
    }

    std::vector<utility::DataWeight> v_copy(first, last);
    std::sort(v_copy.begin(), v_copy.end());
    return calculate(v_copy.begin(), v_copy.end(), true, tag);
  }


}}} // of namespace statistics, yat, and theplu

#endif