Changeset 1736

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:

• test/normalization_test.cc (modified) (3 diffs)
• yat/normalizer/qQuantileNormalizer.cc (modified) (4 diffs)
• yat/normalizer/qQuantileNormalizer.h (modified) (7 diffs)

trunk/test/normalization_test.cc

@@ -2 +2 @@
 
 /*
-  Copyright (C) 2008 Jari Häkkinen, Peter Johansson
-  Copyright (C) 2009 Jari Häkkinen
+  Copyright (C) 2008, 2009 Jari Häkkinen, Peter Johansson
 
   This file is part of the yat library, http://dev.thep.lu.se/yat
@@ -137 +136 @@
 
   suite.err() << "testing number of parts (Q) boundary conditions\n";
-  qQuantileNormalizer(m.column_const_view(0),m.rows());
-  qQuantileNormalizer(m.column_const_view(0),3);
-
-  qQuantileNormalizer qqn(m.column_const_view(0),9);  // first column as target
+  qQuantileNormalizer(m.begin_column(0), m.end_column(0), m.rows());
+  qQuantileNormalizer(m.begin_column(0), m.end_column(0), 3);
+
+  // first column as target
+  qQuantileNormalizer qqn(m.begin_column(0), m.end_column(0) ,9);  
+  ColumnNormalizer<qQuantileNormalizer> cn(qqn);
   utility::Matrix result(m.rows(),m.columns());
-  qqn(m,result);
+  cn(m, result);
 
   suite.err() << "test that result can be stored in the source matrix\n";
-  qqn(m,m);
+  cn(m,m);
   suite.add(result==m);
 
@@ -156 +157 @@
   m2(2,0) = 1; m2(2,1) = 0;
   m2(3,0) = 3; m2(3,1) = 7;
-  qQuantileNormalizer qqn2(m2.column_const_view(0),m2.rows());
+  qQuantileNormalizer qqn2(m2.begin_column(0), m2.end_column(0), m2.rows());
+  ColumnNormalizer<qQuantileNormalizer> cn2(qqn2);
   utility::Matrix result2(m2.rows(),m2.columns());
-  qqn2(m2,result2);
+  cn2(m2,result2);
   suite.add( suite.equal_fix(m2(0,0),result2(2,1),1.0e-12) &&
              suite.equal_fix(m2(1,0),result2(3,1),1.0e-12) &&

trunk/yat/normalizer/qQuantileNormalizer.cc

@@ -2 +2 @@
 
 /*
-  Copyright (C) 2009 Jari Häkkinen
+  Copyright (C) 2009 Jari Häkkinen, Peter Johansson
 
   This file is part of the yat library, http://dev.thep.lu.se/yat
@@ -24 +24 @@
 #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"
@@ -36 +35 @@
 
 
-  qQuantileNormalizer::Partitioner::Partitioner(const utility::VectorBase& vec,
-                                                unsigned int N)
-    : average_(utility::Vector(N)), index_(utility::Vector(N))
+  void 
+  qQuantileNormalizer::Partitioner::init(const utility::VectorBase& sortedvec,
+                                         unsigned int N)
   {
     assert(N>1);
-    assert(N<=vec.size());
-    double range=static_cast<double>(vec.size())/N;
+    assert(N<=sortedvec.size());
+    double range=static_cast<double>(sortedvec.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) {
@@ -77 +74 @@
   }
 
-
-  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));
-      // 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 larege 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[column][row];
-        result(srow,column) = matrix(srow,column) - 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==(matrix.rows()-1))
-        --end;
-      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=end+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

trunk/yat/normalizer/qQuantileNormalizer.h

@@ -3 +3 @@
 
 /*
-  Copyright (C) 2009 Jari Häkkinen
+  Copyright (C) 2009 Jari Häkkinen, Peter Johansson
 
   This file is part of the yat library, http://dev.thep.lu.se/yat
@@ -21 +21 @@
 */
 
+#include "yat/regression/CSplineInterpolation.h"
 #include "yat/utility/Vector.h"
+#include "yat/utility/yat_assert.h"
+
+#include <algorithm>
+#include <iterator>
+#include <stdexcept>
 
 namespace theplu {
 namespace yat {
 namespace utility {
-  class Matrix;
   class VectorBase;
 }
@@ -78 +83 @@
        normalization.
     */
-    qQuantileNormalizer(const utility::VectorBase& target, unsigned int Q);
+    template<typename BidirectionalIterator>
+    qQuantileNormalizer(BidirectionalIterator first, BidirectionalIterator last,
+                        unsigned int Q);
 
     /**
@@ -88 +95 @@
        \note dimensions of \a matrix and \a result must match.
      */
-    void operator()(const utility::Matrix& matrix,
-                    utility::Matrix& result) const;
+    template<typename RandomAccessIterator1, typename RandomAccessIterator2>
+    RandomAccessIterator2 operator()(RandomAccessIterator1 first, 
+                                     RandomAccessIterator1 last,
+                                     RandomAccessIterator2 result) const;
 
   private:
@@ -106 +115 @@
        \brief Create the partition and perform required calculations.
     */
-    Partitioner(const utility::VectorBase& vec, unsigned int N);
+    template<typename BidirectionalIterator>
+    Partitioner(BidirectionalIterator first, BidirectionalIterator last, 
+                unsigned int N);
 
     /**
@@ -128 +139 @@
 
   private:
+    void init(const utility::VectorBase&, unsigned int N);
+
     utility::Vector average_;
     utility::Vector index_;
@@ -136 +149 @@
   };
 
+
+  // 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>
+  RandomAccessIterator2 
+  qQuantileNormalizer::operator()(RandomAccessIterator1 first, 
+                                  RandomAccessIterator1 last,
+                                  RandomAccessIterator2 result) 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);
+
+    Partitioner source(first, last, 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));
+    // 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);
+    }
+    return result + N;
+  }
+
+
+  template<typename BidirectionalIterator>
+  qQuantileNormalizer::Partitioner::Partitioner(BidirectionalIterator first, 
+                                                BidirectionalIterator last, 
+                                                unsigned int N)
+    : average_(utility::Vector(N)), index_(utility::Vector(N))
+  {
+    utility::Vector vec(std::distance(first, last));
+    std::copy(first, last, vec.begin());
+    std::sort(vec.begin(), vec.end());
+    init(vec, N);
+  }
+
+
 }}} // end of namespace normalizer, yat and thep