yat/normalizer/qQuantileNormalizer.h

Code
Comments
Other
Rev Date Author Line
1708 13 Jan 09 jari 1 #ifndef _theplu_yat_normalizer_qquantile_normalizer_
1708 13 Jan 09 jari 2 #define _theplu_yat_normalizer_qquantile_normalizer_
1708 13 Jan 09 jari 3
1708 13 Jan 09 jari 4 /*
2119 12 Dec 09 peter 5   Copyright (C) 2009 Jari Häkkinen, Peter Johansson
3572 12 Jan 17 peter 6   Copyright (C) 2010, 2016, 2017 Peter Johansson
1708 13 Jan 09 jari 7
1708 13 Jan 09 jari 8   This file is part of the yat library, http://dev.thep.lu.se/yat
1708 13 Jan 09 jari 9
1708 13 Jan 09 jari 10   The yat library is free software; you can redistribute it and/or
1708 13 Jan 09 jari 11   modify it under the terms of the GNU General Public License as
1708 13 Jan 09 jari 12   published by the Free Software Foundation; either version 3 of the
1708 13 Jan 09 jari 13   License, or (at your option) any later version.
1708 13 Jan 09 jari 14
1708 13 Jan 09 jari 15   The yat library is distributed in the hope that it will be useful,
1708 13 Jan 09 jari 16   but WITHOUT ANY WARRANTY; without even the implied warranty of
1708 13 Jan 09 jari 17   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1708 13 Jan 09 jari 18   General Public License for more details.
1708 13 Jan 09 jari 19
1708 13 Jan 09 jari 20   You should have received a copy of the GNU General Public License
1708 13 Jan 09 jari 21   along with yat. If not, see <http://www.gnu.org/licenses/>.
1708 13 Jan 09 jari 22 */
1708 13 Jan 09 jari 23
2158 18 Jan 10 peter 24 #include "utility.h"
2158 18 Jan 10 peter 25
1736 19 Jan 09 peter 26 #include "yat/regression/CSplineInterpolation.h"
2263 26 May 10 peter 27 #include "yat/utility/concept_check.h"
1740 21 Jan 09 peter 28 #include "yat/utility/DataIterator.h"
1737 19 Jan 09 peter 29 #include "yat/utility/DataWeight.h"
1737 19 Jan 09 peter 30 #include "yat/utility/iterator_traits.h"
1777 05 Feb 09 jari 31 #include "yat/utility/sort_index.h"
1708 13 Jan 09 jari 32 #include "yat/utility/Vector.h"
1740 21 Jan 09 peter 33 #include "yat/utility/WeightIterator.h"
1736 19 Jan 09 peter 34 #include "yat/utility/yat_assert.h"
1708 13 Jan 09 jari 35
2150 17 Jan 10 peter 36 #include <boost/concept_check.hpp>
2150 17 Jan 10 peter 37
1736 19 Jan 09 peter 38 #include <algorithm>
1755 27 Jan 09 peter 39 #include <cmath>
1736 19 Jan 09 peter 40 #include <iterator>
2055 08 Sep 09 peter 41 #include <limits>
1778 06 Feb 09 peter 42 #include <numeric>
1736 19 Jan 09 peter 43 #include <stdexcept>
1737 19 Jan 09 peter 44 #include <vector>
1736 19 Jan 09 peter 45
1708 13 Jan 09 jari 46 namespace theplu {
1708 13 Jan 09 jari 47 namespace yat {
1708 13 Jan 09 jari 48 namespace utility {
1712 13 Jan 09 jari 49   class VectorBase;
1708 13 Jan 09 jari 50 }
1708 13 Jan 09 jari 51 namespace normalizer {
1708 13 Jan 09 jari 52
1708 13 Jan 09 jari 53   /**
1708 13 Jan 09 jari 54      \brief Perform Q-quantile normalization
1708 13 Jan 09 jari 55
1810 20 Feb 09 peter 56      Perform a Q-quantile normalization on a \a source range, after
1810 20 Feb 09 peter 57      which it will approximately have the same distribution of data as
1810 20 Feb 09 peter 58      the \a target range (the Q-quantiles are the same). The rank of
1810 20 Feb 09 peter 59      an element in the \a source range is not changed.
1708 13 Jan 09 jari 60
1810 20 Feb 09 peter 61      The class works also with unweighed ranges, and there is no
1810 20 Feb 09 peter 62      restriction that weighted \a source range requires weighted \a
1810 20 Feb 09 peter 63      target range or vice versa.
1810 20 Feb 09 peter 64
1810 20 Feb 09 peter 65      Normalization goes like this:
1810 20 Feb 09 peter 66      - Data are not assumed to be sorted.
1819 23 Feb 09 jari 67      - Partition sorted \a target data in Q parts. Q must be 3 or larger
1954 07 May 09 jari 68        because of requirements from the underlying cubic spline fit
1810 20 Feb 09 peter 69      - Calculate the arithmetic (weighted) mean for each part, the mean is
1712 13 Jan 09 jari 70        assigned to the mid point of each part.
1954 07 May 09 jari 71      - Do the above for the data to be transformed (called \a source
1712 13 Jan 09 jari 72        here).
1810 20 Feb 09 peter 73      - For each part, calculate the difference between the \a target
1819 23 Feb 09 jari 74        and \a the source. Now we have \a Q differences \f$ d_i \f$
1810 20 Feb 09 peter 75        with associated rank (midpoint of each part).
1810 20 Feb 09 peter 76      - Create a cubic spline fit to this difference vector \a d. The
1810 20 Feb 09 peter 77        resulting curve is used to recalculate all values in \a source.
1712 13 Jan 09 jari 78        - Use the cubic spline fit for values within the cubic spline
1712 13 Jan 09 jari 79          fit range [midpoint 1st part, midpoint last part].
1712 13 Jan 09 jari 80        - For data outside the cubic spline fit use linear
1810 20 Feb 09 peter 81          extrapolation, i.e., a constant shift. \f$ d_{first} \f$ for
1811 20 Feb 09 peter 82          points below fit range, and \f$ d_{last} \f$ for points above fit
1810 20 Feb 09 peter 83          range.
1708 13 Jan 09 jari 84
1708 13 Jan 09 jari 85      \since New in yat 0.5
1708 13 Jan 09 jari 86    */
1708 13 Jan 09 jari 87   class qQuantileNormalizer
1708 13 Jan 09 jari 88   {
1708 13 Jan 09 jari 89   public:
1708 13 Jan 09 jari 90     /**
1954 07 May 09 jari 91        \brief Constructor
1708 13 Jan 09 jari 92
1812 20 Feb 09 peter 93        Divides a sorted copy of range [\a first,\a last) into \a Q
1812 20 Feb 09 peter 94        parts. Parts are divided such that the sum of weights is
1812 20 Feb 09 peter 95        approximately the same in the different parts. If a relative
1812 20 Feb 09 peter 96        weight, \f$ w_i / \sum w_i \f$, is larger than 1/Q this might
1812 20 Feb 09 peter 97        be difficult to achieve, in which case a an exception is
1812 20 Feb 09 peter 98        thrown. In the unweighted case this implies that \a Q should be
1812 20 Feb 09 peter 99        smaller (or equal) than number of elements in [\a first, \a
1812 20 Feb 09 peter 100        last).
1812 20 Feb 09 peter 101
1954 07 May 09 jari 102        The range supplied to the constructor sets the target
1954 07 May 09 jari 103        distribution.
1812 20 Feb 09 peter 104
1954 07 May 09 jari 105        As the \a source range is also divided into \a Q parts (when
1954 07 May 09 jari 106        operator() is called), it is recommended to keep \a Q smaller
1954 07 May 09 jari 107        (or equal) than the size of ranges that will be normalized.
1954 07 May 09 jari 108
1812 20 Feb 09 peter 109        Also, \a Q must not be smaller than 3 due to restrictions in
1954 07 May 09 jari 110        the cubic spline fit utilized in the normalization.
1812 20 Feb 09 peter 111
1821 24 Feb 09 peter 112        \b Type \b Requirements:
3547 30 Dec 16 peter 113        - \c ForwardIterator is a model of \forward_traversal_iterator
2141 12 Jan 10 peter 114        - \c ForwardIterator is a \ref concept_data_iterator
1708 13 Jan 09 jari 115     */
1801 17 Feb 09 peter 116     template<typename ForwardIterator>
1801 17 Feb 09 peter 117     qQuantileNormalizer(ForwardIterator first, ForwardIterator last,
1736 19 Jan 09 peter 118                         unsigned int Q);
1708 13 Jan 09 jari 119
1708 13 Jan 09 jari 120     /**
1821 24 Feb 09 peter 121        \brief Perform the Q-quantile normalization.
1708 13 Jan 09 jari 122
1812 20 Feb 09 peter 123        Elements in [\a first, \a last) are normalized as described
1812 20 Feb 09 peter 124        above and the result is assigned to [\a result, \a result + \a
1814 20 Feb 09 peter 125        last-\a first). Input range [\a first, \a last) is not
1814 20 Feb 09 peter 126        modified. If ranges are weighted, the weights are copied from
1814 20 Feb 09 peter 127        [\a first, \a last) to \a result range.
1812 20 Feb 09 peter 128
1708 13 Jan 09 jari 129        It is possible to normalize "in place"; it is permissible for
1812 20 Feb 09 peter 130        \a first and \a result to be the same. However, as assignment
1812 20 Feb 09 peter 131        occurs sequentially, the operation is undefined if \a result is
1812 20 Feb 09 peter 132        the same as any in range [\a first, \a last).
1812 20 Feb 09 peter 133
1821 24 Feb 09 peter 134        \b Type Requirements:
2141 12 Jan 10 peter 135        - \c RandomAccessIterator1 is a \ref concept_data_iterator
3547 30 Dec 16 peter 136        - \c RandomAccessIterator1 is a \random_access_traversal_iterator
3547 30 Dec 16 peter 137        - \c RandomAccessIterator2 is a \ref concept_data_iterator
3547 30 Dec 16 peter 138        - \c RandomAccessIterator2 is a \random_access_traversal_iterator
3547 30 Dec 16 peter 139        - \c RandomAccessIterator2 is a \writable_iterator
1814 20 Feb 09 peter 140
1708 13 Jan 09 jari 141      */
1736 19 Jan 09 peter 142     template<typename RandomAccessIterator1, typename RandomAccessIterator2>
1739 21 Jan 09 peter 143     void operator()(RandomAccessIterator1 first, RandomAccessIterator1 last,
1739 21 Jan 09 peter 144                     RandomAccessIterator2 result) const;
1708 13 Jan 09 jari 145
1708 13 Jan 09 jari 146   private:
1716 13 Jan 09 jari 147
1716 13 Jan 09 jari 148   /**
1751 26 Jan 09 peter 149      \brief Partition a range of data into equal sizes.
1716 13 Jan 09 jari 150
1814 20 Feb 09 peter 151      Copy the range [first, last), sort the copy, and divide the
1819 23 Feb 09 jari 152      sorted copy in Q parts. The parts are created such that the total
1819 23 Feb 09 jari 153      weight in a part is approximately W/Q where W is the total weight
1814 20 Feb 09 peter 154      (over all parts). The class calculates the average value in each
3545 23 Dec 16 peter 155      part and also the "quantile".
1716 13 Jan 09 jari 156   */
1716 13 Jan 09 jari 157   class Partitioner
1716 13 Jan 09 jari 158   {
1716 13 Jan 09 jari 159   public:
1716 13 Jan 09 jari 160     /**
1716 13 Jan 09 jari 161        \brief Create the partition and perform required calculations.
1716 13 Jan 09 jari 162     */
1801 17 Feb 09 peter 163     template<typename ForwardIterator>
3545 23 Dec 16 peter 164     Partitioner(ForwardIterator first, ForwardIterator last,
1736 19 Jan 09 peter 165                 unsigned int N);
1716 13 Jan 09 jari 166
1716 13 Jan 09 jari 167     /**
1716 13 Jan 09 jari 168        \brief Return the averages for each part.
1716 13 Jan 09 jari 169
1716 13 Jan 09 jari 170        \return The average vector.
1716 13 Jan 09 jari 171     */
1716 13 Jan 09 jari 172     const utility::Vector& averages(void) const;
1716 13 Jan 09 jari 173
1716 13 Jan 09 jari 174     /**
1814 20 Feb 09 peter 175        The quantile (here) is defined as (w_lower + w_upper) / 2W,
1814 20 Feb 09 peter 176        where w_lower is the total weight of elements smaller than the
1814 20 Feb 09 peter 177        smallest element in the part, and w_upper is the total weight
1814 20 Feb 09 peter 178        of elements smaller (or equal) than the largest value in the
3545 23 Dec 16 peter 179        part.
1716 13 Jan 09 jari 180
1814 20 Feb 09 peter 181        In the unweighted case all weights are 1.0, which implies q_0 =
1814 20 Feb 09 peter 182        n_0/N, q_1 = (n_0+n_1/2)/N, q_2 = (n_0+n_1+n_2/2)/N where n_i
1814 20 Feb 09 peter 183        is number of elements in ith part.
1814 20 Feb 09 peter 184
1813 20 Feb 09 peter 185        \return The quantiles vector.
1716 13 Jan 09 jari 186     */
1813 20 Feb 09 peter 187     const utility::Vector& quantiles(void) const;
1716 13 Jan 09 jari 188
1716 13 Jan 09 jari 189     /**
1716 13 Jan 09 jari 190        \return The number of parts.
1716 13 Jan 09 jari 191     */
1716 13 Jan 09 jari 192     size_t size(void) const;
1716 13 Jan 09 jari 193
1716 13 Jan 09 jari 194   private:
1737 19 Jan 09 peter 195     // unweighted "constructor"
1801 17 Feb 09 peter 196     template<typename ForwardIterator>
3545 23 Dec 16 peter 197     void build(ForwardIterator first, ForwardIterator last, unsigned int N,
1737 19 Jan 09 peter 198                utility::unweighted_iterator_tag);
1738 20 Jan 09 peter 199     // weighted "constructor"
1801 17 Feb 09 peter 200     template<typename ForwardIterator>
3545 23 Dec 16 peter 201     void build(ForwardIterator first, ForwardIterator last, unsigned int N,
1738 20 Jan 09 peter 202                utility::weighted_iterator_tag);
1736 19 Jan 09 peter 203     void init(const utility::VectorBase&, unsigned int N);
1738 20 Jan 09 peter 204     void init(const std::vector<utility::DataWeight>&, unsigned int N);
1736 19 Jan 09 peter 205
1716 13 Jan 09 jari 206     utility::Vector average_;
1813 20 Feb 09 peter 207     utility::Vector quantiles_;
1716 13 Jan 09 jari 208   };
1716 13 Jan 09 jari 209
1708 13 Jan 09 jari 210     Partitioner target_;
1740 21 Jan 09 peter 211
1740 21 Jan 09 peter 212     // unweighted version
1740 21 Jan 09 peter 213     template<typename RandomAccessIterator1, typename RandomAccessIterator2>
3545 23 Dec 16 peter 214     void normalize(const Partitioner& source,RandomAccessIterator1 first,
1740 21 Jan 09 peter 215                    RandomAccessIterator1 last, RandomAccessIterator2 result,
1740 21 Jan 09 peter 216                    utility::unweighted_iterator_tag tag) const;
1740 21 Jan 09 peter 217
1740 21 Jan 09 peter 218     // weighted version
1740 21 Jan 09 peter 219     template<typename RandomAccessIterator1, typename RandomAccessIterator2>
3545 23 Dec 16 peter 220     void normalize(const Partitioner& source,RandomAccessIterator1 first,
1740 21 Jan 09 peter 221                    RandomAccessIterator1 last, RandomAccessIterator2 result,
1740 21 Jan 09 peter 222                    utility::weighted_iterator_tag tag) const;
1708 13 Jan 09 jari 223   };
1708 13 Jan 09 jari 224
1736 19 Jan 09 peter 225
1736 19 Jan 09 peter 226   // template implementations
1736 19 Jan 09 peter 227
1801 17 Feb 09 peter 228   template<typename ForwardIterator>
3545 23 Dec 16 peter 229   qQuantileNormalizer::qQuantileNormalizer(ForwardIterator first,
1801 17 Feb 09 peter 230                                            ForwardIterator last,
1736 19 Jan 09 peter 231                                            unsigned int Q)
1778 06 Feb 09 peter 232     : target_(Partitioner(first, last, Q))
1736 19 Jan 09 peter 233   {
2210 05 Mar 10 peter 234     YAT_ASSERT(Q>2);
1736 19 Jan 09 peter 235   }
1736 19 Jan 09 peter 236
1736 19 Jan 09 peter 237
1736 19 Jan 09 peter 238   template<typename RandomAccessIterator1, typename RandomAccessIterator2>
3545 23 Dec 16 peter 239   void qQuantileNormalizer::operator()(RandomAccessIterator1 first,
1739 21 Jan 09 peter 240                                        RandomAccessIterator1 last,
1739 21 Jan 09 peter 241                                        RandomAccessIterator2 result) const
1736 19 Jan 09 peter 242   {
3547 30 Dec 16 peter 243     BOOST_CONCEPT_ASSERT((utility::DataIterator<RandomAccessIterator1>));
3547 30 Dec 16 peter 244     using boost_concepts::RandomAccessTraversal;
3547 30 Dec 16 peter 245     BOOST_CONCEPT_ASSERT((RandomAccessTraversal<RandomAccessIterator1>));
3547 30 Dec 16 peter 246     BOOST_CONCEPT_ASSERT((utility::DataIterator<RandomAccessIterator2>));
3547 30 Dec 16 peter 247     BOOST_CONCEPT_ASSERT((RandomAccessTraversal<RandomAccessIterator2>));
3547 30 Dec 16 peter 248     using boost_concepts::WritableIterator;
3547 30 Dec 16 peter 249     BOOST_CONCEPT_ASSERT((WritableIterator<RandomAccessIterator2>));
3547 30 Dec 16 peter 250
1740 21 Jan 09 peter 251     Partitioner source(first, last, target_.size());
3572 12 Jan 17 peter 252     typename utility::weighted_iterator_traits<RandomAccessIterator1>::type tag;
1740 21 Jan 09 peter 253     normalize(source, first, last, result, tag);
1740 21 Jan 09 peter 254   }
1740 21 Jan 09 peter 255
1740 21 Jan 09 peter 256
1740 21 Jan 09 peter 257   template<typename RandomAccessIterator1, typename RandomAccessIterator2>
3545 23 Dec 16 peter 258   void
1740 21 Jan 09 peter 259   qQuantileNormalizer::normalize(const qQuantileNormalizer::Partitioner& source,
3545 23 Dec 16 peter 260                                  RandomAccessIterator1 first,
3545 23 Dec 16 peter 261                                  RandomAccessIterator1 last,
1740 21 Jan 09 peter 262                                  RandomAccessIterator2 result,
1740 21 Jan 09 peter 263                                  utility::unweighted_iterator_tag tag) const
1740 21 Jan 09 peter 264   {
1750 26 Jan 09 peter 265     utility::check_iterator_is_unweighted(first);
3572 12 Jan 17 peter 266     // copy the weights if needed
3572 12 Jan 17 peter 267     detail::copy_weight_if_weighted(first, last, result);
3572 12 Jan 17 peter 268
1736 19 Jan 09 peter 269     size_t N = last-first;
2210 05 Mar 10 peter 270     YAT_ASSERT(N >= target_.size());
1736 19 Jan 09 peter 271
1736 19 Jan 09 peter 272     std::vector<size_t> sorted_index(last-first);
1736 19 Jan 09 peter 273     utility::sort_index(first, last, sorted_index);
1736 19 Jan 09 peter 274
1736 19 Jan 09 peter 275     utility::Vector diff(source.averages());
1736 19 Jan 09 peter 276     diff-=target_.averages();
1813 20 Feb 09 peter 277     const utility::Vector& idx=target_.quantiles();
1736 19 Jan 09 peter 278     regression::CSplineInterpolation cspline(idx,diff);
1736 19 Jan 09 peter 279
1768 03 Feb 09 jari 280     // linear extrapolation for first part, i.e., use first diff for
1736 19 Jan 09 peter 281     // all points in the first part.
1736 19 Jan 09 peter 282     size_t start=0;
1778 06 Feb 09 peter 283     size_t end = static_cast<size_t>(std::ceil(N*idx(0) - 0.5));
1755 27 Jan 09 peter 284     // take care of limiting case number of parts approximately equal
1755 27 Jan 09 peter 285     // to the number of elements in range.
3572 12 Jan 17 peter 286     utility::iterator_traits<RandomAccessIterator1> traits1;
3572 12 Jan 17 peter 287     utility::iterator_traits<RandomAccessIterator2> traits2;
1755 27 Jan 09 peter 288     if (end==0)
1736 19 Jan 09 peter 289       ++end;
1778 06 Feb 09 peter 290     for (size_t i=start; i<end; ++i) {
1778 06 Feb 09 peter 291       size_t si = sorted_index[i];
3572 12 Jan 17 peter 292       traits2.data(result+si) = traits1.data(first+si) - diff(0);
1736 19 Jan 09 peter 293     }
1778 06 Feb 09 peter 294
1778 06 Feb 09 peter 295     using utility::yat_assert;
3545 23 Dec 16 peter 296
1736 19 Jan 09 peter 297     // cspline interpolation for all data between the mid points of
1736 19 Jan 09 peter 298     // the first and last part
1736 19 Jan 09 peter 299     start=end;
1778 06 Feb 09 peter 300     end = static_cast<size_t>(std::ceil(N*idx(idx.size()-1) - 0.5));
1778 06 Feb 09 peter 301     if (end>=N)
1778 06 Feb 09 peter 302       end = N-1;
1778 06 Feb 09 peter 303     for ( size_t i=start; i<end; ++i) {
1778 06 Feb 09 peter 304       size_t si = sorted_index[i];
3545 23 Dec 16 peter 305
2210 05 Mar 10 peter 306       YAT_ASSERT((i+0.5)/N>idx(0));
3572 12 Jan 17 peter 307       traits2.data(result+si) =
3572 12 Jan 17 peter 308         traits1.data(first+si) - cspline.evaluate((i+0.5)/N);
1736 19 Jan 09 peter 309     }
3545 23 Dec 16 peter 310
1768 03 Feb 09 jari 311     // linear extrapolation for last part, i.e., use last diff for
1736 19 Jan 09 peter 312     // all points in the last part.
1778 06 Feb 09 peter 313     for (size_t i=end ; i<N; ++i) {
1778 06 Feb 09 peter 314       size_t si = sorted_index[i];
3572 12 Jan 17 peter 315       traits2.data(result+si) = traits1.data(first+si) - diff(diff.size()-1);
1736 19 Jan 09 peter 316     }
1736 19 Jan 09 peter 317   }
1736 19 Jan 09 peter 318
1736 19 Jan 09 peter 319
1740 21 Jan 09 peter 320   template<typename RandomAccessIterator1, typename RandomAccessIterator2>
1740 21 Jan 09 peter 321   void qQuantileNormalizer::normalize(const Partitioner& source,
3545 23 Dec 16 peter 322                                       RandomAccessIterator1 first,
3545 23 Dec 16 peter 323                                       RandomAccessIterator1 last,
1740 21 Jan 09 peter 324                                       RandomAccessIterator2 result,
1740 21 Jan 09 peter 325                                       utility::weighted_iterator_tag tag) const
1740 21 Jan 09 peter 326   {
1740 21 Jan 09 peter 327     // copy the weights
2158 18 Jan 10 peter 328     detail::copy_weight_if_weighted(first, last, result);
1778 06 Feb 09 peter 329
1778 06 Feb 09 peter 330     double total_w = std::accumulate(utility::weight_iterator(first),
1778 06 Feb 09 peter 331                                      utility::weight_iterator(last), 0.0);
1778 06 Feb 09 peter 332
1778 06 Feb 09 peter 333     std::vector<size_t> sorted_index(last-first);
3571 10 Jan 17 peter 334     utility::sort_index(first, last, sorted_index);
1778 06 Feb 09 peter 335
1778 06 Feb 09 peter 336     utility::Vector diff(source.averages());
1778 06 Feb 09 peter 337     diff-=target_.averages();
1813 20 Feb 09 peter 338     const utility::Vector& idx=target_.quantiles();
1778 06 Feb 09 peter 339     regression::CSplineInterpolation cspline(idx,diff);
1778 06 Feb 09 peter 340
1778 06 Feb 09 peter 341     double sum_w=0;
1803 18 Feb 09 peter 342     utility::iterator_traits<RandomAccessIterator1> traits1;
1778 06 Feb 09 peter 343     utility::iterator_traits<RandomAccessIterator2> traits2;
1778 06 Feb 09 peter 344     for (size_t i=0; i<sorted_index.size(); ++i) {
1778 06 Feb 09 peter 345       size_t si = sorted_index[i];
1778 06 Feb 09 peter 346       double w = (sum_w + 0.5*traits1.weight(first+si))/total_w;
1778 06 Feb 09 peter 347       double correction = 0;
1778 06 Feb 09 peter 348       if (w <= idx(0)) {
1778 06 Feb 09 peter 349         // linear extrapolation for first part, i.e., use first diff for
1778 06 Feb 09 peter 350         // all points in the first part.
1778 06 Feb 09 peter 351         correction = diff(0);
1778 06 Feb 09 peter 352       }
1778 06 Feb 09 peter 353       else if (w < idx(idx.size()-1) ) {
1778 06 Feb 09 peter 354         // cspline interpolation for all data between the mid points of
1778 06 Feb 09 peter 355         // the first and last part
1778 06 Feb 09 peter 356         correction = cspline.evaluate(w);
1778 06 Feb 09 peter 357       }
1778 06 Feb 09 peter 358       else {
1778 06 Feb 09 peter 359         // linear extrapolation for last part, i.e., use last diff for
1778 06 Feb 09 peter 360         // all points in the last part.
1778 06 Feb 09 peter 361         correction = diff(diff.size()-1);
1778 06 Feb 09 peter 362       }
1778 06 Feb 09 peter 363       traits2.data(result+si) = traits1.data(first+si) - correction;
1778 06 Feb 09 peter 364       sum_w += traits1.weight(first+si);
1778 06 Feb 09 peter 365     }
1740 21 Jan 09 peter 366   }
1740 21 Jan 09 peter 367
1740 21 Jan 09 peter 368
1801 17 Feb 09 peter 369   template<typename ForwardIterator>
3545 23 Dec 16 peter 370   qQuantileNormalizer::Partitioner::Partitioner(ForwardIterator first,
3545 23 Dec 16 peter 371                                                 ForwardIterator last,
1736 19 Jan 09 peter 372                                                 unsigned int N)
1813 20 Feb 09 peter 373     : average_(utility::Vector(N)), quantiles_(utility::Vector(N))
1736 19 Jan 09 peter 374   {
3547 30 Dec 16 peter 375     BOOST_CONCEPT_ASSERT((boost_concepts::ForwardTraversal<ForwardIterator>));
2263 26 May 10 peter 376     BOOST_CONCEPT_ASSERT((utility::DataIterator<ForwardIterator>));
3545 23 Dec 16 peter 377     build(first, last, N,
1801 17 Feb 09 peter 378           typename utility::weighted_iterator_traits<ForwardIterator>::type());
1737 19 Jan 09 peter 379   }
1737 19 Jan 09 peter 380
1737 19 Jan 09 peter 381
1801 17 Feb 09 peter 382   template<typename ForwardIterator>
3545 23 Dec 16 peter 383   void qQuantileNormalizer::Partitioner::build(ForwardIterator first,
3545 23 Dec 16 peter 384                                                ForwardIterator last,
3545 23 Dec 16 peter 385                                                unsigned int N,
1737 19 Jan 09 peter 386                                                utility::unweighted_iterator_tag)
1737 19 Jan 09 peter 387   {
1736 19 Jan 09 peter 388     utility::Vector vec(std::distance(first, last));
1736 19 Jan 09 peter 389     std::copy(first, last, vec.begin());
1736 19 Jan 09 peter 390     std::sort(vec.begin(), vec.end());
1736 19 Jan 09 peter 391     init(vec, N);
1736 19 Jan 09 peter 392   }
1736 19 Jan 09 peter 393
1736 19 Jan 09 peter 394
1801 17 Feb 09 peter 395   template<typename ForwardIterator>
3545 23 Dec 16 peter 396   void qQuantileNormalizer::Partitioner::build(ForwardIterator first,
3545 23 Dec 16 peter 397                                                ForwardIterator last,
3545 23 Dec 16 peter 398                                                unsigned int N,
1738 20 Jan 09 peter 399                                                utility::weighted_iterator_tag)
1738 20 Jan 09 peter 400   {
1954 07 May 09 jari 401     using namespace utility;
3546 29 Dec 16 peter 402     std::vector<DataWeight> vec(first, last);
1738 20 Jan 09 peter 403     std::sort(vec.begin(), vec.end());
1738 20 Jan 09 peter 404     init(vec, N);
1738 20 Jan 09 peter 405   }
1738 20 Jan 09 peter 406
1738 20 Jan 09 peter 407
1708 13 Jan 09 jari 408 }}} // end of namespace normalizer, yat and thep
1708 13 Jan 09 jari 409
1708 13 Jan 09 jari 410 #endif

Generated on Sat Oct 14 02:22:36 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)