1 #ifndef _theplu_yat_normalizer_qquantile_normalizer_ 2 #define _theplu_yat_normalizer_qquantile_normalizer_ 26 #include "yat/regression/CSplineInterpolation.h" 27 #include "yat/utility/concept_check.h" 28 #include "yat/utility/DataIterator.h" 29 #include "yat/utility/DataWeight.h" 30 #include "yat/utility/iterator_traits.h" 31 #include "yat/utility/sort_index.h" 32 #include "yat/utility/Vector.h" 33 #include "yat/utility/WeightIterator.h" 34 #include "yat/utility/yat_assert.h" 36 #include <boost/concept_check.hpp> 51 namespace normalizer {
116 template<
typename ForwardIterator>
142 template<
typename RandomAccessIterator1,
typename RandomAccessIterator2>
143 void operator()(RandomAccessIterator1 first, RandomAccessIterator1 last,
144 RandomAccessIterator2 result)
const;
163 template<
typename ForwardIterator>
164 Partitioner(ForwardIterator first, ForwardIterator last,
192 size_t size(
void)
const;
196 template<
typename ForwardIterator>
197 void build(ForwardIterator first, ForwardIterator last,
unsigned int N,
200 template<
typename ForwardIterator>
201 void build(ForwardIterator first, ForwardIterator last,
unsigned int N,
204 void init(
const std::vector<utility::DataWeight>&,
unsigned int N);
213 template<
typename RandomAccessIterator1,
typename RandomAccessIterator2>
214 void normalize(
const Partitioner& source,RandomAccessIterator1 first,
215 RandomAccessIterator1 last, RandomAccessIterator2 result,
219 template<
typename RandomAccessIterator1,
typename RandomAccessIterator2>
220 void normalize(
const Partitioner& source,RandomAccessIterator1 first,
221 RandomAccessIterator1 last, RandomAccessIterator2 result,
228 template<
typename ForwardIterator>
230 ForwardIterator last,
232 : target_(Partitioner(first, last, Q))
238 template<
typename RandomAccessIterator1,
typename RandomAccessIterator2>
240 RandomAccessIterator1 last,
241 RandomAccessIterator2 result)
const 244 using boost_concepts::RandomAccessTraversal;
245 BOOST_CONCEPT_ASSERT((RandomAccessTraversal<RandomAccessIterator1>));
247 BOOST_CONCEPT_ASSERT((RandomAccessTraversal<RandomAccessIterator2>));
248 using boost_concepts::WritableIterator;
249 BOOST_CONCEPT_ASSERT((WritableIterator<RandomAccessIterator2>));
251 Partitioner source(first, last, target_.size());
253 normalize(source, first, last, result, tag);
257 template<
typename RandomAccessIterator1,
typename RandomAccessIterator2>
259 qQuantileNormalizer::normalize(
const qQuantileNormalizer::Partitioner& source,
260 RandomAccessIterator1 first,
261 RandomAccessIterator1 last,
262 RandomAccessIterator2 result,
267 detail::copy_weight_if_weighted(first, last, result);
269 size_t N = last-first;
270 YAT_ASSERT(N >= target_.size());
272 std::vector<size_t> sorted_index(last-first);
276 diff-=target_.averages();
283 size_t end =
static_cast<size_t>(std::ceil(N*idx(0) - 0.5));
290 for (
size_t i=start; i<end; ++i) {
291 size_t si = sorted_index[i];
292 traits2.
data(result+si) = traits1.
data(first+si) - diff(0);
295 using utility::yat_assert;
300 end =
static_cast<size_t>(std::ceil(N*idx(idx.
size()-1) - 0.5));
303 for (
size_t i=start; i<end; ++i) {
304 size_t si = sorted_index[i];
306 YAT_ASSERT((i+0.5)/N>idx(0));
307 traits2.
data(result+si) =
308 traits1.
data(first+si) - cspline.evaluate((i+0.5)/N);
313 for (
size_t i=end ; i<N; ++i) {
314 size_t si = sorted_index[i];
315 traits2.
data(result+si) = traits1.
data(first+si) - diff(diff.size()-1);
320 template<
typename RandomAccessIterator1,
typename RandomAccessIterator2>
321 void qQuantileNormalizer::normalize(
const Partitioner& source,
322 RandomAccessIterator1 first,
323 RandomAccessIterator1 last,
324 RandomAccessIterator2 result,
325 utility::weighted_iterator_tag tag)
const 328 detail::copy_weight_if_weighted(first, last, result);
330 double total_w = std::accumulate(utility::weight_iterator(first),
331 utility::weight_iterator(last), 0.0);
333 std::vector<size_t> sorted_index(last-first);
336 utility::Vector diff(source.averages());
337 diff-=target_.averages();
338 const utility::Vector& idx=target_.quantiles();
339 regression::CSplineInterpolation cspline(idx,diff);
342 utility::iterator_traits<RandomAccessIterator1> traits1;
343 utility::iterator_traits<RandomAccessIterator2> traits2;
344 for (
size_t i=0; i<sorted_index.size(); ++i) {
345 size_t si = sorted_index[i];
346 double w = (sum_w + 0.5*traits1.weight(first+si))/total_w;
347 double correction = 0;
351 correction = diff(0);
353 else if (w < idx(idx.size()-1) ) {
356 correction = cspline.evaluate(w);
361 correction = diff(diff.size()-1);
363 traits2.data(result+si) = traits1.data(first+si) - correction;
364 sum_w += traits1.weight(first+si);
369 template<
typename ForwardIterator>
370 qQuantileNormalizer::Partitioner::Partitioner(ForwardIterator first,
371 ForwardIterator last,
373 : average_(utility::Vector(N)), quantiles_(utility::Vector(N))
375 BOOST_CONCEPT_ASSERT((boost_concepts::ForwardTraversal<ForwardIterator>));
376 BOOST_CONCEPT_ASSERT((utility::DataIterator<ForwardIterator>));
377 build(first, last, N,
382 template<
typename ForwardIterator>
383 void qQuantileNormalizer::Partitioner::build(ForwardIterator first,
384 ForwardIterator last,
386 utility::unweighted_iterator_tag)
388 utility::Vector vec(std::distance(first, last));
389 std::copy(first, last, vec.begin());
390 std::sort(vec.begin(), vec.end());
395 template<
typename ForwardIterator>
396 void qQuantileNormalizer::Partitioner::build(ForwardIterator first,
397 ForwardIterator last,
399 utility::weighted_iterator_tag)
401 using namespace utility;
402 std::vector<DataWeight> vec(first, last);
403 std::sort(vec.begin(), vec.end());
detail::weighted_iterator_traits_detail< value >::type type
Definition: iterator_traits.h:114
Perform Q-quantile normalization.
Definition: qQuantileNormalizer.h:87
Cubic spline with natural boundary conditions.
Definition: CSplineInterpolation.h:47
data_reference data(Iter iter) const
Definition: iterator_traits.h:440
Definition: iterator_traits.h:412
The Department of Theoretical Physics namespace as we define it.
Definition: iterator_traits.h:47
Definition: iterator_traits.h:55
DataIterator.
Definition: DataIterator.h:63
void operator()(RandomAccessIterator1 first, RandomAccessIterator1 last, RandomAccessIterator2 result) const
Perform the Q-quantile normalization.
Definition: qQuantileNormalizer.h:239
This is the yat interface to GSL vector.
Definition: Vector.h:59
qQuantileNormalizer(ForwardIterator first, ForwardIterator last, unsigned int Q)
Constructor.
Definition: qQuantileNormalizer.h:229
This is the yat interface to GSL vector.
Definition: VectorBase.h:55
void check_iterator_is_unweighted(Iter iter)
check (at compile time) that iterator is unweighted.
Definition: iterator_traits.h:200
void sort_index(InputIterator first, InputIterator last, std::vector< size_t > &sort_index)
Definition: sort_index.h:146