Changeset 1736


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

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

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/normalization_test.cc

    r1735 r1736  
    22
    33/*
    4   Copyright (C) 2008 Jari Häkkinen, Peter Johansson
    5   Copyright (C) 2009 Jari Häkkinen
     4  Copyright (C) 2008, 2009 Jari Häkkinen, Peter Johansson
    65
    76  This file is part of the yat library, http://dev.thep.lu.se/yat
     
    137136
    138137  suite.err() << "testing number of parts (Q) boundary conditions\n";
    139   qQuantileNormalizer(m.column_const_view(0),m.rows());
    140   qQuantileNormalizer(m.column_const_view(0),3);
    141 
    142   qQuantileNormalizer qqn(m.column_const_view(0),9);  // first column as target
     138  qQuantileNormalizer(m.begin_column(0), m.end_column(0), m.rows());
     139  qQuantileNormalizer(m.begin_column(0), m.end_column(0), 3);
     140
     141  // first column as target
     142  qQuantileNormalizer qqn(m.begin_column(0), m.end_column(0) ,9); 
     143  ColumnNormalizer<qQuantileNormalizer> cn(qqn);
    143144  utility::Matrix result(m.rows(),m.columns());
    144   qqn(m,result);
     145  cn(m, result);
    145146
    146147  suite.err() << "test that result can be stored in the source matrix\n";
    147   qqn(m,m);
     148  cn(m,m);
    148149  suite.add(result==m);
    149150
     
    156157  m2(2,0) = 1; m2(2,1) = 0;
    157158  m2(3,0) = 3; m2(3,1) = 7;
    158   qQuantileNormalizer qqn2(m2.column_const_view(0),m2.rows());
     159  qQuantileNormalizer qqn2(m2.begin_column(0), m2.end_column(0), m2.rows());
     160  ColumnNormalizer<qQuantileNormalizer> cn2(qqn2);
    159161  utility::Matrix result2(m2.rows(),m2.columns());
    160   qqn2(m2,result2);
     162  cn2(m2,result2);
    161163  suite.add( suite.equal_fix(m2(0,0),result2(2,1),1.0e-12) &&
    162164             suite.equal_fix(m2(1,0),result2(3,1),1.0e-12) &&
  • 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
  • trunk/yat/normalizer/qQuantileNormalizer.h

    r1733 r1736  
    33
    44/*
    5   Copyright (C) 2009 Jari Häkkinen
     5  Copyright (C) 2009 Jari Häkkinen, Peter Johansson
    66
    77  This file is part of the yat library, http://dev.thep.lu.se/yat
     
    2121*/
    2222
     23#include "yat/regression/CSplineInterpolation.h"
    2324#include "yat/utility/Vector.h"
     25#include "yat/utility/yat_assert.h"
     26
     27#include <algorithm>
     28#include <iterator>
     29#include <stdexcept>
    2430
    2531namespace theplu {
    2632namespace yat {
    2733namespace utility {
    28   class Matrix;
    2934  class VectorBase;
    3035}
     
    7883       normalization.
    7984    */
    80     qQuantileNormalizer(const utility::VectorBase& target, unsigned int Q);
     85    template<typename BidirectionalIterator>
     86    qQuantileNormalizer(BidirectionalIterator first, BidirectionalIterator last,
     87                        unsigned int Q);
    8188
    8289    /**
     
    8895       \note dimensions of \a matrix and \a result must match.
    8996     */
    90     void operator()(const utility::Matrix& matrix,
    91                     utility::Matrix& result) const;
     97    template<typename RandomAccessIterator1, typename RandomAccessIterator2>
     98    RandomAccessIterator2 operator()(RandomAccessIterator1 first,
     99                                     RandomAccessIterator1 last,
     100                                     RandomAccessIterator2 result) const;
    92101
    93102  private:
     
    106115       \brief Create the partition and perform required calculations.
    107116    */
    108     Partitioner(const utility::VectorBase& vec, unsigned int N);
     117    template<typename BidirectionalIterator>
     118    Partitioner(BidirectionalIterator first, BidirectionalIterator last,
     119                unsigned int N);
    109120
    110121    /**
     
    128139
    129140  private:
     141    void init(const utility::VectorBase&, unsigned int N);
     142
    130143    utility::Vector average_;
    131144    utility::Vector index_;
     
    136149  };
    137150
     151
     152  // template implementations
     153
     154  template<typename BidirectionalIterator>
     155  qQuantileNormalizer::qQuantileNormalizer(BidirectionalIterator first,
     156                                           BidirectionalIterator last,
     157                                           unsigned int Q)
     158    : target_(Partitioner(first, last, Q))
     159  {
     160    utility::yat_assert<std::runtime_error>(Q>2,
     161                                            "qQuantileNormalizer: Q too small");
     162  }
     163
     164
     165  template<typename RandomAccessIterator1, typename RandomAccessIterator2>
     166  RandomAccessIterator2
     167  qQuantileNormalizer::operator()(RandomAccessIterator1 first,
     168                                  RandomAccessIterator1 last,
     169                                  RandomAccessIterator2 result) const
     170  {
     171    size_t N = last-first;
     172    utility::yat_assert<std::runtime_error>
     173      (N >= target_.size(), "qQuantileNormalizer: Input range too small");
     174
     175    std::vector<size_t> sorted_index(last-first);
     176    utility::sort_index(first, last, sorted_index);
     177
     178    Partitioner source(first, last, target_.size());
     179    utility::Vector diff(source.averages());
     180    diff-=target_.averages();
     181    const utility::Vector& idx=target_.index();
     182    regression::CSplineInterpolation cspline(idx,diff);
     183
     184    // linear interpolation for first part, i.e., use first diff for
     185    // all points in the first part.
     186    size_t start=0;
     187    size_t end=static_cast<unsigned int>(idx(0));
     188    // The first condition below takes care of limiting case number
     189    // of parts approximately equal to the number of matrix rows and
     190    // the second condition makes sure that index is large enough
     191    // when using cspline below ... the static cast above takes the
     192    // floor whereas we want to take the "roof" forcing next index
     193    // range to be within interpolation range for the cspline.
     194    if ((end==0) || (end<idx(0)))
     195      ++end;
     196    for (size_t row=start; row<end; ++row) {
     197      size_t srow=sorted_index[row];
     198      result[srow] = first[srow] - diff(0);
     199    }
     200   
     201    // cspline interpolation for all data between the mid points of
     202    // the first and last part
     203    start=end;
     204    end=static_cast<unsigned int>(idx(target_.size()-1));
     205    // take care of limiting case number of parts approximately
     206    // equal to the number of matrix rows
     207    if (end==(static_cast<size_t>(N-1)) )
     208      --end;
     209    for (size_t row=start; row<=end; ++row) {
     210      size_t srow=sorted_index[row];
     211      result[srow] = first[srow] - cspline.evaluate(row) ;
     212    }
     213   
     214    // linear interpolation for last part, i.e., use last diff for
     215    // all points in the last part.
     216    start=end+1;
     217    end=N;
     218    for (size_t row=start; row<end; ++row) {
     219      size_t srow=sorted_index[row];
     220      result[srow] = first[srow] - diff(diff.size()-1);
     221    }
     222    return result + N;
     223  }
     224
     225
     226  template<typename BidirectionalIterator>
     227  qQuantileNormalizer::Partitioner::Partitioner(BidirectionalIterator first,
     228                                                BidirectionalIterator last,
     229                                                unsigned int N)
     230    : average_(utility::Vector(N)), index_(utility::Vector(N))
     231  {
     232    utility::Vector vec(std::distance(first, last));
     233    std::copy(first, last, vec.begin());
     234    std::sort(vec.begin(), vec.end());
     235    init(vec, N);
     236  }
     237
     238
    138239}}} // end of namespace normalizer, yat and thep
    139240
Note: See TracChangeset for help on using the changeset viewer.