yat/normalizer/Spearman.h

Code
Comments
Other
Rev Date Author Line
1497 12 Sep 08 peter 1 #ifndef _theplu_yat_normalizer_spearman_
1497 12 Sep 08 peter 2 #define _theplu_yat_normalizer_spearman_
1496 12 Sep 08 peter 3
1575 14 Oct 08 jari 4 // $Id$
1575 14 Oct 08 jari 5
1496 12 Sep 08 peter 6 /*
2119 12 Dec 09 peter 7   Copyright (C) 2008 Jari H√§kkinen, Peter Johansson
3550 03 Jan 17 peter 8   Copyright (C) 2009, 2010, 2014, 2016 Peter Johansson
1496 12 Sep 08 peter 9
1496 12 Sep 08 peter 10   This file is part of the yat library, http://dev.thep.lu.se/yat
1496 12 Sep 08 peter 11
1496 12 Sep 08 peter 12   The yat library is free software; you can redistribute it and/or
1496 12 Sep 08 peter 13   modify it under the terms of the GNU General Public License as
1496 12 Sep 08 peter 14   published by the Free Software Foundation; either version 3 of the
1496 12 Sep 08 peter 15   License, or (at your option) any later version.
1496 12 Sep 08 peter 16
1496 12 Sep 08 peter 17   The yat library is distributed in the hope that it will be useful,
1496 12 Sep 08 peter 18   but WITHOUT ANY WARRANTY; without even the implied warranty of
1496 12 Sep 08 peter 19   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1496 12 Sep 08 peter 20   General Public License for more details.
1496 12 Sep 08 peter 21
1496 12 Sep 08 peter 22   You should have received a copy of the GNU General Public License
1496 12 Sep 08 peter 23   along with yat. If not, see <http://www.gnu.org/licenses/>.
1496 12 Sep 08 peter 24 */
1496 12 Sep 08 peter 25
2155 17 Jan 10 peter 26 #include "utility.h"
2155 17 Jan 10 peter 27
3289 12 Jul 14 peter 28 #include "yat/utility/concept_check.h"
1509 17 Sep 08 peter 29 #include "yat/utility/DataIterator.h"
1496 12 Sep 08 peter 30 #include "yat/utility/iterator_traits.h"
1496 12 Sep 08 peter 31 #include "yat/utility/sort_index.h"
3289 12 Jul 14 peter 32 #include "yat/utility/stl_utility.h"
3289 12 Jul 14 peter 33 #include "yat/utility/utility.h"
1509 17 Sep 08 peter 34 #include "yat/utility/WeightIterator.h"
1496 12 Sep 08 peter 35
2149 17 Jan 10 peter 36 #include <boost/concept_check.hpp>
3543 23 Dec 16 peter 37 #include <boost/iterator/iterator_concepts.hpp>
2149 17 Jan 10 peter 38
1509 17 Sep 08 peter 39 #include <algorithm>
1509 17 Sep 08 peter 40 #include <functional>
3289 12 Jul 14 peter 41 #include <numeric>
1496 12 Sep 08 peter 42 #include <vector>
1496 12 Sep 08 peter 43
1496 12 Sep 08 peter 44 namespace theplu {
1496 12 Sep 08 peter 45 namespace yat {
1497 12 Sep 08 peter 46 namespace normalizer {
1496 12 Sep 08 peter 47
1496 12 Sep 08 peter 48   /**
1496 12 Sep 08 peter 49      \brief Replace elements with normalized rank
1496 12 Sep 08 peter 50
1717 13 Jan 09 peter 51      Each element x is replaced by \f$ \frac{\sum I(x_i-x) w_i}{\sum
1717 13 Jan 09 peter 52      w_i} \f$ where \f$ I(x) = 1 \f$ for \f$ x>0 \f$, \f$I(x) = 0.5
1717 13 Jan 09 peter 53      \f$ for \f$ x=0 \f$, and \f$ I(x) = 0 \f$ for \f$ x<0 \f$.
1717 13 Jan 09 peter 54
1496 12 Sep 08 peter 55      \since New in yat 0.5
1496 12 Sep 08 peter 56    */
1496 12 Sep 08 peter 57   class Spearman
1496 12 Sep 08 peter 58   {
1496 12 Sep 08 peter 59   public:
1496 12 Sep 08 peter 60     /**
1496 12 Sep 08 peter 61        \brief default constructor
1496 12 Sep 08 peter 62      */
1496 12 Sep 08 peter 63     Spearman(void){}
1496 12 Sep 08 peter 64
1496 12 Sep 08 peter 65     /**
1496 12 Sep 08 peter 66        It is possible to centralize a range "in place"; it is
1496 12 Sep 08 peter 67        permissible for the iterators \a first and \a result to be the
1496 12 Sep 08 peter 68        same.
2149 17 Jan 10 peter 69
2149 17 Jan 10 peter 70        Type Requirements:
3289 12 Jul 14 peter 71        - \c RandomAccessIter1 is \random_access_traversal_iterator
3289 12 Jul 14 peter 72        - \c RandomAccessIter1 is \ref concept_data_iterator
3289 12 Jul 14 peter 73        - \c RandomAccessIter2 is \writable_iterator
3289 12 Jul 14 peter 74        - \c RandomAccessIter2 is \random_access_traversal_iterator
3289 12 Jul 14 peter 75        - \c RandomAccessIter2 is \ref concept_data_iterator
1496 12 Sep 08 peter 76      */
2149 17 Jan 10 peter 77     template<typename RandomAccessIter1, typename RandomAccessIter2>
2149 17 Jan 10 peter 78     void operator()(RandomAccessIter1 first, RandomAccessIter1 last,
2149 17 Jan 10 peter 79                     RandomAccessIter2 result) const
1496 12 Sep 08 peter 80     {
3289 12 Jul 14 peter 81       using boost_concepts::RandomAccessTraversal;
3289 12 Jul 14 peter 82       BOOST_CONCEPT_ASSERT((RandomAccessTraversal<RandomAccessIter1>));
3289 12 Jul 14 peter 83       using utility::DataIteratorConcept;
3289 12 Jul 14 peter 84       BOOST_CONCEPT_ASSERT((DataIteratorConcept<RandomAccessIter1>));
3289 12 Jul 14 peter 85
3289 12 Jul 14 peter 86       using boost_concepts::WritableIterator;
3289 12 Jul 14 peter 87       BOOST_CONCEPT_ASSERT((WritableIterator<RandomAccessIter2>));
3289 12 Jul 14 peter 88       BOOST_CONCEPT_ASSERT((RandomAccessTraversal<RandomAccessIter2>));
3289 12 Jul 14 peter 89       BOOST_CONCEPT_ASSERT((DataIteratorConcept<RandomAccessIter2>));
3289 12 Jul 14 peter 90
2155 17 Jan 10 peter 91       using utility::weighted_if_any2;
2155 17 Jan 10 peter 92       typename weighted_if_any2<RandomAccessIter1,RandomAccessIter2>::type tag;
1739 21 Jan 09 peter 93       normalize(first, last, result, tag);
1507 16 Sep 08 peter 94     }
1507 16 Sep 08 peter 95
1507 16 Sep 08 peter 96
1507 16 Sep 08 peter 97   private:
1507 16 Sep 08 peter 98     // unweighted version
2148 17 Jan 10 peter 99     template<typename RandomAccessIter1, typename RandomAccessIter2>
2148 17 Jan 10 peter 100     void normalize(RandomAccessIter1 first, RandomAccessIter1 last,
3289 12 Jul 14 peter 101                    RandomAccessIter2 result,
1739 21 Jan 09 peter 102                    utility::unweighted_iterator_tag) const
1507 16 Sep 08 peter 103     {
1496 12 Sep 08 peter 104       std::vector<size_t> perm;
1496 12 Sep 08 peter 105       utility::sort_index(first, last, perm);
1496 12 Sep 08 peter 106       double n = perm.size();
3289 12 Jul 14 peter 107       size_t i=0;
1512 19 Sep 08 peter 108       while ( i<perm.size() ) {
1512 19 Sep 08 peter 109         size_t min_i = i;
1512 19 Sep 08 peter 110         while (i<perm.size() && first[perm[i]]<=first[perm[min_i]])
1512 19 Sep 08 peter 111            ++i;
3289 12 Jul 14 peter 112         double res = (i + min_i)/(2.0*n);
1512 19 Sep 08 peter 113         for ( ; min_i < i; ++min_i)
1512 19 Sep 08 peter 114           result[perm[min_i]] = res;
1512 19 Sep 08 peter 115       }
1496 12 Sep 08 peter 116     }
1496 12 Sep 08 peter 117
1508 17 Sep 08 peter 118
1508 17 Sep 08 peter 119     // weighted version
2148 17 Jan 10 peter 120     template<typename RandomAccessIter1, typename RandomAccessIter2>
2148 17 Jan 10 peter 121     void normalize(RandomAccessIter1 first, RandomAccessIter1 last,
3289 12 Jul 14 peter 122                    RandomAccessIter2 result,
1739 21 Jan 09 peter 123                    utility::weighted_iterator_tag) const
1508 17 Sep 08 peter 124     {
2155 17 Jan 10 peter 125       detail::copy_weight_if_weighted(first, last, result);
3289 12 Jul 14 peter 126       std::vector<size_t> perm;
3289 12 Jul 14 peter 127       utility::sort_index(first, last, perm);
2148 17 Jan 10 peter 128       utility::iterator_traits<RandomAccessIter1> trait;
2148 17 Jan 10 peter 129       utility::iterator_traits<RandomAccessIter2> rtrait;
1512 19 Sep 08 peter 130
3289 12 Jul 14 peter 131       double total_w = utility::sum_weight(first, last);
1512 19 Sep 08 peter 132       double sum_w=0;
3289 12 Jul 14 peter 133       size_t i=0;
1512 19 Sep 08 peter 134       while ( i<perm.size() ) {
1512 19 Sep 08 peter 135         double w=0;
1512 19 Sep 08 peter 136         size_t min_i = i;
1512 19 Sep 08 peter 137         while (i<perm.size() && (trait.weight(first+perm[i])==0 ||
3289 12 Jul 14 peter 138                                  trait.data(first+perm[i]) <=
1512 19 Sep 08 peter 139                                  trait.data(first+perm[min_i])) ) {
1512 19 Sep 08 peter 140           w += trait.weight(first+perm[i]);
1512 19 Sep 08 peter 141           ++i;
1512 19 Sep 08 peter 142         }
3289 12 Jul 14 peter 143         double res = (sum_w + 0.5*w)/total_w;
1512 19 Sep 08 peter 144         for ( size_t j=min_i; j<i; ++j)
1512 19 Sep 08 peter 145           rtrait.data(result+perm[j]) = res;
1512 19 Sep 08 peter 146         sum_w += w;
3289 12 Jul 14 peter 147       }
1508 17 Sep 08 peter 148     }
1496 12 Sep 08 peter 149   };
1496 12 Sep 08 peter 150
1497 12 Sep 08 peter 151 }}} // end of namespace normalizer, yat and thep
1496 12 Sep 08 peter 152 #endif