distribution.tcc

Go to the documentation of this file.

 
#if !defined(ULTRA_DISTRIBUTION_H)
#  error "Don't include this file directly, include the specific .h instead"
#endif
 
#if !defined(ULTRA_DISTRIBUTION_TCC)
#define      ULTRA_DISTRIBUTION_TCC

template<ArithmeticFloatingType T>
void distribution<T>::clear()
{
  *this = distribution();
}

template<ArithmeticFloatingType T>
bool distribution<T>::empty() const noexcept
{
  return size() == 0;
}

template<ArithmeticFloatingType T>
std::size_t distribution<T>::size() const noexcept
{
  return size_;
}

template<ArithmeticFloatingType T>
T distribution<T>::max() const noexcept
{
  Expects(size());
  return max_;
}

template<ArithmeticFloatingType T>
T distribution<T>::min() const noexcept
{
  Expects(size());
  return min_;
}

template<ArithmeticFloatingType T>
T distribution<T>::mean() const noexcept
{
  Expects(size());
  return mean_;
}

template<ArithmeticFloatingType T>
T distribution<T>::variance() const
{
  Expects(size());
  return m2_ / static_cast<T>(size());
}

template<ArithmeticFloatingType T>
void distribution<T>::add(T val)
{
  using std::isnan;
  if (isnan(val))
    return;
 
  if (!size())
    min_ = max_ = mean_ = val;
  else if (val < min())
    min_ = val;
  else if (val > max())
    max_ = val;
 
  ++size_;
 
  // Values are rounded before being recorded in `seen_`.
  // The histogram (and derived quantities such as entropy) is therefore
  // computed over _rounded bins_, not raw sample values.
  //
  // Rationale:
  // - avoids pathological growth of histogram keys for continuous domains;
  // - treats numerically close values as belonging to the same bin.
  //
  // Consequences:
  // - entropy reflects the distribution of rounded values;
  // - distinct raw values may collapse into the same bin.
  ++seen_[round_to(val)];
 
  update_variance(val);
}
 
template<ArithmeticFloatingType T>
const std::map<T, std::uintmax_t> &distribution<T>::seen() const
{
  return seen_;
}

template<ArithmeticFloatingType T>
double distribution<T>::entropy() const
{
  if (empty())
    return 0.0;
 
  const double c(1.0 / std::log(2.0));
 
  double h(0.0);
  for (const auto &f : seen())  // f.first: fitness, f.second: sightings
  {
    const auto p(static_cast<double>(f.second) / static_cast<double>(size()));
 
    h -= p * std::log(p) * c;
  }
 
  return h;
}

template<ArithmeticFloatingType T>
void distribution<T>::update_variance(T val)
{
  Expects(size());
 
  // Common Kahan summation helper.
  static const auto kahan_add([](T &acc, T &comp, T x)
  {
    // `comp` is zero the first time around.
    const T y(x - comp);
 
    // Alas, `acc` is big, `y` small, so low-order digits of `y` are lost.
    const T t(acc + y);
 
    // `(t - acc)` cancels the high-order part of y;
    // subtracting `y` recovers negative (low part of `y`)
    comp = (t - acc) - y;
    // Algebraically, `comp` should always be zero.
 
    acc = t;
  });
 
  const T delta(val - mean_);
 
  kahan_add(mean_, mean_comp_, delta / static_cast<T>(size_));
  kahan_add(m2_, m2_comp_, delta * (val - mean_));
}

template<ArithmeticFloatingType T>
T distribution<T>::standard_deviation() const
{
  // This way, for "regular" types we'll use std::sqrt ("taken in" by the
  // using statement), while for our types the overload will prevail due to
  // Koenig lookup (http://www.gotw.ca/gotw/030.htm).
  using std::sqrt;
 
  return sqrt(variance());
}

template<ArithmeticFloatingType T>
bool distribution<T>::save(std::ostream &out) const
{
  out << size() << '\n';
 
  if (size())
  {
    SAVE_FLAGS(out);
 
    out << std::scientific
        << std::setprecision(std::numeric_limits<T>::digits10 + 1)
        << mean() << '\n'
        << min() << '\n'
        << max() << '\n'
        << m2_ << '\n';
 
    out << seen().size() << '\n';
    for (const auto &elem : seen())
      out << elem.first << ' ' << elem.second << '\n';
  }
 
  return out.good();
}

