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

#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

     After a Q-quantile normalization each column has approximately
     the same distribution of data (the Q-quantiles are the
     same). Also, within each column the rank of an element is not
     changed.

     The normalization goes like this
     - Data is not assumed to be sorted.
     - Partition sorted target data in N parts. N must be 3 larger
       because of requirements from the underlying cspline fit
     - Calculate the arithmetic 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 source
       here).
     - For each part, calculate the difference between the target and
       the source. Now we have N differences d_i with associated rank
       (midpoint of each part).
     - Create a cubic spline fit to this difference vector d. The
       resulting curve is used to recalculate all column values.
       - 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. d_first for points
         below fit range, and d_last for points above fit range.

     \since New in yat 0.5
   */
  class qQuantileNormalizer
  {
  public:
    /**
       \brief Documentation please.

       \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