libsequence  1.9.5
hprime_faywuh_aggregator.hpp
1 #ifndef SEQUENCE_DETAIL_HPRIME_FAYWUH_AGGREGATOR_HPP
2 #define SEQUENCE_DETAIL_HPRIME_FAYWUH_AGGREGATOR_HPP
3 
4 #include <Sequence/AlleleCountMatrix.hpp>
5 #include <stdexcept>
6 #include <cmath>
7 
8 namespace Sequence
9 {
10  namespace detail
11  {
12  struct stat_is_faywuh
13  {
14  };
15 
16  struct stat_is_hprime
17  {
18  };
19 
20  class hprime_faywuh_row_processor
21  {
22  private:
23  inline double
24  update_stat(const double d, const stat_is_faywuh)
25  {
26  return std::pow(d, 2.0);
27  }
28 
29  inline double
30  update_stat(const double d, const stat_is_hprime)
31  {
32  return d;
33  }
34 
35  inline double
36  denominator(const std::size_t nsam, const stat_is_faywuh)
37  {
38  return 2./static_cast<double>(nsam * (nsam - 1));
39  }
40 
41  inline double
42  denominator(const std::size_t nsam, const stat_is_hprime)
43  {
44  return 1./static_cast<double>(nsam - 1);
45  }
46 
47  template <typename T>
48  inline void
49  stat_details(const Sequence::AlleleCountMatrix &ac,
50  const std::size_t i, const std::size_t refindex,
51  const T t)
52  {
53  unsigned nstates = 0;
54  bool refseen = false;
55  double temp = 0.0;
56  double homozygosity = 0.0;
57  for (std::size_t j = i; j < i + ac.ncol; ++j)
58  {
59  auto ci = ac.counts[j];
60  if (ci > 0)
61  {
62  homozygosity
63  += static_cast<double>(ci * (ci - 1));
64  ++nstates;
65  if (j - i != refindex)
66  {
67  temp += update_stat(ci, t);
68  }
69  else
70  {
71  refseen = true;
72  }
73  }
74  }
75  if (nstates > 2)
76  {
77  throw std::runtime_error(
78  "site has more than one derived state");
79  }
80 
81  if (nstates > 1)
82  {
83  ++S;
84  }
85  if (refseen)
86  {
87  double nnm1
88  = static_cast<double>(ac.nsam * (ac.nsam - 1));
89  pi += 1.0 - homozygosity / nnm1;
90  double x = denominator(ac.nsam, t);
91  theta += temp * x;
92  }
93  }
94 
95  public:
96  unsigned S;
97  double pi, theta;
98  hprime_faywuh_row_processor()
99  : S{ 0 }, pi{ 0.0 }, theta{ 0.0 }
100  {
101  }
102 
103  inline void
105  const std::size_t i, const std::size_t refindex,
106  const stat_is_faywuh dispatch)
107  {
108  stat_details(ac, i, refindex, dispatch);
109  }
110 
111  inline void
113  const std::size_t i, const std::size_t refindex,
114  const stat_is_hprime dispatch)
115  {
116  stat_details(ac, i, refindex, dispatch);
117  }
118  };
119  } // namespace detail
120 } // namespace Sequence
121 
122 #endif
bool operator()(const std::pair< key, value > &l, const std::pair< key, value > &r) const
The namespace in which this library resides.
Matrix representation of allele counts in a VariantMatrix To be constructed.