template<ArithmeticFloatingType T>
bool distribution<T>::load(std::istream &in)
{
  distribution<T> temp;
 
  if (!(in >> temp.size_))
    return false;
 
  if (!temp.size_)
  {
    clear();
    return true;
  }
 
  SAVE_FLAGS(in);
 
  in >> std::scientific
     >> std::setprecision(std::numeric_limits<T>::digits10 + 1);
 
  if (!(in >> temp.mean_))
    return false;
 
  if (!(in >> temp.min_))
    return false;
 
  if (!(in >> temp.max_))
    return false;
 
  if (!(in >> temp.m2_))
    return false;
 
  typename decltype(seen_)::size_type n;
  if (!(in >> n))
    return false;
 
  for (decltype(n) i(0); i < n; ++i)
  {
    typename decltype(seen_)::key_type key;
    typename decltype(seen_)::mapped_type val;
    if (!(in >> key >> val))
      return false;
 
    temp.seen_[key] = val;
  }
 
  if (!temp.is_valid())
    return false;
 
  *this = std::move(temp);
  return true;
}

template<ArithmeticFloatingType T>
bool distribution<T>::is_valid() const
{
  if (!size())
  {
    distribution<T> empty_dist;
 
    return seen_ == empty_dist.seen_ && issmall(m2_ - empty_dist.m2_)
           && issmall(m2_comp_ - empty_dist.m2_comp_)
           && issmall(mean_ - empty_dist.mean_)
           && issmall(mean_comp_ - empty_dist.mean_comp_)
           && issmall(max_ - empty_dist.max_) && issmall(min_ - empty_dist.min_);
  }
 
  // This way, for "regular" types we'll use `std::infinite` / `std::isnan`
  // ("taken in" by the using statement), while for our types the overload
  // will prevail due to Koenig lookup (http://www.gotw.ca/gotw/030.htm).
  using std::isfinite;
  using std::isnan;
 
  if (isfinite(min()) && isfinite(mean()) && min() > mean())
  {
    ultraERROR("Distribution: min={} > mean={}",
               internal::streamed(min()),
               internal::streamed(mean()));
    return false;
  }
 
  if (isfinite(max()) && isfinite(mean()) && max() < mean())
  {
    ultraERROR("Distribution: max={} < mean={}",
               internal::streamed(max()),
               internal::streamed(mean()));
    return false;
  }
 
  if (isnan(variance()) || !isnonnegative(variance()))
  {
    ultraERROR("Distribution: negative variance");
    return false;
  }
 
  std::size_t total(0);
  for (const auto values(seen_ | std::views::values); auto v : values)
  {
    if (!v)
    {
      ultraERROR("Seen histogram contains zero-count bin");
      return false;
    }
 
    if (v > std::numeric_limits<std::size_t>::max() - total)
    {
      ultraERROR("Seen histogram total overflow");
      return false;
    }
 
    total += v;
  }
 
  if (total != size())
  {
    ultraERROR("Seen histogram total {} does not match distribution size {}",
                 total, size_);
    return false;
  }
 
  return true;
}

template<ArithmeticFloatingType T>
void distribution<T>::merge(const distribution<T> &d2)
{
  if (empty())
  {
    *this = d2;
    return;
  }
 
  if (!d2.size())
    return;
 
  const auto max1(max());
  const auto max2(d2.max());
 
  const auto min1(min());
  const auto min2(d2.min());
 
  const auto size1(size());
  const auto size2(d2.size());
  const auto new_size(size1 + size2);
 
  const auto n1(static_cast<T>(size1));
  const auto n2(static_cast<T>(size2));
  const auto n(static_cast<T>(new_size));
 
  const auto delta(d2.mean() - mean());
  const auto new_mean(size1 >= size2
    ? mean() + delta * (n2 / n)
    : d2.mean() - delta * (n1 / n));
 
  m2_ += d2.m2_ + delta * delta * (n1 * n2 / n);
  mean_ = new_mean;
 
  mean_comp_ = T{};
  m2_comp_ = T{};
 
  max_ = std::max(max1, max2);
  min_ = std::min(min1, min2);
 
  size_ = new_size;
 
  for (auto &[k, v] : d2.seen_)
    seen_[k] += v;
}
 
#endif  // include guard