utility.h

#ifndef _theplu_yat_statistics_utility_
#define _theplu_yat_statistics_utility_

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

/*
  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
  Copyright (C) 2005 Peter Johansson
  Copyright (C) 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
  Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
  Copyright (C) 2009, 2010, 2011, 2013, 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 "Percentiler.h"

#include "yat/classifier/DataLookupWeighted1D.h"
#include "yat/classifier/Target.h"
#include "yat/normalizer/utility.h"
#include "yat/utility/concept_check.h"
#include "yat/utility/DataIterator.h"
#include "yat/utility/deprecate.h"
#include "yat/utility/iterator_traits.h"
#include "yat/utility/sort_index.h"
#include "yat/utility/Vector.h"
#include "yat/utility/VectorBase.h"
#include "yat/utility/yat_assert.h"

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

#include <gsl/gsl_math.h>
#include <gsl/gsl_statistics_double.h>

#include <algorithm>
#include <cmath>
#include <iterator>
#include <stdexcept>
#include <vector>

namespace theplu {
namespace yat {
namespace statistics {

  template <typename T, typename ForwardIterator>
  void add(T& o, ForwardIterator first, ForwardIterator last,
           const classifier::Target& target);

  template<typename BidirectionalIterator1, typename BidirectionalIterator2>
  void benjamini_hochberg(BidirectionalIterator1 first,
                          BidirectionalIterator1 last,
                          BidirectionalIterator2 result);


  template<typename RandomAccessIterator, typename MutableRandomAccessIterator>
  void benjamini_hochberg_unsorted(RandomAccessIterator first,
                                   RandomAccessIterator last,
                                   MutableRandomAccessIterator result);


  double cdf_hypergeometric_P(unsigned int k, unsigned int n1,
                              unsigned int n2, unsigned int t);


  template<typename InputIterator>
  double entropy(InputIterator first, InputIterator last);

  double pearson_p_value(double r, unsigned int n);

  double kurtosis(const utility::VectorBase&);


  template <class RandomAccessIterator>
  double mad(RandomAccessIterator first, RandomAccessIterator last,
             bool sorted=false);


  template <class RandomAccessIterator>
  double median(RandomAccessIterator first, RandomAccessIterator last,
                bool sorted=false);


  template<class T>
  double mutual_information(const T& A);

  template <class RandomAccessIterator>
  double percentile(RandomAccessIterator first, RandomAccessIterator last,
                    double p, bool sorted=false) YAT_DEPRECATE;

  template <class RandomAccessIterator>
  double percentile2(RandomAccessIterator first, RandomAccessIterator last,
                     double p, bool sorted=false)
  {
    Percentiler percentiler(p, sorted);
    return percentiler(first, last);
  }

  double skewness(const utility::VectorBase&);



  template <typename T, typename ForwardIterator>
  void add(T& o, ForwardIterator first, ForwardIterator last,
           const classifier::Target& target)
  {
    BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<ForwardIterator>));
    BOOST_CONCEPT_ASSERT((boost_concepts::SinglePassIterator<ForwardIterator>));
    utility::iterator_traits<ForwardIterator> traits;
    for (size_t i=0; first!=last; ++i, ++first)
      o.add(traits.data(first), target.binary(i), traits.weight(first));
  }


  template<typename BidirectionalIterator1, typename BidirectionalIterator2>
  void benjamini_hochberg(BidirectionalIterator1 first,
                          BidirectionalIterator1 last,
                          BidirectionalIterator2 result)
  {
    BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<BidirectionalIterator1>));
    BOOST_CONCEPT_ASSERT((boost_concepts::BidirectionalTraversal<BidirectionalIterator1>));

    BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<BidirectionalIterator2>));
    BOOST_CONCEPT_ASSERT((boost_concepts::WritableIterator<BidirectionalIterator2>));
    BOOST_CONCEPT_ASSERT((boost_concepts::BidirectionalTraversal<BidirectionalIterator2>));

    typedef typename boost::iterator_difference<BidirectionalIterator1>::type
      diff_type;

    diff_type n = std::distance(first, last);
    if (!n)
      return;
    std::advance(result, n-1);
    first = last;
    --first;
    diff_type rank = n;

    double x = 1.0;
    while (rank) {
      x = std::min(*first * n/static_cast<double>(rank), x);
      *result = x;
      --rank;
      --first;
      --result;
    }
  }


  template<typename RandomAccessIterator, typename MutableRandomAccessIterator>
  void benjamini_hochberg_unsorted(RandomAccessIterator first,
                                   RandomAccessIterator last,
                                   MutableRandomAccessIterator result)
  {
    BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<RandomAccessIterator>));
    BOOST_CONCEPT_ASSERT((boost_concepts::RandomAccessTraversal<RandomAccessIterator>));

    BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<MutableRandomAccessIterator>));
    BOOST_CONCEPT_ASSERT((boost_concepts::WritableIterator<MutableRandomAccessIterator>));
    BOOST_CONCEPT_ASSERT((boost_concepts::RandomAccessTraversal<MutableRandomAccessIterator>));

    std::vector<size_t> idx;
    utility::sort_index(first, last, idx);
    benjamini_hochberg(boost::make_permutation_iterator(first, idx.begin()),
                       boost::make_permutation_iterator(first, idx.end()),
                       boost::make_permutation_iterator(result, idx.begin()));
  }


  template<typename InputIterator>
  double entropy(InputIterator first, InputIterator last)
  {
    BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<InputIterator>));
    BOOST_CONCEPT_ASSERT((boost_concepts::SinglePassIterator<InputIterator>));
    using boost::Convertible;
    typedef typename boost::iterator_value<InputIterator>::type T;
    BOOST_CONCEPT_ASSERT((Convertible<T,double>));
    double sum = 0;
    double N = 0;
    for (; first != last; ++first) {
      if (*first) {
        N += *first;
        sum += *first * std::log(static_cast<double>(*first));
      }
    }
    return -sum / N + log(N);
  }


  template <class RandomAccessIterator>
  double mad(RandomAccessIterator first, RandomAccessIterator last,
             bool sorted)
  {
    BOOST_CONCEPT_ASSERT((utility::DataIteratorConcept<RandomAccessIterator>));
    BOOST_CONCEPT_ASSERT((boost_concepts::RandomAccessTraversal<RandomAccessIterator>));
    double m = median(first, last, sorted);
    typedef typename boost::iterator_value<RandomAccessIterator>::type T;
    std::vector<T> ad(std::distance(first, last));
    // copy weights if needed
    normalizer::detail::copy_weight_if_weighted(first, last, ad.begin());
    // assign data (absolute deviation from median)
    utility::iterator_traits<RandomAccessIterator> traits;
    utility::DataIterator<typename std::vector<T>::iterator> first2(ad.begin());
    while (first!=last) {
      *first2 = std::abs(traits.data(first)-m);
      ++first;
      ++first2;
    }
    return median(ad.begin(), ad.end(), false);
  }


  template <class RandomAccessIterator>
  double median(RandomAccessIterator first, RandomAccessIterator last,
                bool sorted)
  {
    BOOST_CONCEPT_ASSERT((utility::DataIteratorConcept<RandomAccessIterator>));
    BOOST_CONCEPT_ASSERT((boost_concepts::RandomAccessTraversal<RandomAccessIterator>));
    return percentile2(first, last, 50.0, sorted);
  }


  template<class T>
  double mutual_information(const T& n)
  {
    BOOST_CONCEPT_ASSERT((utility::Container2D<T>));
    using boost::Convertible;
    BOOST_CONCEPT_ASSERT((Convertible<typename T::value_type,double>));

    // p_x = \sum_y p_xy

    // Mutual Information is defined as
    // \sum_xy p_xy * log (p_xy / (p_x p_y)) =
    // \sum_xy p_xy * [log p_xy - log p_x - log p_y]
    // \sum_xy p_xy log p_xy - p_xy log p_x - p_xy log p_y
    // \sum_xy p_xy log p_xy - \sum_x p_x log p_x - \sum_y p_y log p_y
    // - entropy_xy + entropy_x + entropy_y

    utility::Vector rowsum(n.columns(), 0);
    for (size_t c = 0; c<n.columns(); ++c)
      rowsum(c) = std::accumulate(n.begin_column(c), n.end_column(c), 0);

    utility::Vector colsum(n.rows(), 0);
    for (size_t r = 0; r<n.rows(); ++r)
      colsum(r) = std::accumulate(n.begin_row(r), n.end_row(r), 0);

    double mi  = - entropy(n.begin(), n.end());
    mi += entropy(rowsum.begin(), rowsum.end());
    mi += entropy(colsum.begin(), colsum.end());

    return mi/M_LN2;
  }


  template <class RandomAccessIterator>
  double percentile(RandomAccessIterator first, RandomAccessIterator last,
                    double p, bool sorted)
  {
    BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<RandomAccessIterator>));
    BOOST_CONCEPT_ASSERT((utility::DataIteratorConcept<RandomAccessIterator>));
    // range is one value only is a special case
    if (first+1 == last)
      return utility::iterator_traits<RandomAccessIterator>().data(first);
    if (sorted) {
      // have to take care of this special case
      if (p>=100)
        return utility::iterator_traits<RandomAccessIterator>().data(--last);
      double j = p/100 * (std::distance(first,last)-1);
      int i = static_cast<int>(j);
      return (1-j+floor(j))*first[i] + (j-floor(j))*first[i+1];
    }
    using std::iterator_traits;
    typedef typename iterator_traits<RandomAccessIterator>::value_type value_t;
    std::vector<value_t> v_copy(first, last);
    size_t i = static_cast<size_t>(p/100 * (v_copy.size()-1));
    if (i+2 < v_copy.size()) {
      std::partial_sort(v_copy.begin(), v_copy.begin()+i+2, v_copy.end());
    }
    else
      std::sort(v_copy.begin(), v_copy.end());
    return percentile(v_copy.begin(), v_copy.end(), p, true);
  }


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

#endif