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

#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 <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 N parts. N must be 3 or larger
       because of requirements from the underlying cspline 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 tranformed (called \a source
       here).
     - For each part, calculate the difference between the \a target
       and \a the source. Now we have \a N 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 Contructor

       \a Q is the number of parts and must be within \f$ [3,N] \f$
       where \f$ N \f$ is the total number of data points in the
       target. However, if \f$ N \f$ is larger than the number of points
       in the data to be normalized the behaviour of the code is
       undefined. Keep \f$ N \f$ equal to or less than the smallest
       number of data points in the target or each data set to be
       normalized against a given target. The lower bound of three is
       due to restrictions in the cspline fit utilized in the
       normalization.
    */
    template<typename ForwardIterator>
    qQuantileNormalizer(ForwardIterator first, ForwardIterator last,
                        unsigned int Q);

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

       It is possible to normalize "in place"; it is permissible for
       \a first and \a result to be the same.
     */
    template<typename RandomAccessIterator1, typename RandomAccessIterator2>
    void operator()(RandomAccessIterator1 first, RandomAccessIterator1 last,
                    RandomAccessIterator2 result) const;

  private:

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

     The class also calculates the average of each part and assigns
     the average to the mid point of each part. The midpoint is a
     double, i.e., it is not forced to be an integer index.
  */
  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;

    /**
       \brief Return the mid point for each partition.

       \return The index vector.
    */
    const utility::Vector& index(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 index_;
  };

    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_.index();
    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);
    utility::sort_index(utility::data_iterator(first), 
                        utility::data_iterator(last), sorted_index);

    utility::Vector diff(source.averages());
    diff-=target_.averages();
    const utility::Vector& idx=target_.index();
    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)), index_(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)
  {
    std::vector<utility::DataWeight> vec;
    vec.reserve(std::distance(first, last));
    std::back_insert_iterator<std::vector<utility::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