SourceXtractorPlusPlus  0.15
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
KappaSigmaBinning.h
Go to the documentation of this file.
1 
18 #ifndef SOURCEXTRACTORPLUSPLUS_KAPPASIGMABINNING_H
19 #define SOURCEXTRACTORPLUSPLUS_KAPPASIGMABINNING_H
20 
21 #include <iostream>
22 #include "Histogram/Histogram.h"
23 
24 namespace SourceXtractor {
25 
39 template <typename VarType>
41 public:
53  KappaSigmaBinning(float kappa1 = 2., float kappa2 = 5., size_t min_pixels = 4, size_t max_size = 4096)
54  : m_kappa(kappa1), m_kappa2(kappa2), m_min_pixels(min_pixels), m_max_size(max_size),
55  m_scale(1), m_zero(0), m_const(0) {}
56 
68  template<typename Iterator>
69  void computeBins(Iterator begin, Iterator end) {
70  // Compute mean and standard deviation of the original data set
71  double mean, sigma;
72  size_t ndata;
73  Stats stats;
74  for (auto i = begin; i != end; ++i) {
75  stats(*i);
76  }
77  std::tie(mean, sigma, ndata) = stats.get();
78 
79  if (sigma == 0) {
80  sigma = std::abs(*begin - mean);
81  }
82  else {
83  sigma *= m_kappa;
84  }
85 
86  // Cuts
87  auto lcut = mean - sigma;
88  auto hcut = mean + sigma;
89 
90  // Re-compute mean and standard deviation of values within cut
91  stats.reset();
92  for (auto i = begin; i != end; ++i) {
93  if (*i >= lcut && *i <= hcut)
94  stats(*i);
95  }
96  std::tie(mean, sigma, ndata) = stats.get();
97 
98  assert(ndata > 0);
99 
100  // Number of bins
101  m_nbins = computeBinCount(ndata);
102 
103  // Bin size and offset
104  m_scale = 2 * (m_kappa2 * sigma) / m_nbins;
105  if (m_scale == 0)
106  m_scale = 1;
107  m_zero = mean - (m_kappa2 * sigma);
108  m_const = 0.49999 - m_zero / m_scale;
109  }
110 
111  ssize_t getBinIndex(VarType value) const final {
112  return value / m_scale + m_const;
113  }
114 
115  VarType getEdge(size_t e) const final {
116  return (static_cast<float>(e) - 0.5) * m_scale + m_zero;
117  }
118 
119  VarType getBin(size_t i) const final {
120  return i * m_scale + m_zero;
121  }
122 
123 private:
124  VarType m_kappa, m_kappa2;
127  VarType m_scale, m_zero, m_const;
128 
129  size_t computeBinCount(size_t ndata) const {
130  size_t nbins = ndata * std::sqrt(M_2_PI) * m_kappa2 / m_min_pixels + 1;
131  return std::min(nbins, static_cast<size_t>(4096));
132  }
133 
138  struct Stats {
139  double mean = 0, sigma = 0;
140  size_t ndata = 0;
141 
142  void operator() (VarType v) {
143  mean += v;
144  sigma += v * v;
145  ++ndata;
146  }
147 
149  mean /= ndata;
150  sigma = sigma / ndata - mean * mean;
152  sigma = std::sqrt(sigma);
153  else
154  sigma = 0.;
155  return std::make_tuple(mean, sigma, ndata);
156  }
157 
158  void reset() {
159  mean = sigma = 0;
160  ndata = 0;
161  }
162  };
163 };
164 
165 } // end of namespace SourceXtractor
166 
167 #endif //SOURCEXTRACTORPLUSPLUS_KAPPASIGMABINNING_H
T tie(T...args)
constexpr double e
VarType getEdge(size_t e) const final
VarType getBin(size_t i) const final
T make_tuple(T...args)
size_t computeBinCount(size_t ndata) const
ssize_t getBinIndex(VarType value) const final
T min(T...args)
std::tuple< double, double, size_t > get()
T sqrt(T...args)
KappaSigmaBinning(float kappa1=2., float kappa2=5., size_t min_pixels=4, size_t max_size=4096)
void computeBins(Iterator begin, Iterator end)