source: trunk/yat/normalizer/qQuantileNormalizer.cc @ 1718

// $Id: qQuantileNormalizer.cc 1718 2009-01-14 15:42:57Z jari $

/*
  Copyright (C) 2009 Jari Häkkinen

  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 "qQuantileNormalizer.h"

#include "yat/regression/CSplineInterpolation.h"
#include "yat/statistics/Averager.h"
#include "yat/utility/Matrix.h"
#include "yat/utility/Vector.h"
#include "yat/utility/VectorBase.h"

#include <algorithm>
#include <cassert>

namespace theplu {
namespace yat {
namespace normalizer {


  qQuantileNormalizer::Partitioner::Partitioner(const utility::VectorBase& vec,
                                                unsigned int N)
    : average_(utility::Vector(N)), index_(utility::Vector(N))
  {
    assert(N>1);
    assert(N<=vec.size());
    double range=static_cast<double>(vec.size())/N;
    assert(range);
    utility::Vector sortedvec(vec);
    std::sort(sortedvec.begin(),sortedvec.end());
    unsigned int start=0;
    for (unsigned int i=0; i<N; ++i) {
      unsigned int end = ( i==(N-1) ? sortedvec.size() :
                           static_cast<unsigned int>((i+1)*range) );
      statistics::Averager av;
      for (unsigned int r=start; r<end; ++r)
        av.add(sortedvec(r));
      average_(i)=av.mean();
      index_(i)= static_cast<double>(end+start)/2;
      start=end;
    }
  }


  const utility::Vector& qQuantileNormalizer::Partitioner::averages(void) const
  {
    return average_;
  }


  const utility::Vector& qQuantileNormalizer::Partitioner::index(void) const
  {
    return index_;
  }


  size_t qQuantileNormalizer::Partitioner::size(void) const
  {
    return average_.size();
  }


  qQuantileNormalizer::qQuantileNormalizer(const utility::VectorBase& target,
                                           unsigned int Q)
    : target_(Partitioner(target,Q))
  {
    assert(Q>2); // required by cspline fit
  }


  void qQuantileNormalizer::operator()(const utility::Matrix& matrix,
                                       utility::Matrix& result) const
  {
    assert(matrix.rows()    == result.rows());
    assert(matrix.columns() == result.columns());
    assert(matrix.rows()    >= target_.size());

    std::vector<std::vector<size_t> > sorted_index(matrix.rows());
    for (size_t column=0; column<matrix.columns(); ++column)
      utility::sort_index(sorted_index[column],
                          matrix.column_const_view(column));

    for (size_t column=0; column<matrix.columns(); ++column) {
      Partitioner source(matrix.column_const_view(column),target_.size());
      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));
      for (size_t row=start; row<end; ++row) {
        size_t srow=sorted_index[column][row];
        result(srow,column) = matrix(srow,column) - diff(0);
      }

      // cspline interpolation for all data between the first and last
      // parts
      start=static_cast<unsigned int>(idx(0));
      end=static_cast<unsigned int>(idx(target_.size()-1));
      for (size_t row=start; row<=end; ++row) {
        size_t srow=sorted_index[column][row];
        result(srow,column) = matrix(srow,column) - cspline.evaluate(row) ;
      }

      // linear interpolation for last part, i.e., use last diff for
      // all points in the last part.
      start=static_cast<unsigned int>(idx(target_.size()-1)+1);
      end=result.rows();
      for (size_t row=start; row<end; ++row) {
        size_t srow=sorted_index[column][row];
        result(srow,column) = matrix(srow,column) - diff(diff.size()-1);
      }
    }
  }

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