Ignore:
Timestamp:
Jan 19, 2009, 3:26:49 PM (13 years ago)
Author:
Peter
Message:

Changing the interface to work on ranges rather than Matrix. This
allows usage within ColumnNormalizer? and RowNormalizer?. refs #425.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/yat/normalizer/qQuantileNormalizer.cc

    r1735 r1736  
    22
    33/*
    4   Copyright (C) 2009 Jari Häkkinen
     4  Copyright (C) 2009 Jari Häkkinen, Peter Johansson
    55
    66  This file is part of the yat library, http://dev.thep.lu.se/yat
     
    2424#include "yat/regression/CSplineInterpolation.h"
    2525#include "yat/statistics/Averager.h"
    26 #include "yat/utility/Matrix.h"
    2726#include "yat/utility/Vector.h"
    2827#include "yat/utility/VectorBase.h"
     
    3635
    3736
    38   qQuantileNormalizer::Partitioner::Partitioner(const utility::VectorBase& vec,
    39                                                 unsigned int N)
    40     : average_(utility::Vector(N)), index_(utility::Vector(N))
     37  void
     38  qQuantileNormalizer::Partitioner::init(const utility::VectorBase& sortedvec,
     39                                         unsigned int N)
    4140  {
    4241    assert(N>1);
    43     assert(N<=vec.size());
    44     double range=static_cast<double>(vec.size())/N;
     42    assert(N<=sortedvec.size());
     43    double range=static_cast<double>(sortedvec.size())/N;
    4544    assert(range);
    46     utility::Vector sortedvec(vec);
    47     std::sort(sortedvec.begin(),sortedvec.end());
    4845    unsigned int start=0;
    4946    for (unsigned int i=0; i<N; ++i) {
     
    7774  }
    7875
    79 
    80   qQuantileNormalizer::qQuantileNormalizer(const utility::VectorBase& target,
    81                                            unsigned int Q)
    82     : target_(Partitioner(target,Q))
    83   {
    84     assert(Q>2); // required by cspline fit
    85   }
    86 
    87 
    88   void qQuantileNormalizer::operator()(const utility::Matrix& matrix,
    89                                        utility::Matrix& result) const
    90   {
    91     assert(matrix.rows()    == result.rows());
    92     assert(matrix.columns() == result.columns());
    93     assert(matrix.rows()    >= target_.size());
    94 
    95     std::vector<std::vector<size_t> > sorted_index(matrix.rows());
    96     for (size_t column=0; column<matrix.columns(); ++column)
    97       utility::sort_index(sorted_index[column],
    98                           matrix.column_const_view(column));
    99 
    100     for (size_t column=0; column<matrix.columns(); ++column) {
    101       Partitioner source(matrix.column_const_view(column),target_.size());
    102       utility::Vector diff(source.averages());
    103       diff-=target_.averages();
    104       const utility::Vector& idx=target_.index();
    105       regression::CSplineInterpolation cspline(idx,diff);
    106 
    107       // linear interpolation for first part, i.e., use first diff for
    108       // all points in the first part.
    109       size_t start=0;
    110       size_t end=static_cast<unsigned int>(idx(0));
    111       // The first condition below takes care of limiting case number
    112       // of parts approximately equal to the number of matrix rows and
    113       // the second condition makes sure that index is larege enough
    114       // when using cspline below ... the static cast above takes the
    115       // floor whereas we want to take the "roof" forcing next index
    116       // range to be within interpolation range for the cspline.
    117       if ((end==0) || (end<idx(0)))
    118         ++end;
    119       for (size_t row=start; row<end; ++row) {
    120         size_t srow=sorted_index[column][row];
    121         result(srow,column) = matrix(srow,column) - diff(0);
    122       }
    123 
    124       // cspline interpolation for all data between the mid points of
    125       // the first and last part
    126       start=end;
    127       end=static_cast<unsigned int>(idx(target_.size()-1));
    128       // take care of limiting case number of parts approximately
    129       // equal to the number of matrix rows
    130       if (end==(matrix.rows()-1))
    131         --end;
    132       for (size_t row=start; row<=end; ++row) {
    133         size_t srow=sorted_index[column][row];
    134         result(srow,column) = matrix(srow,column) - cspline.evaluate(row) ;
    135       }
    136 
    137       // linear interpolation for last part, i.e., use last diff for
    138       // all points in the last part.
    139       start=end+1;
    140       end=result.rows();
    141       for (size_t row=start; row<end; ++row) {
    142         size_t srow=sorted_index[column][row];
    143         result(srow,column) = matrix(srow,column) - diff(diff.size()-1);
    144       }
    145     }
    146   }
    147 
    14876}}} // end of namespace normalizer, yat and thep
Note: See TracChangeset for help on using the changeset viewer.