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

#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/Vector.h"
#include "yat/utility/WeightIterator.h"
#include "yat/utility/yat_assert.h"

#include <algorithm>
#include <iterator>
#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.

     There is currently no weighted version of qQuantileNormalizer

     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 BidirectionalIterator>
    qQuantileNormalizer(BidirectionalIterator first, BidirectionalIterator last,
                        unsigned int Q);

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

       It is possible to normalize "in place"; it is permissible for
       \a matrix and \a result to reference the same Matrix.

       \note dimensions of \a matrix and \a result must match.
     */
    template<typename RandomAccessIterator1, typename RandomAccessIterator2>
    void operator()(RandomAccessIterator1 first, RandomAccessIterator1 last,
                    RandomAccessIterator2 result) const;

  private:

  /**
     \brief Partition a vector 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 BidirectionalIterator>
    Partitioner(BidirectionalIterator first, BidirectionalIterator 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 Iterator>
    void build(Iterator first, Iterator last, unsigned int N, 
               utility::unweighted_iterator_tag);
    // weighted "constructor"
    template<typename Iterator>
    void build(Iterator first, Iterator 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
    // copy weights and apply unweighted version on data part
    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 BidirectionalIterator>
  qQuantileNormalizer::qQuantileNormalizer(BidirectionalIterator first, 
                                           BidirectionalIterator 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
  {
    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 interpolation for first part, i.e., use first diff for
    // all points in the first part.
    size_t start=0;
    size_t end=static_cast<unsigned int>(idx(0));
    // The first condition below takes care of limiting case number
    // of parts approximately equal to the number of matrix rows and
    // the second condition makes sure that index is large enough
    // when using cspline below ... the static cast above takes the
    // floor whereas we want to take the "roof" forcing next index
    // range to be within interpolation range for the cspline.
    if ((end==0) || (end<idx(0)))
      ++end;
    for (size_t row=start; row<end; ++row) {
      size_t srow=sorted_index[row];
      result[srow] = first[srow] - diff(0);
    }
    
    // cspline interpolation for all data between the mid points of
    // the first and last part
    start=end;
    end=static_cast<unsigned int>(idx(target_.size()-1));
    // take care of limiting case number of parts approximately
    // equal to the number of matrix rows
    if (end==(static_cast<size_t>(N-1)) )
      --end;
    for (size_t row=start; row<=end; ++row) {
      size_t srow=sorted_index[row];
      result[srow] = first[srow] - cspline.evaluate(row) ;
    }
    
    // linear interpolation for last part, i.e., use last diff for
    // all points in the last part.
    start=end+1;
    end=N;
    for (size_t row=start; row<end; ++row) {
      size_t srow=sorted_index[row];
      result[srow] = first[srow] - 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<RandomAccessIterator1>(first), 
              utility::weight_iterator<RandomAccessIterator1>(last), 
              utility::weight_iterator<RandomAccessIterator2>(result)); 
    // apply algorithm on data part of range
    normalize(source, utility::data_iterator<RandomAccessIterator1>(first),
              utility::data_iterator<RandomAccessIterator1>(last),
              utility::data_iterator<RandomAccessIterator2>(result),
              utility::unweighted_iterator_tag());
  }


  template<typename BidirectionalIterator>
  qQuantileNormalizer::Partitioner::Partitioner(BidirectionalIterator first, 
                                                BidirectionalIterator last, 
                                                unsigned int N)
    : average_(utility::Vector(N)), index_(utility::Vector(N))
  {
    typedef typename 
      utility::weighted_iterator_traits<BidirectionalIterator>::type tag;
    build(first, last, N, tag());
          
  }


  template<typename Iterator>
  void qQuantileNormalizer::Partitioner::build(Iterator first, Iterator 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 Iterator>
  void qQuantileNormalizer::Partitioner::build(Iterator first, 
                                               Iterator 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