libsequence  1.9.5
Snn.tcc
1 // -*- C++ -*-
2 /*
3 
4 Copyright (C) 2003-2009 Kevin Thornton, krthornt[]@[]uci.edu
5 
6 Remove the brackets to email me.
7 
8 This file is part of libsequence.
9 
10 libsequence is free software: you can redistribute it and/or modify
11 it under the terms of the GNU General Public License as published by
12 the Free Software Foundation, either version 3 of the License, or
13 (at your option) any later version.
14 
15 libsequence is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU General Public License for more details.
19 
20 You should have received a copy of the GNU General Public License
21 long with libsequence. If not, see <http://www.gnu.org/licenses/>.
22 
23 */
24 
25 #ifndef __SEQUENCE_BITS_SNN_TCC__
26 #define __SEQUENCE_BITS_SNN_TCC__
27 
28 namespace Sequence
29 {
30  template< typename shuffler >
31  std::pair<double,double>
32  Snn_test_details(const PolyTable & snpTable,
33  const unsigned config[],
34  const size_t & npop,
35  shuffler & s,
36  const unsigned & nperms)
37  {
38  std::vector < std::vector<double> > dkj(snpTable.size(),
39  std::vector<double>(snpTable.size(),0.));
40  for(unsigned i=0;i<snpTable.size()-1;++i)
41  {
42  for(unsigned j=i+1;j<snpTable.size();++j)
43  {
44  dkj[i][j] = double(NumDiffs(snpTable[i],snpTable[j]));
45  }
46  }
47  std::vector <unsigned> __individuals(snpTable.size(),0u);
48  for(unsigned i=0;i<snpTable.size();++i)
49  {
50  __individuals[i]=i;
51  }
52  const double observed = Snn_statistic(&__individuals[0],
53  dkj,
54  config,
55  npop,unsigned(__individuals.size()));
56  unsigned perm = nperms;
57  unsigned pv = 0;
58  while( perm-- )
59  {
60  //std::random_shuffle(__individuals.begin(),__individuals.end(),uni_int);
61  s(__individuals.begin(),__individuals.end());
62  double permuted_stat = Snn_statistic(&__individuals[0],
63  dkj,
64  config,
65  npop,unsigned(__individuals.size()));
66  pv += (permuted_stat >= observed) ? 1 : 0;
67  }
68  return std::make_pair(observed,double(pv)/double(nperms));
69  }
70 
71  template< typename shuffler >
72  std::vector< std::vector<double> >
73  Snn_test_pairwise_details(const PolyTable & snpTable,
74  const unsigned config[],
75  const size_t & npop,
76  shuffler & s,
77  const unsigned & nperms)
78  {
79  std::vector<unsigned> offsets(npop+1);
80  for(unsigned i=1;i<=npop;++i)
81  {
82  offsets[i] = offsets[i-1]+config[i-1];
83  }
84  std::vector < std::vector<double> > dkj(snpTable.size(),
85  std::vector<double>(snpTable.size(),0.));
86  for(unsigned i=0;i<snpTable.size()-1;++i)
87  {
88  for(unsigned j=i+1;j<snpTable.size();++j)
89  {
90  dkj[i][j] = Sequence::NumDiffs(snpTable[i],snpTable[j]);
91  }
92  }
93  std::vector <unsigned> __individuals(snpTable.size(),0u),individuals;
94  unsigned _config[2];
95 
96  for(unsigned i=0;i<snpTable.size();++i)
97  {
98  __individuals[i]=i;
99  }
100  std::vector< std::vector<double> > rv;
101  for(unsigned i=0;i< npop-1;++i)
102  {
103  for(unsigned j=i+1;j<npop;++j)
104  {
105  individuals.assign(__individuals.begin()+offsets[i],
106  __individuals.begin()+offsets[i+1]);
107  individuals.insert(individuals.end(),
108  __individuals.begin()+offsets[j],
109  __individuals.begin()+offsets[j+1]);
110  _config[0] = config[i];
111  _config[1] = config[j];
112  const double observed = Snn_statistic( &individuals[0],dkj,_config,2,
113  _config[0]+_config[1]);
114  unsigned p =0, _nperms = nperms;
115  while(_nperms>0)
116  {
117  //std::random_shuffle(__individuals.begin(),__individuals.end(),uni_int);
118  s(__individuals.begin(),__individuals.end());
119  double perm = Snn_statistic( &individuals[0],dkj,_config,2,
120  _config[0]+_config[1]);
121  if(perm>=observed) ++p;
122  _nperms--;
123  }
124  std::vector<double> result(4);
125  result[0]=double(i+1);
126  result[1]=double(j+1);
127  result[2]=double(observed);
128  result[3]=double(p)/double(nperms);
129  rv.push_back(result);
130  }
131  }
132  return rv;
133  }
134 
135  template< typename shuffler >
136  std::pair<double,double>
137  Snn_test(const PolyTable & snpTable,
138  const unsigned config[],
139  const size_t & npop,
140  shuffler & s,
141  const unsigned & nperms)
142  /*!
143  @brief Conducts a permutation-test of Hudson's Snn (sequence nearest-neighbor) statistic
144  \param snpTable The data on which we wish to perform the test.
145  \param config An array of the sample sizes in each deme.
146  \param npop The number of populations. For example, npop could equal config.size() if config were a vector
147  \param uni_int A random number generator whose operator() takes one argument, n, and returns a value uniformly-distributed on the half-open interval [0,N)
148  \param nperms The number of permutations to do for the test
149  \return A pair of doubles (std::pair<double,double>). the first member of the pair is
150  the observed value of the statistic, and the second member is the estimated p-value
151  \ingroup popgenanalysis
152  */
153  {
154  return Snn_test_details(snpTable,config,npop,s,nperms);
155  }
156 
157  template< typename shuffler >
158  std::vector< std::vector<double> >
159  Snn_test_pairwise(const PolyTable & snpTable,
160  const unsigned config[],
161  const size_t & npop,
162  shuffler & s,
163  const unsigned & nperms)
164  /*!
165  @brief Conducts a permutation-test of Hudson's Snn (sequence nearest-neighbor) statistic,
166  for all pairwise combinations of populations
167  \param snpTable The data on which we wish to perform the test.
168  \param config An array of the sample sizes in each deme.
169  \param npop The number of populations. For example, npop could equal config.size() if config were a vector
170  \param uni_int A random number generator whose operator() takes one argument, n,
171  and returns a value uniformly-distributed on the half-open interval [0,N)
172  \param nperms The number of permutations to do for the test
173  \return A vector of vector<double>. Each vector contains 4 elements, indexed 0 to 3. Elements 0 and 1
174  are the indexes of the i-th and j-th population. Element 2 is the observed Snn between the i-th and j-th
175  population, and element 3 is the p-value estimated by permutation
176  \ingroup popgenanalysis
177  */
178  {
179  return Snn_test_pairwise_details(snpTable,config,npop,s,nperms);
180  }
181 
182  // template< typename shuffler >
183  // std::pair<double,double>
184  // Snn_test(const PolyTable & snpTable,
185  // const unsigned config[],
186  // const size_t & npop,
187  // shuffler & s,
188  // const unsigned & nperms)
189  // /*!
190  // @brief Conducts a permutation-test of Hudson's Snn (sequence nearest-neighbor) statistic
191  // \param snpTable The data on which we wish to perform the test.
192  // \param config An array of the sample sizes in each deme.
193  // \param npop The number of populations. For example, npop could equal config.size() if config were a vector
194  // \param uni_int A random number generator whose operator() takes one argument, n, and returns a value uniformly-distributed on the half-open interval [0,N)
195  // \param nperms The number of permutations to do for the test
196  // \return A pair of doubles (std::pair<double,double>). the first member of the pair is
197  // the observed value of the statistic, and the second member is the estimated p-value
198  // \ingroup popgenanalysis
199  // */
200  // {
201  // return Snn_test_details(snpTable,config,npop,s,nperms);
202  // }
203 
204  // template< typename shuffler >
205  // std::vector< std::vector<double> >
206  // Snn_test_pairwise(const PolyTable & snpTable,
207  // const unsigned config[],
208  // const size_t & npop,
209  // shuffler & s,
210  // const unsigned & nperms)
211  // /*!
212  // @brief Conducts a permutation-test of Hudson's Snn (sequence nearest-neighbor) statistic,
213  // for all pairwise combinations of populations
214  // \param snpTable The data on which we wish to perform the test.
215  // \param config An array of the sample sizes in each deme.
216  // \param npop The number of populations. For example, npop could equal config.size() if config were a vector
217  // \param uni_int A random number generator whose operator() takes one argument, n,
218  // and returns a value uniformly-distributed on the half-open interval [0,N)
219  // \param nperms The number of permutations to do for the test
220  // \return A vector of vector<double>. Each vector contains 4 elements, indexed 0 to 3. Elements 0 and 1
221  // are the indexes of the i-th and j-th population. Element 2 is the observed Snn between the i-th and j-th
222  // population, and element 3 is the p-value estimated by permutation
223  // \ingroup popgenanalysis
224  // */
225  // {
226  // return Snn_test_pairwise_details(snpTable,config,npop,s,nperms);
227  // }
228 }//ns Sequence
229 #endif