source: trunk/yat/statistics/Percentiler.h @ 2076

#ifndef _theplu_yat_statistics_percentiler_ 
#define _theplu_yat_statistics_percentiler_ 

// $Id: Percentiler.h 2076 2009-10-06 22:17:22Z peter $

/*
  Copyright (C) 2008 Jari Häkkinen, Peter Johansson
  Copyright (C) 2009 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/DataWeight.h"
#include "yat/utility/iterator_traits.h"
#include "yat/utility/yat_assert.h"
#include "yat/utility/WeightIterator.h"

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

namespace theplu {
namespace yat {
namespace statistics {  

  /**
     \brief Functor to calculate percentile of a range

     \since New in yat 0.5
   */
  class Percentiler
  {
  public:
    /**
       \param perc percentile to calculate [0,100]. Default value is
       50, which implies class will calculate median.
       \param sorted if true class assumes that ranges are already
       sorted, if false the range will copied to a new range which is
       sorted.
     */
    Percentiler(double perc=50, bool sorted=false);

    /**
       Function is a non-mutable function, i.e., \a first and \a last
       can be const_iterators.

       The \a Nth percentile is defined such that, for example, when
       having four numbers \f$ 0.69 < 1.41 < 3.14 < 28 \f$, the \a Nth
       percentile is:

       - \f$ 0.69 \f$ if \f$ N < 25 \f$
       - \f$ (0.69+1.41)/2 \f$ if \f$ N=25 \f$
       - \f$ 1.41 \f$ if \f$ 25 < N < 50 \f$
       - \f$ (1.41+3.14)/2 \f$ if \f$ N=50 \f$
       - \f$ 3.14 \f$ if \f$ 50 < N < 75 \f$
       - \f$ (3.14+28)/2 \f$ if \f$ N=75 \f$
       - \f$ 28 \f$ if \f$ 75 < N \f$

       Similarily, if we have a weighted range \f$ x_0=0.69, w_0=1 ;
       x_1=1.41, w_1=0 ; x_2=3.14, w_2=0.5 ; x_3=28, w_3=0.5 \f$, we
       calculate the accumulated normalized weight \f$ S_k = \frac
       {100}{\sum w_i}\sum_{i=0}^k w_i \f$ and the percentile is

       - \f$ 0.69 \f$ if \f$ N < S_0 \f$
       - \f$ (0.69+3.14)/2 \f$ if \f$ N=S_0  \f$
       - \f$ 3.14 \f$ if \f$ S_0 < N < S_2 \f$
       - \f$ (3.14+28)/2 \f$ if \f$ N=S_2 \f$
       - \f$ 28 \f$ if \f$ S_2 < N \f$
       
       Note, that data point with weight zero is completely ignored. 

       \return percentile of range
     */
    template<typename RandomAccessIterator>
    double operator()(RandomAccessIterator first, 
                      RandomAccessIterator last) const
    {
      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 <= w_bound+margin && upper!=accum_w.end())
        ++upper;
      if (upper==accum_w.end())
        --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