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>
229 qQuantileNormalizer::qQuantileNormalizer(ForwardIterator first,
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) =
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,
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);
337 diff-=target_.averages();
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) ) {
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,
375 BOOST_CONCEPT_ASSERT((boost_concepts::ForwardTraversal<ForwardIterator>));
377 build(first, last, N,
382 template<
typename ForwardIterator>
383 void qQuantileNormalizer::Partitioner::build(ForwardIterator first,
384 ForwardIterator last,
389 std::copy(first, last, vec.
begin());
395 template<
typename ForwardIterator>
396 void qQuantileNormalizer::Partitioner::build(ForwardIterator first,
397 ForwardIterator last,
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
double evaluate(double x)
Calculate the interpolated value for x.
Definition: iterator_traits.h:412
The Department of Theoretical Physics namespace as we define it.
Definition: iterator_traits.h:47
void operator()(RandomAccessIterator1 first, RandomAccessIterator1 last, RandomAccessIterator2 result) const
Perform the Q-quantile normalization.
Definition: qQuantileNormalizer.h:239
Definition: iterator_traits.h:55
DataIterator.
Definition: DataIterator.h:63
This is the yat interface to GSL vector.
Definition: Vector.h:59
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
weight_reference weight(Iter iter) const
Definition: iterator_traits.h:446