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

// $Id: qQuantileNormalizer.cc 2103 2009-11-06 12:40:26Z peter $

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

#include "yat/regression/CSplineInterpolation.h"
#include "yat/statistics/Averager.h"
#include "yat/statistics/AveragerWeighted.h"
#include "yat/utility/DataWeight.h"
#include "yat/utility/Vector.h"
#include "yat/utility/VectorBase.h"
#include "yat/utility/WeightIterator.h"

#include <algorithm>
#include <cassert>
#include <numeric>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>

namespace theplu {
namespace yat {
namespace normalizer {


  void 
  qQuantileNormalizer::Partitioner::init(const utility::VectorBase& sortedvec,
                                         unsigned int N)
  {
    assert(N>1);
    assert(N<=sortedvec.size());
    double range=static_cast<double>(sortedvec.size())/N;
    assert(range);
    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();
      quantiles_(i)   = 0.5*(end+start);
      start=end;
    }
    // rescale quantiles to be in range (0,1)
    quantiles_ *= 1.0/sortedvec.size();
  }


  void qQuantileNormalizer::Partitioner::init
  (const std::vector<utility::DataWeight>& sortedvec, unsigned int N)
  {
    assert(N>1);
    assert(N<=sortedvec.size());
    double total_w = std::accumulate(utility::weight_iterator(sortedvec.begin()),
                                     utility::weight_iterator(sortedvec.end()), 
                                     0.0);

    assert(total_w);
    double sum_w = 0;
    std::vector<utility::DataWeight>::const_iterator iter(sortedvec.begin());
    for (unsigned int i=0; i<N; ++i) {
      statistics::AveragerWeighted av;
      double end_sum_w = (i+1) * total_w / N - sum_w; 
      if (i!=N-1) {
        while(av.sum_w() + iter->weight() <= end_sum_w) {
          av.add(iter->data(), iter->weight());
          ++iter;
        }
      }
      // use all remaining data for last bin (to avoid problems
      // due to rounding errors)
      else 
        add(av, iter, sortedvec.end());

      if (av.sum_w() == 0) {
        std::stringstream ss;
        ss << "yat::normalizer::qQuantileNormalizer: relative weight too "
           << "large. See qQuantileNormalizer constructor documentation "
           << "for details on weights.\n";
        throw std::runtime_error(ss.str());
      }
      average_(i) = av.mean();
      quantiles_(i)   = (sum_w + 0.5*av.sum_w())/total_w;
      sum_w += av.sum_w();
    }
  }


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


  const utility::Vector& qQuantileNormalizer::Partitioner::quantiles(void) const
  {
    return quantiles_;
  }


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

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