Changeset 1736
- Timestamp:
- Jan 19, 2009, 3:26:49 PM (13 years ago)
- Location:
- trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/test/normalization_test.cc
r1735 r1736 2 2 3 3 /* 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 6 5 7 6 This file is part of the yat library, http://dev.thep.lu.se/yat … … 137 136 138 137 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); 143 144 utility::Matrix result(m.rows(),m.columns()); 144 qqn(m,result);145 cn(m, result); 145 146 146 147 suite.err() << "test that result can be stored in the source matrix\n"; 147 qqn(m,m);148 cn(m,m); 148 149 suite.add(result==m); 149 150 … … 156 157 m2(2,0) = 1; m2(2,1) = 0; 157 158 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); 159 161 utility::Matrix result2(m2.rows(),m2.columns()); 160 qqn2(m2,result2);162 cn2(m2,result2); 161 163 suite.add( suite.equal_fix(m2(0,0),result2(2,1),1.0e-12) && 162 164 suite.equal_fix(m2(1,0),result2(3,1),1.0e-12) && -
trunk/yat/normalizer/qQuantileNormalizer.cc
r1735 r1736 2 2 3 3 /* 4 Copyright (C) 2009 Jari Häkkinen 4 Copyright (C) 2009 Jari Häkkinen, Peter Johansson 5 5 6 6 This file is part of the yat library, http://dev.thep.lu.se/yat … … 24 24 #include "yat/regression/CSplineInterpolation.h" 25 25 #include "yat/statistics/Averager.h" 26 #include "yat/utility/Matrix.h"27 26 #include "yat/utility/Vector.h" 28 27 #include "yat/utility/VectorBase.h" … … 36 35 37 36 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) 41 40 { 42 41 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; 45 44 assert(range); 46 utility::Vector sortedvec(vec);47 std::sort(sortedvec.begin(),sortedvec.end());48 45 unsigned int start=0; 49 46 for (unsigned int i=0; i<N; ++i) { … … 77 74 } 78 75 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 fit85 }86 87 88 void qQuantileNormalizer::operator()(const utility::Matrix& matrix,89 utility::Matrix& result) const90 {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 for108 // 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 number112 // of parts approximately equal to the number of matrix rows and113 // the second condition makes sure that index is larege enough114 // when using cspline below ... the static cast above takes the115 // floor whereas we want to take the "roof" forcing next index116 // 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 of125 // the first and last part126 start=end;127 end=static_cast<unsigned int>(idx(target_.size()-1));128 // take care of limiting case number of parts approximately129 // equal to the number of matrix rows130 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 for138 // 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 148 76 }}} // end of namespace normalizer, yat and thep -
trunk/yat/normalizer/qQuantileNormalizer.h
r1733 r1736 3 3 4 4 /* 5 Copyright (C) 2009 Jari Häkkinen 5 Copyright (C) 2009 Jari Häkkinen, Peter Johansson 6 6 7 7 This file is part of the yat library, http://dev.thep.lu.se/yat … … 21 21 */ 22 22 23 #include "yat/regression/CSplineInterpolation.h" 23 24 #include "yat/utility/Vector.h" 25 #include "yat/utility/yat_assert.h" 26 27 #include <algorithm> 28 #include <iterator> 29 #include <stdexcept> 24 30 25 31 namespace theplu { 26 32 namespace yat { 27 33 namespace utility { 28 class Matrix;29 34 class VectorBase; 30 35 } … … 78 83 normalization. 79 84 */ 80 qQuantileNormalizer(const utility::VectorBase& target, unsigned int Q); 85 template<typename BidirectionalIterator> 86 qQuantileNormalizer(BidirectionalIterator first, BidirectionalIterator last, 87 unsigned int Q); 81 88 82 89 /** … … 88 95 \note dimensions of \a matrix and \a result must match. 89 96 */ 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; 92 101 93 102 private: … … 106 115 \brief Create the partition and perform required calculations. 107 116 */ 108 Partitioner(const utility::VectorBase& vec, unsigned int N); 117 template<typename BidirectionalIterator> 118 Partitioner(BidirectionalIterator first, BidirectionalIterator last, 119 unsigned int N); 109 120 110 121 /** … … 128 139 129 140 private: 141 void init(const utility::VectorBase&, unsigned int N); 142 130 143 utility::Vector average_; 131 144 utility::Vector index_; … … 136 149 }; 137 150 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 138 239 }}} // end of namespace normalizer, yat and thep 139 240
Note: See TracChangeset
for help on using the changeset viewer.