libsequence: Correlations.tcc Source File

00001 // -*- C++ -*-
00002 /*
00003 
00004 Copyright (C) 2003-2009 Kevin Thornton, krthornt[]@[]uci.edu
00005 
00006 Remove the brackets to email me.
00007 
00008 This file is part of libsequence.
00009 
00010 libsequence is free software: you can redistribute it and/or modify
00011 it under the terms of the GNU General Public License as published by
00012 the Free Software Foundation, either version 3 of the License, or
00013 (at your option) any later version.
00014 
00015 libsequence is distributed in the hope that it will be useful,
00016 but WITHOUT ANY WARRANTY; without even the implied warranty of
00017 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018 GNU General Public License for more details.
00019 
00020 You should have received a copy of the GNU General Public License
00021 long with libsequence.  If not, see <http://www.gnu.org/licenses/>.
00022 
00023 */
00024 
00025 
00030 #ifndef __CORRELATIONS_TCC__
00031 #define __CORRELATIONS_TCC__
00032 
00033 #include <Sequence/Portability/random_shuffle.hpp>
00034 #include <algorithm>
00035 #include <vector>
00036 #include <map>
00037 #include <cmath>
00038 #include <limits>
00039 
00040 namespace Sequence 
00041 {
00052   template<typename iter1, typename iter2>
00053   typename ensureFloating<typename std::iterator_traits<iter1>::value_type,
00054                           typename std::iterator_traits<iter2>::value_type>::type
00055   ProductMoment::operator()(iter1 beg_x, iter1 end_x, iter2 beg_y) const
00056   {
00057     typedef typename ensureFloating<typename std::iterator_traits<iter1>::value_type,
00058       typename std::iterator_traits<iter2>::value_type>::type rtype;
00059 
00060     if (beg_x >= end_x) return std::numeric_limits<rtype>::min();
00061 
00062     unsigned nsam = 0;
00063     rtype x_bar=rtype();
00064     rtype y_bar=rtype();
00065     rtype x_squared=rtype();
00066     rtype y_squared=rtype();
00067     rtype xy=rtype();
00068     
00069     for ( ; beg_x != end_x  ; ++beg_x,++beg_y)
00070       {
00071         x_bar += *beg_x;
00072         x_squared += (*beg_x)*(*beg_x);
00073         y_bar += *beg_y;
00074         y_squared += (*beg_y)*(*beg_y); 
00075         xy += (*beg_x)*(*beg_y);
00076         ++nsam;
00077       }
00078     x_bar /= rtype(nsam);
00079     y_bar /= rtype(nsam);
00080     x_squared /= rtype(nsam);
00081     y_squared /= rtype(nsam);
00082     xy -= x_bar*rtype(nsam)*y_bar;
00083     xy -= y_bar*rtype(nsam)*x_bar;
00084     xy += rtype(nsam)*x_bar*y_bar;
00085     rtype sum_x_sq = rtype(nsam)*(x_squared-x_bar*x_bar);
00086     rtype sum_y_sq = rtype(nsam)*(y_squared-y_bar*y_bar);
00087     return(xy/std::pow((sum_x_sq*sum_y_sq),0.5));
00088   }
00089 
00090   template<typename iter1, typename iter2>
00091   typename ensureFloating<typename std::iterator_traits<iter1>::value_type,
00092                           typename std::iterator_traits<iter2>::value_type>::type
00093   SpearmansRank::operator()(iter1 beg_x, iter1 end_x, iter2 beg_y) const
00094   {
00095     typedef typename ensureFloating<typename std::iterator_traits<iter1>::value_type,
00096       typename std::iterator_traits<iter2>::value_type>::type rtype;
00097   
00098   if(beg_x>=end_x) return std::numeric_limits<rtype>::min();
00099   
00100   typedef typename std::iterator_traits<iter1>::value_type t1;
00101   typedef typename std::iterator_traits<iter2>::value_type t2;
00102   typedef typename std::iterator_traits<iter1>::difference_type dtype1;
00103   typedef typename std::iterator_traits<iter2>::difference_type dtype2;
00104   dtype1 d = end_x-beg_x;
00105   std::vector<t1> xranks;
00106   std::vector<t2> yranks;
00107   //get the unique elements in each range
00108   std::unique_copy(beg_x,end_x,std::back_inserter(xranks));
00109   std::unique_copy(beg_y,beg_y+size_t(d),std::back_inserter(yranks));
00110   //sort them
00111   std::sort(xranks.begin(),xranks.end());
00112   std::sort(yranks.begin(),yranks.end());
00113 
00114   //Recode the data using a map<value,rank>,
00115   //as it's a quick way to do the rank lookups
00116   std::map<t1,unsigned> mx;
00117   std::map<t2,unsigned> my;
00118   unsigned rank=0;
00119   for (typename std::vector<t1>::const_iterator i =
00120          xranks.begin() ;
00121        i != xranks.end() ;
00122        ++i)
00123     {
00124       mx[*i] = rank++;
00125     }
00126   rank=0;
00127   for (typename std::vector<t2>::const_iterator i =
00128          yranks.begin() ;
00129        i != yranks.end() ;
00130        ++i)
00131     {
00132       my[*i] = rank++;
00133     }
00134   //calculate Pearson's rho on the ranks
00135   unsigned nsam = 0;
00136   rtype x_bar=rtype();
00137   rtype y_bar=rtype();
00138   rtype x_squared=rtype();
00139   rtype y_squared=rtype();
00140   rtype xy=rtype();
00141   typename std::vector<t1>::const_iterator a;
00142   typename std::vector<t2>::const_iterator b;
00143   rtype ra,rb;
00144   for( ; beg_x != end_x ; ++beg_x,++beg_y )
00145     {
00146       //find where each element is in the ranks vector
00147       ra = mx[*beg_x];
00148       rb = my[*beg_y];
00149       x_bar += ra;
00150       x_squared += ra*ra;
00151       y_bar += rb;
00152       y_squared += rb*rb;
00153       xy += ra*rb;
00154       ++nsam;
00155     }
00156     x_bar /= rtype(nsam);
00157     y_bar /= rtype(nsam);
00158     x_squared /= rtype(nsam);
00159     y_squared /= rtype(nsam);
00160     xy -= x_bar*rtype(nsam)*y_bar;
00161     xy -= y_bar*rtype(nsam)*x_bar;
00162     xy += rtype(nsam)*x_bar*y_bar;
00163     rtype sum_x_sq = rtype(nsam)*(x_squared-x_bar*x_bar);
00164     rtype sum_y_sq = rtype(nsam)*(y_squared-y_bar*y_bar);
00165     return(xy/std::pow((sum_x_sq*sum_y_sq),0.5));
00166   }
00167    
00168   template<typename iter1, typename iter2,  
00169            typename correlation_type,
00170            typename comparison_function,
00171            typename UniformIntGenerator>
00172   typename ensureFloating<typename std::iterator_traits<iter1>::value_type,
00173                           typename std::iterator_traits<iter2>::value_type>::type
00174   PermuteCorrelation_details(iter1 beg_x, iter1 end_x, iter2 beg_y,
00175                              const correlation_type &  corr,
00176                              const comparison_function & comp,
00177                              UniformIntGenerator & rand,
00178                              const unsigned & NPERM)
00179   {
00180     typedef typename  ensureFloating<typename std::iterator_traits<iter1>::value_type,
00181       typename std::iterator_traits<iter2>::value_type>::type rtype;
00182     typedef typename std::iterator_traits<iter1>::value_type type1;
00183 
00184     rtype _obs = corr(beg_x,end_x,beg_y);
00185     unsigned _prob=0;
00186 
00187     type1 *copy_x = new type1[end_x-beg_x+1];
00188     std::copy(beg_x,end_x,copy_x);
00189 
00190     for(unsigned i = 0 ; i < NPERM ; ++i)
00191       {
00192         Sequence::random_shuffle(copy_x,copy_x+(end_x-beg_x),rand);
00193         if ( comp(std::fabs(corr(copy_x,copy_x+(end_x-beg_x),beg_y)),
00194                   std::fabs(_obs)) == true )
00195           {
00196             ++_prob;
00197           }
00198       }
00199     delete [] copy_x;
00200     return rtype(_prob)/rtype(NPERM);
00201   }
00202 
00203   template<typename iter1, typename iter2,  
00204            typename correlation_type,
00205            typename comparison_function,
00206            typename UniformIntGenerator>
00207   typename ensureFloating<typename std::iterator_traits<iter1>::value_type,
00208                           typename std::iterator_traits<iter2>::value_type>::type
00209   PermuteCorrelation(iter1 beg_x, iter1 end_x, iter2 beg_y,
00210                      const correlation_type &  corr,
00211                      const comparison_function & comp,
00212                      UniformIntGenerator & rand,
00213                      const unsigned & NPERM)
00228   {
00229     return PermuteCorrelation_details(beg_x,end_x,beg_y,corr,comp,rand,NPERM);
00230   }
00231 
00232   template<typename iter1, typename iter2,  
00233            typename correlation_type,
00234            typename comparison_function,
00235            typename UniformIntGenerator>
00236   typename ensureFloating<typename std::iterator_traits<iter1>::value_type,
00237                           typename std::iterator_traits<iter2>::value_type>::type
00238   PermuteCorrelation(iter1 beg_x, iter1 end_x, iter2 beg_y,
00239                      const correlation_type &  corr,
00240                      const comparison_function & comp,
00241                      const UniformIntGenerator & rand,
00242                      const unsigned & NPERM)
00257   {
00258     return PermuteCorrelation_details(beg_x,end_x,beg_y,corr,comp,rand,NPERM);
00259   }
00260 
00261 }//namespace Sequence
00262 #endif