source: trunk/yat/normalizer/qQuantileNormalizer.h @ 2119

#ifndef _theplu_yat_normalizer_qquantile_normalizer_
#define _theplu_yat_normalizer_qquantile_normalizer_

/*
  Copyright (C) 2009 Jari Häkkinen, 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/regression/CSplineInterpolation.h"
#include "yat/utility/DataIterator.h"
#include "yat/utility/DataWeight.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 <algorithm>
#include <cmath>
#include <iterator>
#include <limits>
#include <numeric>
#include <stdexcept>
#include <vector>

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

  /**
     \brief Perform Q-quantile normalization

     Perform a Q-quantile normalization on a \a source range, after
     which it will approximately have the same distribution of data as
     the \a target range (the Q-quantiles are the same). The rank of
     an element in the \a source range is not changed.

     The class works also with unweighed ranges, and there is no
     restriction that weighted \a source range requires weighted \a
     target range or vice versa.

     Normalization goes like this:
     - Data are not assumed to be sorted.
     - Partition sorted \a target data in Q parts. Q must be 3 or larger
       because of requirements from the underlying cubic spline fit
     - Calculate the arithmetic (weighted) mean for each part, the mean is
       assigned to the mid point of each part.
     - Do the above for the data to be transformed (called \a source
       here).
     - For each part, calculate the difference between the \a target
       and \a the source. Now we have \a Q differences \f$ d_i \f$
       with associated rank (midpoint of each part).
     - Create a cubic spline fit to this difference vector \a d. The
       resulting curve is used to recalculate all values in \a source.
       - Use the cubic spline fit for values within the cubic spline
         fit range [midpoint 1st part, midpoint last part].
       - For data outside the cubic spline fit use linear
         extrapolation, i.e., a constant shift. \f$ d_{first} \f$ for
         points below fit range, and \f$ d_{last} \f$ for points above fit
         range.

     \since New in yat 0.5
   */
  class qQuantileNormalizer
  {
  public:
    /**
       \brief Constructor

       Divides a sorted copy of range [\a first,\a last) into \a Q
       parts. Parts are divided such that the sum of weights is
       approximately the same in the different parts. If a relative
       weight, \f$ w_i / \sum w_i \f$, is larger than 1/Q this might
       be difficult to achieve, in which case a an exception is
       thrown. In the unweighted case this implies that \a Q should be
       smaller (or equal) than number of elements in [\a first, \a
       last).

       The range supplied to the constructor sets the target
       distribution.

       As the \a source range is also divided into \a Q parts (when
       operator() is called), it is recommended to keep \a Q smaller
       (or equal) than the size of ranges that will be normalized.

       Also, \a Q must not be smaller than 3 due to restrictions in
       the cubic spline fit utilized in the normalization.

       \b Type \b Requirements:
       - \c ForwardIterator is a model of \forward_iterator.
       - \c ForwardIterator is either a \ref concept_weighted_iterator
         or unweighted with value type convertible to \c double.

    */
    template<typename ForwardIterator>
    qQuantileNormalizer(ForwardIterator first, ForwardIterator last,
                        unsigned int Q);

    /**
       \brief Perform the Q-quantile normalization.

       Elements in [\a first, \a last) are normalized as described
       above and the result is assigned to [\a result, \a result + \a
       last-\a first). Input range [\a first, \a last) is not
       modified. If ranges are weighted, the weights are copied from
       [\a first, \a last) to \a result range.

       It is possible to normalize "in place"; it is permissible for
       \a first and \a result to be the same. However, as assignment
       occurs sequentially, the operation is undefined if \a result is
       the same as any in range [\a first, \a last).

       \b Type Requirements:
       - \c RandomAccessIterator1 is a model of \random_access_iterator
       - \c RandomAccessIterator1 is either a \ref
         concept_weighted_iterator or unweighted with value type
         convertible to \c double.
       - \c RandomAccessIterator1's value type is convertible to
         \c RandomAccessIterator2's value type
       - \c RandomAccessIterator2 is a model of \random_access_iterator 
       - \c RandomAccessIterator2 is mutable.

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

  private:

  /**
     \brief Partition a range of data into equal sizes.

     Copy the range [first, last), sort the copy, and divide the
     sorted copy in Q parts. The parts are created such that the total
     weight in a part is approximately W/Q where W is the total weight
     (over all parts). The class calculates the average value in each
     part and also the "quantile". 
  */
  class Partitioner
  {
  public:
    /**
       \brief Create the partition and perform required calculations.
    */
    template<typename ForwardIterator>
    Partitioner(ForwardIterator first, ForwardIterator last, 
                unsigned int N);

    /**
       \brief Return the averages for each part.

       \return The average vector.
    */
    const utility::Vector& averages(void) const;

    /**
       The quantile (here) is defined as (w_lower + w_upper) / 2W,
       where w_lower is the total weight of elements smaller than the
       smallest element in the part, and w_upper is the total weight
       of elements smaller (or equal) than the largest value in the
       part.  

       In the unweighted case all weights are 1.0, which implies q_0 =
       n_0/N, q_1 = (n_0+n_1/2)/N, q_2 = (n_0+n_1+n_2/2)/N where n_i
       is number of elements in ith part.

       \return The quantiles vector.
    */
    const utility::Vector& quantiles(void) const;

    /**
       \return The number of parts.
    */
    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))
  {
    utility::yat_assert<std::runtime_error>(Q>2,
                                            "qQuantileNormalizer: Q too small");
  }


  template<typename RandomAccessIterator1, typename RandomAccessIterator2>
  void qQuantileNormalizer::operator()(RandomAccessIterator1 first, 
                                       RandomAccessIterator1 last,
                                       RandomAccessIterator2 result) const
  {
    Partitioner source(first, last, target_.size());
    typename utility::weighted_iterator_traits<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
  {
    utility::check_iterator_is_unweighted(first);
    utility::check_iterator_is_unweighted(result);
    size_t N = last-first;
    utility::yat_assert<std::runtime_error>
      (N >= target_.size(), "qQuantileNormalizer: Input range too small");

    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<std::runtime_error>((i+0.5)/N>idx(0), 
                               "qQuantileNormalizer: invalid input to cspline");
      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
  {
    // copy the weights
    std::copy(utility::weight_iterator(first), 
              utility::weight_iterator(last), 
              utility::weight_iterator(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))
  {
    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