KappaSigmaBinning.h

Go to the documentation of this file.

 
#ifndef SOURCEXTRACTORPLUSPLUS_KAPPASIGMABINNING_H
#define SOURCEXTRACTORPLUSPLUS_KAPPASIGMABINNING_H
 
#include <iostream>
#include "Histogram/Histogram.h"
 
namespace SourceXtractor {

template <typename VarType>
class KappaSigmaBinning: public Euclid::Histogram::BinStrategy<VarType> {
public:
  KappaSigmaBinning(float kappa1 = 2., float kappa2 = 5., size_t min_pixels = 4, size_t max_size = 4096)
    : m_kappa(kappa1), m_kappa2(kappa2), m_min_pixels(min_pixels), m_max_size(max_size),
      m_scale(1), m_zero(0), m_const(0) {}

  template<typename Iterator>
  void computeBins(Iterator begin, Iterator end) {
    // Compute mean and standard deviation of the original data set
    double mean, sigma;
    size_t ndata;
    Stats stats;
    for (auto i = begin; i != end; ++i) {
      stats(*i);
    }
    std::tie(mean, sigma, ndata) = stats.get();
 
    if (sigma == 0) {
      sigma = std::abs(*begin - mean);
    }
    else {
      sigma *= m_kappa;
    }
 
    // Cuts
    auto lcut = mean - sigma;
    auto hcut = mean + sigma;
 
    // Re-compute mean and standard deviation of values within cut
    stats.reset();
    for (auto i = begin; i != end; ++i) {
      if (*i >= lcut && *i <= hcut)
        stats(*i);
    }
    std::tie(mean, sigma, ndata) = stats.get();
 
    assert(ndata > 0);
 
    // Number of bins
    m_nbins = computeBinCount(ndata);
 
    // Bin size and offset
    m_scale = 2 * (m_kappa2 * sigma) / m_nbins;
    if (m_scale == 0)
      m_scale = 1;
    m_zero = mean - (m_kappa2 * sigma);
    m_const = 0.49999 - m_zero / m_scale;
  }
 
  ssize_t getBinIndex(VarType value) const final {
    return value / m_scale + m_const;
  }
 
  VarType getEdge(size_t e) const final {
    return (static_cast<float>(e) - 0.5) * m_scale + m_zero;
  }
 
  VarType getBin(size_t i) const final {
    return i * m_scale + m_zero;
  }
 
private:
  VarType m_kappa, m_kappa2;
  size_t m_min_pixels, m_max_size;
  using Euclid::Histogram::BinStrategy<VarType>::m_nbins;
  VarType m_scale, m_zero, m_const;
 
  size_t computeBinCount(size_t ndata) const {
    size_t nbins = ndata * std::sqrt(M_2_PI) * m_kappa2 / m_min_pixels + 1;
    return std::min(nbins, static_cast<size_t>(4096));
  }

  struct Stats {
    double mean = 0, sigma = 0;
    size_t ndata = 0;
 
    void operator() (VarType v) {
      mean += v;
      sigma += v * v;
      ++ndata;
    }
 
    std::tuple<double, double, size_t> get() {
      mean /= ndata;
      sigma = sigma / ndata - mean * mean;
      if (sigma > std::numeric_limits<double>::epsilon())
        sigma = std::sqrt(sigma);
      else
        sigma = 0.;
      return std::make_tuple(mean, sigma, ndata);
    }
 
    void reset() {
      mean = sigma = 0;
      ndata = 0;
    }
  };
};
 
} // end of namespace SourceXtractor
 
#endif //SOURCEXTRACTORPLUSPLUS_KAPPASIGMABINNING_H