Percentiler.h

#ifndef _theplu_yat_statistics_percentiler_ 
#define _theplu_yat_statistics_percentiler_ 

// $Id: Percentiler.h 2470 2011-04-12 03:21:45Z peter $

/*
  Copyright (C) 2008, 2009 Jari Häkkinen, Peter Johansson
  Copyright (C) 2010, 2011 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 <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((boost::RandomAccessIterator<RandomAccessIterator>));
      BOOST_CONCEPT_ASSERT((utility::DataIteratorConcept<RandomAccessIterator>));
      if (first==last)
        return std::numeric_limits<double>::quiet_NaN();
      return calculate(first, last, sorted_, 
       typename utility::weighted_iterator_traits<RandomAccessIterator>::type());
    }
  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
  {
    // range is one value only is a special case
    if (first+1 == last)
      *first;
    if (sorted) {
      size_t n = last - first;
      // have to take care of this special case
      if (perc_>= 100.0)
        return *(--last);
      if (perc_<= 0.0)
        return *first;
      double j = n * perc_ / 100.0;
      size_t i = static_cast<size_t>(j);
      if (i==j)
        return (first[i]+first[i-1])/2;
      return first[i];
    }

    std::vector<double> v_copy;
    v_copy.reserve(std::distance(first,last));
    std::copy(first, last, std::back_inserter(v_copy));
    std::sort(v_copy.begin(), v_copy.end());
    return calculate(v_copy.begin(), v_copy.end(), true, tag);
  }


  // 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;
    v_copy.reserve(last-first);
    std::copy(first, last, std::back_inserter(v_copy));
    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