qQuantileNormalizer.h

#ifndef _theplu_yat_normalizer_qquantile_normalizer_
#define _theplu_yat_normalizer_qquantile_normalizer_

/*
  Copyright (C) 2009 Jari Häkkinen, Peter Johansson
  Copyright (C) 2010 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 "utility.h"

#include "yat/regression/CSplineInterpolation.h"
#include "yat/utility/concept_check.h"
#include "yat/utility/DataIterator.h"
#include "yat/utility/DataWeight.h"
#include "yat/utility/Exception.h"
#include "yat/utility/iterator_traits.h"
#include "yat/utility/sort_index.h"
#include "yat/utility/Vector.h"
#include "yat/utility/WeightIterator.h"
#include "yat/utility/yat_assert.h"

#include <boost/concept_check.hpp>

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

namespace theplu {
namespace yat {
namespace utility {
  class VectorBase;
}
namespace normalizer {

  class qQuantileNormalizer
  {
  public:
    template<typename ForwardIterator>
    qQuantileNormalizer(ForwardIterator first, ForwardIterator last,
                        unsigned int Q);

    template<typename RandomAccessIterator1, typename RandomAccessIterator2>
    void operator()(RandomAccessIterator1 first, RandomAccessIterator1 last,
                    RandomAccessIterator2 result) const;

  private:

  class Partitioner
  {
  public:
    template<typename ForwardIterator>
    Partitioner(ForwardIterator first, ForwardIterator last, 
                unsigned int N);

    const utility::Vector& averages(void) const;

    const utility::Vector& quantiles(void) const;

    size_t size(void) const;

  private:
    // unweighted "constructor"
    template<typename ForwardIterator>
    void build(ForwardIterator first, ForwardIterator last, unsigned int N, 
               utility::unweighted_iterator_tag);
    // weighted "constructor"
    template<typename ForwardIterator>
    void build(ForwardIterator first, ForwardIterator last, unsigned int N, 
               utility::weighted_iterator_tag);
    void init(const utility::VectorBase&, unsigned int N);
    void init(const std::vector<utility::DataWeight>&, unsigned int N);

    utility::Vector average_;
    utility::Vector quantiles_;
  };

    Partitioner target_;

    // unweighted version
    template<typename RandomAccessIterator1, typename RandomAccessIterator2>
    void normalize(const Partitioner& source,RandomAccessIterator1 first, 
                   RandomAccessIterator1 last, RandomAccessIterator2 result,
                   utility::unweighted_iterator_tag tag) const;

    // weighted version
    template<typename RandomAccessIterator1, typename RandomAccessIterator2>
    void normalize(const Partitioner& source,RandomAccessIterator1 first, 
                   RandomAccessIterator1 last, RandomAccessIterator2 result,
                   utility::weighted_iterator_tag tag) const;
  };


  // template implementations

  template<typename ForwardIterator>
  qQuantileNormalizer::qQuantileNormalizer(ForwardIterator first, 
                                           ForwardIterator last,
                                           unsigned int Q)
    : target_(Partitioner(first, last, Q))
  {
    YAT_ASSERT(Q>2);
  }


  template<typename RandomAccessIterator1, typename RandomAccessIterator2>
  void qQuantileNormalizer::operator()(RandomAccessIterator1 first, 
                                       RandomAccessIterator1 last,
                                       RandomAccessIterator2 result) const
  {
    Partitioner source(first, last, target_.size());
    typename utility::weighted_if_any2<RandomAccessIterator1, RandomAccessIterator2>::type tag;
    normalize(source, first, last, result, tag);
  }


  template<typename RandomAccessIterator1, typename RandomAccessIterator2>
  void 
  qQuantileNormalizer::normalize(const qQuantileNormalizer::Partitioner& source,
                                 RandomAccessIterator1 first, 
                                 RandomAccessIterator1 last, 
                                 RandomAccessIterator2 result,
                                 utility::unweighted_iterator_tag tag) const
  {
    BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<RandomAccessIterator1>));
    BOOST_CONCEPT_ASSERT((boost::Mutable_RandomAccessIterator<RandomAccessIterator2>));
    BOOST_CONCEPT_ASSERT((utility::DataIteratorConcept<RandomAccessIterator1>));
    BOOST_CONCEPT_ASSERT((utility::DataIteratorConcept<RandomAccessIterator2>));
    utility::check_iterator_is_unweighted(first);
    utility::check_iterator_is_unweighted(result);
    size_t N = last-first;
    YAT_ASSERT(N >= target_.size());

    std::vector<size_t> sorted_index(last-first);
    utility::sort_index(first, last, sorted_index);

    utility::Vector diff(source.averages());
    diff-=target_.averages();
    const utility::Vector& idx=target_.quantiles();
    regression::CSplineInterpolation cspline(idx,diff);

    // linear extrapolation for first part, i.e., use first diff for
    // all points in the first part.
    size_t start=0;
    size_t end = static_cast<size_t>(std::ceil(N*idx(0) - 0.5));
    // take care of limiting case number of parts approximately equal
    // to the number of elements in range.
    if (end==0)
      ++end;
    for (size_t i=start; i<end; ++i) {
      size_t si = sorted_index[i];
      result[si] = first[si] - diff(0);
    }

    using utility::yat_assert;
    
    // cspline interpolation for all data between the mid points of
    // the first and last part
    start=end;
    end = static_cast<size_t>(std::ceil(N*idx(idx.size()-1) - 0.5));
    if (end>=N)
      end = N-1;
    for ( size_t i=start; i<end; ++i) {
      size_t si = sorted_index[i];
      
      YAT_ASSERT((i+0.5)/N>idx(0));
      result[si] = first[si] - cspline.evaluate((i+0.5)/N);
    }
    
    // linear extrapolation for last part, i.e., use last diff for
    // all points in the last part.
    for (size_t i=end ; i<N; ++i) {
      size_t si = sorted_index[i];
      result[si] = first[si] - diff(diff.size()-1);
    }
  }


  template<typename RandomAccessIterator1, typename RandomAccessIterator2>
  void qQuantileNormalizer::normalize(const Partitioner& source,
                                      RandomAccessIterator1 first, 
                                      RandomAccessIterator1 last, 
                                      RandomAccessIterator2 result,
                                      utility::weighted_iterator_tag tag) const
  {
    BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<RandomAccessIterator1>));
    BOOST_CONCEPT_ASSERT((boost::Mutable_RandomAccessIterator<RandomAccessIterator2>));
    BOOST_CONCEPT_ASSERT((utility::DataIterator<RandomAccessIterator1>));
    BOOST_CONCEPT_ASSERT((utility::DataIterator<RandomAccessIterator2>));
    // copy the weights
    detail::copy_weight_if_weighted(first, last, result);

    double total_w = std::accumulate(utility::weight_iterator(first),
                                     utility::weight_iterator(last), 0.0);

    std::vector<size_t> sorted_index(last-first);
    // Code to avoid problems with NaN (ticket:535)
    // utility::sort_index(utility::data_iterator(first),
    //                     utility::data_iterator(last), sorted_index);
    // ... above statement replaced below code block
    {
      using namespace utility;
      std::vector<double> vec;
      vec.reserve(std::distance(first, last));
      std::copy(utility::data_iterator(first), utility::data_iterator(last), 
                std::back_inserter(vec));
      for (std::vector<double>::iterator i(vec.begin()); i!=vec.end(); ++i)
        if (std::isnan(*i))
          *i = std::numeric_limits<double>::infinity();
      utility::sort_index(vec.begin(), vec.end(), sorted_index);
    }
    // end Code to avoid problems with NaN (ticket:535)

    utility::Vector diff(source.averages());
    diff-=target_.averages();
    const utility::Vector& idx=target_.quantiles();
    regression::CSplineInterpolation cspline(idx,diff);

    double sum_w=0;
    utility::iterator_traits<RandomAccessIterator1> traits1;
    utility::iterator_traits<RandomAccessIterator2> traits2;
    for (size_t i=0; i<sorted_index.size(); ++i) {
      size_t si = sorted_index[i];
      double w = (sum_w + 0.5*traits1.weight(first+si))/total_w;
      double correction = 0;
      if (w <= idx(0)) {
        // linear extrapolation for first part, i.e., use first diff for
        // all points in the first part.
        correction = diff(0);
      }
      else if (w < idx(idx.size()-1) ) {
        // cspline interpolation for all data between the mid points of
        // the first and last part
        correction = cspline.evaluate(w);
      }
      else {
        // linear extrapolation for last part, i.e., use last diff for
        // all points in the last part.
        correction = diff(diff.size()-1);
      }
      traits2.data(result+si) = traits1.data(first+si) - correction;
      sum_w += traits1.weight(first+si);
    }
  }


  template<typename ForwardIterator>
  qQuantileNormalizer::Partitioner::Partitioner(ForwardIterator first, 
                                                ForwardIterator last, 
                                                unsigned int N)
    : average_(utility::Vector(N)), quantiles_(utility::Vector(N))
  {
    BOOST_CONCEPT_ASSERT((boost::ForwardIterator<ForwardIterator>));
    BOOST_CONCEPT_ASSERT((utility::DataIterator<ForwardIterator>));
    build(first, last, N, 
          typename utility::weighted_iterator_traits<ForwardIterator>::type());
  }


  template<typename ForwardIterator>
  void qQuantileNormalizer::Partitioner::build(ForwardIterator first, 
                                               ForwardIterator last, 
                                               unsigned int N, 
                                               utility::unweighted_iterator_tag)
  {
    utility::Vector vec(std::distance(first, last));
    std::copy(first, last, vec.begin());
    std::sort(vec.begin(), vec.end());
    init(vec, N);
  }


  template<typename ForwardIterator>
  void qQuantileNormalizer::Partitioner::build(ForwardIterator first, 
                                               ForwardIterator last, 
                                               unsigned int N, 
                                               utility::weighted_iterator_tag)
  {
    using namespace utility;
    std::vector<DataWeight> vec;
    vec.reserve(std::distance(first, last));
    std::back_insert_iterator<std::vector<DataWeight> > inserter(vec);
    std::copy(first, last, inserter);
    std::sort(vec.begin(), vec.end());
    init(vec, N);
  }


}}} // end of namespace normalizer, yat and thep

#endif