From 59e728c7cfc86b73f8503227cd2ecb263f6db070 Mon Sep 17 00:00:00 2001
From: Ryan <ryan.elandt@gmail.com>
Date: Sun, 4 Jun 2023 21:48:24 -0400
Subject: [PATCH 1/8] added files for weighted_fraction accumulator

---
 include/boost/histogram/accumulators.hpp      |   1 +
 .../boost/histogram/accumulators/fraction.hpp |  26 ++-
 .../boost/histogram/accumulators/ostream.hpp  |   8 +
 .../accumulators/weighted_fraction.hpp        | 178 ++++++++++++++++++
 include/boost/histogram/fwd.hpp               |   3 +
 .../histogram/utility/wilson_interval.hpp     |  72 +++++++
 test/CMakeLists.txt                           |   3 +-
 test/Jamfile                                  |   1 +
 test/accumulators_weighted_fraction_test.cpp  | 151 +++++++++++++++
 9 files changed, 433 insertions(+), 10 deletions(-)
 create mode 100644 include/boost/histogram/accumulators/weighted_fraction.hpp
 create mode 100644 test/accumulators_weighted_fraction_test.cpp

diff --git a/include/boost/histogram/accumulators.hpp b/include/boost/histogram/accumulators.hpp
index af449d32..d930ad19 100644
--- a/include/boost/histogram/accumulators.hpp
+++ b/include/boost/histogram/accumulators.hpp
@@ -19,6 +19,7 @@
 
 #include <boost/histogram/accumulators/count.hpp>
 #include <boost/histogram/accumulators/fraction.hpp>
+#include <boost/histogram/accumulators/weighted_fraction.hpp>
 #include <boost/histogram/accumulators/mean.hpp>
 #include <boost/histogram/accumulators/sum.hpp>
 #include <boost/histogram/accumulators/weighted_mean.hpp>
diff --git a/include/boost/histogram/accumulators/fraction.hpp b/include/boost/histogram/accumulators/fraction.hpp
index 81fd41c5..53757d25 100644
--- a/include/boost/histogram/accumulators/fraction.hpp
+++ b/include/boost/histogram/accumulators/fraction.hpp
@@ -9,6 +9,7 @@
 
 #include <boost/core/nvp.hpp>
 #include <boost/histogram/fwd.hpp> // for fraction<>
+#include <boost/histogram/weight.hpp>
 #include <boost/histogram/utility/wilson_interval.hpp>
 #include <type_traits> // for std::common_type
 
@@ -36,7 +37,8 @@ class fraction {
   using const_reference = const value_type&;
   using real_type = typename std::conditional<std::is_floating_point<value_type>::value,
                                               value_type, double>::type;
-  using interval_type = typename utility::wilson_interval<real_type>::interval_type;
+  using score_type = typename utility::wilson_interval<real_type>;
+  using interval_type = typename score_type::interval_type;
 
   fraction() noexcept = default;
 
@@ -51,11 +53,14 @@ class fraction {
                  static_cast<value_type>(e.failures())} {}
 
   /// Insert boolean sample x.
-  void operator()(bool x) noexcept {
+  void operator()(bool x) noexcept { operator()(weight(1), x); }
+
+  /// Insert boolean sample x with weight w.
+  void operator()(const weight_type<value_type>& w, bool x) noexcept {
     if (x)
-      ++succ_;
+      succ_ += w.value;
     else
-      ++fail_;
+      fail_ += w.value;
   }
 
   /// Add another accumulator.
@@ -78,19 +83,22 @@ class fraction {
   real_type value() const noexcept { return static_cast<real_type>(succ_) / count(); }
 
   /// Return variance of the success fraction.
-  real_type variance() const noexcept {
+  real_type variance() const noexcept { return variance_for_p_and_n_eff(value(), count()); }
+
+  /// Calculate the variance for a given success fraction and effective number of samples.
+  template <class T>
+  static real_type variance_for_p_and_n_eff(const real_type& p, const T& n_eff ) noexcept {
     // We want to compute Var(p) for p = X / n with Var(X) = n p (1 - p)
     // For Var(X) see
     // https://en.wikipedia.org/wiki/Binomial_distribution#Expected_value_and_variance
     // Error propagation: Var(p) = p'(X)^2 Var(X) = p (1 - p) / n
-    const real_type p = value();
-    return p * (1 - p) / count();
+    return p * (1 - p) / n_eff;
   }
 
   /// Return standard interval with 68.3 % confidence level (Wilson score interval).
   interval_type confidence_interval() const noexcept {
-    return utility::wilson_interval<real_type>()(static_cast<real_type>(successes()),
-                                                 static_cast<real_type>(failures()));
+    return score_type()(static_cast<real_type>(successes()),
+                        static_cast<real_type>(failures()));
   }
 
   bool operator==(const fraction& rhs) const noexcept {
diff --git a/include/boost/histogram/accumulators/ostream.hpp b/include/boost/histogram/accumulators/ostream.hpp
index 46b914c9..a36ae1c6 100644
--- a/include/boost/histogram/accumulators/ostream.hpp
+++ b/include/boost/histogram/accumulators/ostream.hpp
@@ -101,6 +101,14 @@ std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>&
   return detail::handle_nonzero_width(os, x);
 }
 
+template <class CharT, class Traits, class U>
+std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& os,
+                                              const weighted_fraction<U>& x) {
+  if (os.width() == 0)
+    return os << "weighted_fraction(" << x.get_fraction() << ", " << x.wsum2() << ")";
+  return detail::handle_nonzero_width(os, x);
+}
+
 } // namespace accumulators
 } // namespace histogram
 } // namespace boost
diff --git a/include/boost/histogram/accumulators/weighted_fraction.hpp b/include/boost/histogram/accumulators/weighted_fraction.hpp
new file mode 100644
index 00000000..71579dcd
--- /dev/null
+++ b/include/boost/histogram/accumulators/weighted_fraction.hpp
@@ -0,0 +1,178 @@
+// Copyright 2022 Jay Gohil, Hans Dembinski
+//
+// Distributed under the Boost Software License, version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_HISTOGRAM_ACCUMULATORS_WEIGHTED_FRACTION_HPP
+#define BOOST_HISTOGRAM_ACCUMULATORS_WEIGHTED_FRACTION_HPP
+
+#include <boost/core/nvp.hpp>
+#include <boost/histogram/detail/square.hpp>
+#include <boost/histogram/weight.hpp>
+#include <boost/histogram/accumulators/fraction.hpp>
+
+namespace boost {
+namespace histogram {
+namespace accumulators {
+
+namespace internal {
+
+// Accumulates the sum of weights squared.
+template <class ValueType>
+class sum_of_weights_squared {
+public:
+  using value_type = ValueType;
+  using const_reference = const value_type&;
+
+  sum_of_weights_squared() = default;
+
+  // Allow implicit conversion from sum_of_weights_squared<T>
+  template <class T>
+  sum_of_weights_squared(const sum_of_weights_squared<T>& o) noexcept
+      : sum_of_weights_squared(o.sum_of_weights_squared_) {}
+
+  // Initialize to external sum of weights squared.
+  sum_of_weights_squared(const_reference wsum2) noexcept : sum_of_weights_squared_(wsum2) {}
+
+  // Increment by one.
+  sum_of_weights_squared& operator++() {
+    ++sum_of_weights_squared_;
+    return *this;
+  }
+
+  // Increment by weight.
+  sum_of_weights_squared& operator+=(const weight_type<value_type>& w) {
+    sum_of_weights_squared_ += detail::square(w.value);
+    return *this;
+  }
+
+  // Added another sum_of_weights_squared.
+  sum_of_weights_squared& operator+=(const sum_of_weights_squared& rhs) {
+    sum_of_weights_squared_ += rhs.sum_of_weights_squared_;
+    return *this;
+  }
+
+  bool operator==(const sum_of_weights_squared& rhs) const noexcept {
+    return sum_of_weights_squared_ == rhs.sum_of_weights_squared_;
+  }
+
+  bool operator!=(const sum_of_weights_squared& rhs) const noexcept { return !operator==(rhs); }
+
+  // Return sum of weights squared.
+  const_reference value() const noexcept { return sum_of_weights_squared_; }
+
+private:
+  ValueType sum_of_weights_squared_{};
+};
+
+}  // namespace internal
+
+/// Accumulates weighted boolean samples and computes the fraction of true samples.
+template <class ValueType>
+class weighted_fraction {
+public:
+  using value_type = ValueType;
+  using const_reference = const value_type&;
+  using fraction_type = fraction<ValueType>;
+  using real_type = typename fraction_type::real_type;
+  using score_type = typename fraction_type::score_type;
+  using interval_type = typename score_type::interval_type;
+
+  weighted_fraction() noexcept = default;
+
+  /// Initialize to external fraction and sum of weights squared.
+  weighted_fraction(const fraction_type& f, const_reference wsum2) noexcept
+      : f_(f), wsum2_(wsum2) {}
+
+  /// Convert the weighted_fraction class to a different type T.
+  template <class T>
+  operator weighted_fraction<T>() const noexcept {
+    return weighted_fraction<T>(static_cast<fraction<T>>(f_),
+                                static_cast<T>(wsum2_.value()));
+  }
+
+  /// Insert boolean sample x with weight 1.
+  void operator()(bool x) noexcept { operator()(weight(1), x); }
+
+  /// Insert boolean sample x with weight w.
+  void operator()(const weight_type<value_type>& w, bool x) noexcept {
+    f_(w, x);
+    wsum2_ += w;
+  }
+
+  /// Add another weighted_fraction.
+  weighted_fraction& operator+=(const weighted_fraction& rhs) noexcept {
+    f_ += rhs.f_;
+    wsum2_ += rhs.wsum2_;
+    return *this;
+  }
+
+  bool operator==(const weighted_fraction& rhs) const noexcept {
+    return f_ == rhs.f_ && wsum2_ == rhs.wsum2_;
+  }
+
+  bool operator!=(const weighted_fraction& rhs) const noexcept { return !operator==(rhs); }
+
+  /// Return number of boolean samples that were true.
+  const_reference successes() const noexcept { return f_.successes(); }
+
+  /// Return number of boolean samples that were false.
+  const_reference failures() const noexcept { return f_.failures(); }
+
+  /// Return effective number of boolean samples.
+  real_type count() const noexcept {
+    return static_cast<real_type>(detail::square(f_.count())) / wsum2_.value();
+  }
+
+  /// Return success weighted_fraction of boolean samples.
+  real_type value() const noexcept { return f_.value(); }
+
+  /// Return variance of the success weighted_fraction.
+  real_type variance() const noexcept {
+    return fraction_type::variance_for_p_and_n_eff(value(), count());
+  }
+
+  /// Return the sum of weights squared.
+  value_type sum_of_weights_squared() const noexcept { return wsum2_.value(); } 
+
+  /// Return standard interval with 68.3 % confidence level (Wilson score interval).
+  interval_type confidence_interval() const noexcept {
+    const real_type n_eff = count();
+    const real_type correction = score_type::third_order_correction(n_eff);
+
+    return score_type().wilson_solve({
+      .n_eff=n_eff,
+      .p_hat=value(),
+      .correction=correction});
+  }
+
+  /// Return the fraction.
+  const fraction_type& get_fraction() const noexcept { return f_; }
+
+  /// Return the sum of weights squared.
+  const value_type& wsum2() const noexcept { return wsum2_.value(); }
+
+private:
+  fraction_type f_;
+  internal::sum_of_weights_squared<ValueType> wsum2_;
+};
+
+} // namespace accumulators
+} // namespace histogram
+} // namespace boost
+
+#ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
+
+namespace std {
+template <class T, class U>
+/// Specialization for boost::histogram::accumulators::weighted_fraction.
+struct common_type<boost::histogram::accumulators::weighted_fraction<T>,
+                   boost::histogram::accumulators::weighted_fraction<U>> {
+  using type = boost::histogram::accumulators::weighted_fraction<common_type_t<T, U>>;
+};
+} // namespace std
+
+#endif
+
+#endif
diff --git a/include/boost/histogram/fwd.hpp b/include/boost/histogram/fwd.hpp
index 3c48f784..068c8c7f 100644
--- a/include/boost/histogram/fwd.hpp
+++ b/include/boost/histogram/fwd.hpp
@@ -97,6 +97,9 @@ class count;
 template <class ValueType = double>
 class fraction;
 
+template <class ValueType = double>
+class weighted_fraction;
+
 template <class ValueType = double>
 class sum;
 
diff --git a/include/boost/histogram/utility/wilson_interval.hpp b/include/boost/histogram/utility/wilson_interval.hpp
index 98fbf123..e97c0234 100644
--- a/include/boost/histogram/utility/wilson_interval.hpp
+++ b/include/boost/histogram/utility/wilson_interval.hpp
@@ -16,6 +16,17 @@ namespace boost {
 namespace histogram {
 namespace utility {
 
+namespace internal {
+
+template <class T>
+struct wilson_interval_problem {
+  T n_eff;
+  T p_hat;
+  T correction;
+};
+
+}  // namespace internal
+
 /**
   Wilson interval.
 
@@ -75,7 +86,68 @@ class wilson_interval : public binomial_proportion_interval<ValueType> {
     return {t1 - t2, t1 + t2};
   }
 
+  /// Returns the third order correction for n_eff.
+  static value_type third_order_correction(value_type n_eff) noexcept {
+    // The approximate formula reads:
+    //   f(n) = (n³ + n² + 2n + 6) / n³
+    //
+    // Applying the substitution x = 1 / n gives:
+    //   f(n) = 1 + x + 2x² + 6x³
+    //
+    // Using Horner's method to evaluate this polynomial gives:
+    //   f(n) = 1 + x (1 + x (2 + 6x))
+    if (n_eff == 0) return 1;
+    const value_type x = 1 / n_eff;
+    return 1 + x * (1 + x * (2 + 6 * x));
+  }
+
+  /** Computer the confidence interval for the provided problem.
+
+    @param p The problem to solve.
+  */
+  interval_type wilson_solve(const internal::wilson_interval_problem<value_type>& p) const noexcept {
+    // Equation 41 from this paper: https://arxiv.org/abs/2110.00294
+    //      (p̂ - p)² = p (1 - p) (z² f(n) / n)
+    //    n (p̂ - p)² = p (1 - p) z² f(n)                 (avoid floating point error when n = 0)
+    //    np² - 2np̂p + np̂² = pz²f(n) - p²z²f(n)          (expand)
+    //    p²(n + z²f(n)) + p(-2np̂ - z²f(n)) + (np̂²) = 0  (collect terms of p)
+    //
+    // This is a quadratic equation ap² + bp + c = 0 where
+    //    a = n + z²f(n)
+    //    b = -2np̂ - z²f(n)
+    //    c = np̂²
+
+    const value_type& n = p.n_eff;
+    const value_type& p_hat = p.p_hat;
+    const value_type zz_correction = (z_ * z_) * p.correction;
+
+    const value_type a = n + zz_correction;
+    const value_type b = -2 * n * p_hat - zz_correction;
+    const value_type c = n * (p_hat * p_hat);
+
+    return quadratic_roots(a, b, c);
+  }
+
 private:
+  // Finds the roots of the quadratic equation ax² + bx + c = 0.
+  static interval_type quadratic_roots(const value_type& a, const value_type& b, const value_type& c) noexcept {
+    // https://people.csail.mit.edu/bkph/articles/Quadratics.pdf
+    
+    const value_type two_a = 2 * a;
+    const value_type two_c = 2 * c;
+    const value_type sqrt_bb_4ac = std::sqrt(b * b - two_a * two_c);
+
+    if (b >= 0) {
+      const value_type root1 = (-b - sqrt_bb_4ac) / two_a;
+      const value_type root2 = two_c / (-b - sqrt_bb_4ac);
+      return {root1, root2};
+    } else {
+      const value_type root1 = two_c / (-b + sqrt_bb_4ac);
+      const value_type root2 = (-b + sqrt_bb_4ac) / two_a;
+      return {root1, root2};
+    }
+  }
+
   value_type z_;
 };
 
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 337da698..103ccdfa 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -41,6 +41,8 @@ boost_test(TYPE compile-fail SOURCES histogram_fail4.cpp)
 set(BOOST_TEST_LINK_LIBRARIES Boost::histogram Boost::core)
 
 boost_test(TYPE run SOURCES accumulators_count_test.cpp)
+boost_test(TYPE run SOURCES accumulators_fraction_test.cpp)
+boost_test(TYPE run SOURCES accumulators_weighted_fraction_test.cpp)
 boost_test(TYPE run SOURCES accumulators_mean_test.cpp)
 boost_test(TYPE run SOURCES accumulators_sum_test.cpp)
 boost_test(TYPE run SOURCES accumulators_weighted_mean_test.cpp)
@@ -96,7 +98,6 @@ boost_test(TYPE run SOURCES unlimited_storage_test.cpp)
 boost_test(TYPE run SOURCES tools_test.cpp)
 boost_test(TYPE run SOURCES issue_327_test.cpp)
 boost_test(TYPE run SOURCES issue_353_test.cpp)
-boost_test(TYPE run SOURCES accumulators_fraction_test.cpp)
 boost_test(TYPE run SOURCES utility_binomial_proportion_interval_test.cpp)
 boost_test(TYPE run SOURCES utility_wald_interval_test.cpp)
 boost_test(TYPE run SOURCES utility_wilson_interval_test.cpp)
diff --git a/test/Jamfile b/test/Jamfile
index d38b5398..2f6e6551 100644
--- a/test/Jamfile
+++ b/test/Jamfile
@@ -41,6 +41,7 @@ alias odr :
 alias cxx14 :
     [ run accumulators_count_test.cpp ]
     [ run accumulators_fraction_test.cpp ]
+    [ run accumulators_weighted_fraction_test.cpp ]
     [ run accumulators_mean_test.cpp ]
     [ run accumulators_sum_test.cpp : : :
       # make sure sum accumulator works even with -ffast-math and optimizations
diff --git a/test/accumulators_weighted_fraction_test.cpp b/test/accumulators_weighted_fraction_test.cpp
new file mode 100644
index 00000000..e963da95
--- /dev/null
+++ b/test/accumulators_weighted_fraction_test.cpp
@@ -0,0 +1,151 @@
+// Copyright 2015-2018 Hans Dembinski
+//
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <boost/core/lightweight_test.hpp>
+#include <boost/histogram/accumulators/weighted_fraction.hpp>
+#include <boost/histogram/accumulators/ostream.hpp>
+#include <boost/histogram/utility/wilson_interval.hpp>
+#include <cmath>
+#include <limits>
+#include "is_close.hpp"
+#include "str.hpp"
+#include "throw_exception.hpp"
+
+using namespace boost::histogram;
+using namespace std::literals;
+
+template <class T>
+void run_tests() {
+  using type_fw_t = accumulators::weighted_fraction<T>;
+  using type_f_t = accumulators::fraction<T>;
+
+  const double eps = std::numeric_limits<typename type_fw_t::real_type>::epsilon();
+
+  {
+    type_fw_t f;
+    BOOST_TEST_EQ(f.successes(), 0);
+    BOOST_TEST_EQ(f.failures(), 0);
+    BOOST_TEST_EQ(f.wsum2(), 0);
+    BOOST_TEST(std::isnan(f.value()));
+    BOOST_TEST(std::isnan(f.variance()));
+
+    const auto ci = f.confidence_interval();
+    BOOST_TEST(std::isnan(ci.first));
+    BOOST_TEST(std::isnan(ci.second));
+  }
+
+  {
+    type_fw_t a(type_f_t(1, 0), 1);
+    type_fw_t b(type_f_t(0, 1), 1);
+    a += b;
+    BOOST_TEST_EQ(a, type_fw_t(type_f_t(1, 1), 2));
+
+    a(weight(2), true);  // adds 2 trues and 2^2 to sum_of_weights_squared
+    a(weight(3), false);  // adds 3 falses and 3^2 to sum_of_weights_squared
+    BOOST_TEST_EQ(a, type_fw_t(type_f_t(3, 4), 15));
+  }
+
+  {
+    type_fw_t f;
+    BOOST_TEST_EQ(f.successes(), 0);
+    BOOST_TEST_EQ(f.failures(), 0);
+    BOOST_TEST_EQ(f.wsum2(), 0);
+
+    f(true);
+    BOOST_TEST_EQ(f.successes(), 1);
+    BOOST_TEST_EQ(f.failures(), 0);
+    BOOST_TEST_EQ(f.wsum2(), 1);
+    BOOST_TEST_EQ(str(f), "weighted_fraction(fraction(1, 0), 1)"s);
+    f(false);
+    BOOST_TEST_EQ(f.successes(), 1);
+    BOOST_TEST_EQ(f.failures(), 1);
+    BOOST_TEST_EQ(f.wsum2(), 2);
+    BOOST_TEST_EQ(str(f), "weighted_fraction(fraction(1, 1), 2)"s);
+    BOOST_TEST_EQ(str(f, 41, false), "     weighted_fraction(fraction(1, 1), 2)"s);
+    BOOST_TEST_EQ(str(f, 41, true), "weighted_fraction(fraction(1, 1), 2)     "s);
+  }
+
+  {
+    type_fw_t f(type_f_t(3, 1), 4);
+    BOOST_TEST_EQ(f.successes(), 3);
+    BOOST_TEST_EQ(f.failures(), 1);
+    BOOST_TEST_EQ(f.value(), 0.75);
+    BOOST_TEST_IS_CLOSE(f.variance(), 0.75 * (1 - 0.75) / 4, eps);
+    const auto ci = f.confidence_interval();
+
+    BOOST_TEST_EQ(f.count(), 4);
+
+    // const auto expected = utility::wilson_interval<double>()(3, 1);
+    const auto wilson = utility::wilson_interval<double>();
+    const auto expected = wilson.wilson_solve({
+      .n_eff = 4,
+      .p_hat = 0.75,
+      .correction = 1.46875  // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=4
+    });
+
+    BOOST_TEST_IS_CLOSE(ci.first, expected.first, eps);
+    BOOST_TEST_IS_CLOSE(ci.second, expected.second, eps);
+  }
+
+  {
+    type_fw_t f(type_f_t(0, 1), 1);
+    BOOST_TEST_EQ(f.successes(), 0);
+    BOOST_TEST_EQ(f.failures(), 1);
+    BOOST_TEST_EQ(f.value(), 0);
+    BOOST_TEST_EQ(f.variance(), 0);
+    const auto ci = f.confidence_interval();
+
+    const auto wilson = utility::wilson_interval<double>();
+    const auto expected = wilson.wilson_solve({
+      .n_eff = 1,
+      .p_hat = 0,
+      .correction = 10  // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=1
+    });
+
+    BOOST_TEST_IS_CLOSE(ci.first, expected.first, eps);
+    BOOST_TEST_IS_CLOSE(ci.second, expected.second, eps);
+  }
+
+  {
+    type_fw_t f(type_f_t(1, 0), 1);
+    BOOST_TEST_EQ(f.successes(), 1);
+    BOOST_TEST_EQ(f.failures(), 0);
+    BOOST_TEST_EQ(f.value(), 1);
+    BOOST_TEST_EQ(f.variance(), 0);
+    const auto ci = f.confidence_interval();
+    
+    const auto wilson = utility::wilson_interval<double>();
+    const auto expected = wilson.wilson_solve({
+      .n_eff = 1,
+      .p_hat = 1,
+      .correction = 10  // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=1
+    });
+
+    BOOST_TEST_IS_CLOSE(ci.first, expected.first, eps);
+    BOOST_TEST_IS_CLOSE(ci.second, expected.second, eps);
+  }
+}
+
+int main() {
+  run_tests<int>();
+  run_tests<double>();
+  run_tests<float>();
+
+  {
+    using type_f_double = accumulators::fraction<double>;
+    using type_fw_double = accumulators::weighted_fraction<double>;
+    using type_fw_int = accumulators::weighted_fraction<int>;
+
+    type_fw_double fw_double(type_f_double(5, 3), 88);
+    type_fw_int fw_int(fw_double);
+
+    BOOST_TEST_EQ(fw_int.successes(), 5);
+    BOOST_TEST_EQ(fw_int.failures(), 3);
+    BOOST_TEST_EQ(fw_int.wsum2(), 88);
+  }
+
+  return boost::report_errors();
+}

From cae119db8563d93982374635f6bd670f397a6530 Mon Sep 17 00:00:00 2001
From: Ryan <ryan.elandt@gmail.com>
Date: Sun, 4 Jun 2023 21:56:34 -0400
Subject: [PATCH 2/8] rename wsum2 to sum_w2

---
 .../boost/histogram/accumulators/ostream.hpp  |  2 +-
 .../accumulators/weighted_fraction.hpp        | 22 +++++++++----------
 test/accumulators_weighted_fraction_test.cpp  | 10 ++++-----
 3 files changed, 17 insertions(+), 17 deletions(-)

diff --git a/include/boost/histogram/accumulators/ostream.hpp b/include/boost/histogram/accumulators/ostream.hpp
index a36ae1c6..3fbf3e45 100644
--- a/include/boost/histogram/accumulators/ostream.hpp
+++ b/include/boost/histogram/accumulators/ostream.hpp
@@ -105,7 +105,7 @@ template <class CharT, class Traits, class U>
 std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& os,
                                               const weighted_fraction<U>& x) {
   if (os.width() == 0)
-    return os << "weighted_fraction(" << x.get_fraction() << ", " << x.wsum2() << ")";
+    return os << "weighted_fraction(" << x.get_fraction() << ", " << x.sum_w2() << ")";
   return detail::handle_nonzero_width(os, x);
 }
 
diff --git a/include/boost/histogram/accumulators/weighted_fraction.hpp b/include/boost/histogram/accumulators/weighted_fraction.hpp
index 71579dcd..8fa9a723 100644
--- a/include/boost/histogram/accumulators/weighted_fraction.hpp
+++ b/include/boost/histogram/accumulators/weighted_fraction.hpp
@@ -33,7 +33,7 @@ class sum_of_weights_squared {
       : sum_of_weights_squared(o.sum_of_weights_squared_) {}
 
   // Initialize to external sum of weights squared.
-  sum_of_weights_squared(const_reference wsum2) noexcept : sum_of_weights_squared_(wsum2) {}
+  sum_of_weights_squared(const_reference sum_w2) noexcept : sum_of_weights_squared_(sum_w2) {}
 
   // Increment by one.
   sum_of_weights_squared& operator++() {
@@ -82,14 +82,14 @@ class weighted_fraction {
   weighted_fraction() noexcept = default;
 
   /// Initialize to external fraction and sum of weights squared.
-  weighted_fraction(const fraction_type& f, const_reference wsum2) noexcept
-      : f_(f), wsum2_(wsum2) {}
+  weighted_fraction(const fraction_type& f, const_reference sum_w2) noexcept
+      : f_(f), sum_w2_(sum_w2) {}
 
   /// Convert the weighted_fraction class to a different type T.
   template <class T>
   operator weighted_fraction<T>() const noexcept {
     return weighted_fraction<T>(static_cast<fraction<T>>(f_),
-                                static_cast<T>(wsum2_.value()));
+                                static_cast<T>(sum_w2_.value()));
   }
 
   /// Insert boolean sample x with weight 1.
@@ -98,18 +98,18 @@ class weighted_fraction {
   /// Insert boolean sample x with weight w.
   void operator()(const weight_type<value_type>& w, bool x) noexcept {
     f_(w, x);
-    wsum2_ += w;
+    sum_w2_ += w;
   }
 
   /// Add another weighted_fraction.
   weighted_fraction& operator+=(const weighted_fraction& rhs) noexcept {
     f_ += rhs.f_;
-    wsum2_ += rhs.wsum2_;
+    sum_w2_ += rhs.sum_w2_;
     return *this;
   }
 
   bool operator==(const weighted_fraction& rhs) const noexcept {
-    return f_ == rhs.f_ && wsum2_ == rhs.wsum2_;
+    return f_ == rhs.f_ && sum_w2_ == rhs.sum_w2_;
   }
 
   bool operator!=(const weighted_fraction& rhs) const noexcept { return !operator==(rhs); }
@@ -122,7 +122,7 @@ class weighted_fraction {
 
   /// Return effective number of boolean samples.
   real_type count() const noexcept {
-    return static_cast<real_type>(detail::square(f_.count())) / wsum2_.value();
+    return static_cast<real_type>(detail::square(f_.count())) / sum_w2_.value();
   }
 
   /// Return success weighted_fraction of boolean samples.
@@ -134,7 +134,7 @@ class weighted_fraction {
   }
 
   /// Return the sum of weights squared.
-  value_type sum_of_weights_squared() const noexcept { return wsum2_.value(); } 
+  value_type sum_of_weights_squared() const noexcept { return sum_w2_.value(); } 
 
   /// Return standard interval with 68.3 % confidence level (Wilson score interval).
   interval_type confidence_interval() const noexcept {
@@ -151,11 +151,11 @@ class weighted_fraction {
   const fraction_type& get_fraction() const noexcept { return f_; }
 
   /// Return the sum of weights squared.
-  const value_type& wsum2() const noexcept { return wsum2_.value(); }
+  const value_type& sum_w2() const noexcept { return sum_w2_.value(); }
 
 private:
   fraction_type f_;
-  internal::sum_of_weights_squared<ValueType> wsum2_;
+  internal::sum_of_weights_squared<ValueType> sum_w2_;
 };
 
 } // namespace accumulators
diff --git a/test/accumulators_weighted_fraction_test.cpp b/test/accumulators_weighted_fraction_test.cpp
index e963da95..25c6c670 100644
--- a/test/accumulators_weighted_fraction_test.cpp
+++ b/test/accumulators_weighted_fraction_test.cpp
@@ -28,7 +28,7 @@ void run_tests() {
     type_fw_t f;
     BOOST_TEST_EQ(f.successes(), 0);
     BOOST_TEST_EQ(f.failures(), 0);
-    BOOST_TEST_EQ(f.wsum2(), 0);
+    BOOST_TEST_EQ(f.sum_w2(), 0);
     BOOST_TEST(std::isnan(f.value()));
     BOOST_TEST(std::isnan(f.variance()));
 
@@ -52,17 +52,17 @@ void run_tests() {
     type_fw_t f;
     BOOST_TEST_EQ(f.successes(), 0);
     BOOST_TEST_EQ(f.failures(), 0);
-    BOOST_TEST_EQ(f.wsum2(), 0);
+    BOOST_TEST_EQ(f.sum_w2(), 0);
 
     f(true);
     BOOST_TEST_EQ(f.successes(), 1);
     BOOST_TEST_EQ(f.failures(), 0);
-    BOOST_TEST_EQ(f.wsum2(), 1);
+    BOOST_TEST_EQ(f.sum_w2(), 1);
     BOOST_TEST_EQ(str(f), "weighted_fraction(fraction(1, 0), 1)"s);
     f(false);
     BOOST_TEST_EQ(f.successes(), 1);
     BOOST_TEST_EQ(f.failures(), 1);
-    BOOST_TEST_EQ(f.wsum2(), 2);
+    BOOST_TEST_EQ(f.sum_w2(), 2);
     BOOST_TEST_EQ(str(f), "weighted_fraction(fraction(1, 1), 2)"s);
     BOOST_TEST_EQ(str(f, 41, false), "     weighted_fraction(fraction(1, 1), 2)"s);
     BOOST_TEST_EQ(str(f, 41, true), "weighted_fraction(fraction(1, 1), 2)     "s);
@@ -144,7 +144,7 @@ int main() {
 
     BOOST_TEST_EQ(fw_int.successes(), 5);
     BOOST_TEST_EQ(fw_int.failures(), 3);
-    BOOST_TEST_EQ(fw_int.wsum2(), 88);
+    BOOST_TEST_EQ(fw_int.sum_w2(), 88);
   }
 
   return boost::report_errors();

From 6b312ddd6f75f2c7e6e17b4ddaf5f342486b7b3d Mon Sep 17 00:00:00 2001
From: Ryan <ryan.elandt@gmail.com>
Date: Sun, 4 Jun 2023 23:00:36 -0400
Subject: [PATCH 3/8] Do Pre-commit lints.

---
 .../boost/histogram/accumulators/fraction.hpp |  8 ++++---
 .../accumulators/weighted_fraction.hpp        | 23 ++++++++++--------
 .../histogram/utility/wilson_interval.hpp     | 19 +++++++++------
 test/accumulators_weighted_fraction_test.cpp  | 24 +++++++++----------
 4 files changed, 42 insertions(+), 32 deletions(-)

diff --git a/include/boost/histogram/accumulators/fraction.hpp b/include/boost/histogram/accumulators/fraction.hpp
index 53757d25..67ee3cb7 100644
--- a/include/boost/histogram/accumulators/fraction.hpp
+++ b/include/boost/histogram/accumulators/fraction.hpp
@@ -9,8 +9,8 @@
 
 #include <boost/core/nvp.hpp>
 #include <boost/histogram/fwd.hpp> // for fraction<>
-#include <boost/histogram/weight.hpp>
 #include <boost/histogram/utility/wilson_interval.hpp>
+#include <boost/histogram/weight.hpp>
 #include <type_traits> // for std::common_type
 
 namespace boost {
@@ -83,11 +83,13 @@ class fraction {
   real_type value() const noexcept { return static_cast<real_type>(succ_) / count(); }
 
   /// Return variance of the success fraction.
-  real_type variance() const noexcept { return variance_for_p_and_n_eff(value(), count()); }
+  real_type variance() const noexcept {
+    return variance_for_p_and_n_eff(value(), count());
+  }
 
   /// Calculate the variance for a given success fraction and effective number of samples.
   template <class T>
-  static real_type variance_for_p_and_n_eff(const real_type& p, const T& n_eff ) noexcept {
+  static real_type variance_for_p_and_n_eff(const real_type& p, const T& n_eff) noexcept {
     // We want to compute Var(p) for p = X / n with Var(X) = n p (1 - p)
     // For Var(X) see
     // https://en.wikipedia.org/wiki/Binomial_distribution#Expected_value_and_variance
diff --git a/include/boost/histogram/accumulators/weighted_fraction.hpp b/include/boost/histogram/accumulators/weighted_fraction.hpp
index 8fa9a723..09370aa7 100644
--- a/include/boost/histogram/accumulators/weighted_fraction.hpp
+++ b/include/boost/histogram/accumulators/weighted_fraction.hpp
@@ -8,9 +8,9 @@
 #define BOOST_HISTOGRAM_ACCUMULATORS_WEIGHTED_FRACTION_HPP
 
 #include <boost/core/nvp.hpp>
+#include <boost/histogram/accumulators/fraction.hpp>
 #include <boost/histogram/detail/square.hpp>
 #include <boost/histogram/weight.hpp>
-#include <boost/histogram/accumulators/fraction.hpp>
 
 namespace boost {
 namespace histogram {
@@ -33,7 +33,8 @@ class sum_of_weights_squared {
       : sum_of_weights_squared(o.sum_of_weights_squared_) {}
 
   // Initialize to external sum of weights squared.
-  sum_of_weights_squared(const_reference sum_w2) noexcept : sum_of_weights_squared_(sum_w2) {}
+  sum_of_weights_squared(const_reference sum_w2) noexcept
+      : sum_of_weights_squared_(sum_w2) {}
 
   // Increment by one.
   sum_of_weights_squared& operator++() {
@@ -57,7 +58,9 @@ class sum_of_weights_squared {
     return sum_of_weights_squared_ == rhs.sum_of_weights_squared_;
   }
 
-  bool operator!=(const sum_of_weights_squared& rhs) const noexcept { return !operator==(rhs); }
+  bool operator!=(const sum_of_weights_squared& rhs) const noexcept {
+    return !operator==(rhs);
+  }
 
   // Return sum of weights squared.
   const_reference value() const noexcept { return sum_of_weights_squared_; }
@@ -66,7 +69,7 @@ class sum_of_weights_squared {
   ValueType sum_of_weights_squared_{};
 };
 
-}  // namespace internal
+} // namespace internal
 
 /// Accumulates weighted boolean samples and computes the fraction of true samples.
 template <class ValueType>
@@ -112,7 +115,9 @@ class weighted_fraction {
     return f_ == rhs.f_ && sum_w2_ == rhs.sum_w2_;
   }
 
-  bool operator!=(const weighted_fraction& rhs) const noexcept { return !operator==(rhs); }
+  bool operator!=(const weighted_fraction& rhs) const noexcept {
+    return !operator==(rhs);
+  }
 
   /// Return number of boolean samples that were true.
   const_reference successes() const noexcept { return f_.successes(); }
@@ -134,17 +139,15 @@ class weighted_fraction {
   }
 
   /// Return the sum of weights squared.
-  value_type sum_of_weights_squared() const noexcept { return sum_w2_.value(); } 
+  value_type sum_of_weights_squared() const noexcept { return sum_w2_.value(); }
 
   /// Return standard interval with 68.3 % confidence level (Wilson score interval).
   interval_type confidence_interval() const noexcept {
     const real_type n_eff = count();
     const real_type correction = score_type::third_order_correction(n_eff);
 
-    return score_type().wilson_solve({
-      .n_eff=n_eff,
-      .p_hat=value(),
-      .correction=correction});
+    return score_type().wilson_solve(
+        {.n_eff = n_eff, .p_hat = value(), .correction = correction});
   }
 
   /// Return the fraction.
diff --git a/include/boost/histogram/utility/wilson_interval.hpp b/include/boost/histogram/utility/wilson_interval.hpp
index e97c0234..350350b4 100644
--- a/include/boost/histogram/utility/wilson_interval.hpp
+++ b/include/boost/histogram/utility/wilson_interval.hpp
@@ -25,7 +25,7 @@ struct wilson_interval_problem {
   T correction;
 };
 
-}  // namespace internal
+} // namespace internal
 
 /**
   Wilson interval.
@@ -105,12 +105,16 @@ class wilson_interval : public binomial_proportion_interval<ValueType> {
 
     @param p The problem to solve.
   */
-  interval_type wilson_solve(const internal::wilson_interval_problem<value_type>& p) const noexcept {
+  interval_type wilson_solve(
+      const internal::wilson_interval_problem<value_type>& p) const noexcept {
     // Equation 41 from this paper: https://arxiv.org/abs/2110.00294
     //      (p̂ - p)² = p (1 - p) (z² f(n) / n)
-    //    n (p̂ - p)² = p (1 - p) z² f(n)                 (avoid floating point error when n = 0)
-    //    np² - 2np̂p + np̂² = pz²f(n) - p²z²f(n)          (expand)
-    //    p²(n + z²f(n)) + p(-2np̂ - z²f(n)) + (np̂²) = 0  (collect terms of p)
+    // Multiply by n to avoid floating point error when n = 0.
+    //    n (p̂ - p)² = p (1 - p) z² f(n)
+    // Expand.
+    //    np² - 2np̂p + np̂² = pz²f(n) - p²z²f(n)
+    // Collect terms of p.
+    //    p²(n + z²f(n)) + p(-2np̂ - z²f(n)) + (np̂²) = 0
     //
     // This is a quadratic equation ap² + bp + c = 0 where
     //    a = n + z²f(n)
@@ -130,9 +134,10 @@ class wilson_interval : public binomial_proportion_interval<ValueType> {
 
 private:
   // Finds the roots of the quadratic equation ax² + bx + c = 0.
-  static interval_type quadratic_roots(const value_type& a, const value_type& b, const value_type& c) noexcept {
+  static interval_type quadratic_roots(const value_type& a, const value_type& b,
+                                       const value_type& c) noexcept {
     // https://people.csail.mit.edu/bkph/articles/Quadratics.pdf
-    
+
     const value_type two_a = 2 * a;
     const value_type two_c = 2 * c;
     const value_type sqrt_bb_4ac = std::sqrt(b * b - two_a * two_c);
diff --git a/test/accumulators_weighted_fraction_test.cpp b/test/accumulators_weighted_fraction_test.cpp
index 25c6c670..0431e0e7 100644
--- a/test/accumulators_weighted_fraction_test.cpp
+++ b/test/accumulators_weighted_fraction_test.cpp
@@ -5,8 +5,8 @@
 // or copy at http://www.boost.org/LICENSE_1_0.txt)
 
 #include <boost/core/lightweight_test.hpp>
-#include <boost/histogram/accumulators/weighted_fraction.hpp>
 #include <boost/histogram/accumulators/ostream.hpp>
+#include <boost/histogram/accumulators/weighted_fraction.hpp>
 #include <boost/histogram/utility/wilson_interval.hpp>
 #include <cmath>
 #include <limits>
@@ -44,7 +44,7 @@ void run_tests() {
     BOOST_TEST_EQ(a, type_fw_t(type_f_t(1, 1), 2));
 
     a(weight(2), true);  // adds 2 trues and 2^2 to sum_of_weights_squared
-    a(weight(3), false);  // adds 3 falses and 3^2 to sum_of_weights_squared
+    a(weight(3), false); // adds 3 falses and 3^2 to sum_of_weights_squared
     BOOST_TEST_EQ(a, type_fw_t(type_f_t(3, 4), 15));
   }
 
@@ -81,9 +81,9 @@ void run_tests() {
     // const auto expected = utility::wilson_interval<double>()(3, 1);
     const auto wilson = utility::wilson_interval<double>();
     const auto expected = wilson.wilson_solve({
-      .n_eff = 4,
-      .p_hat = 0.75,
-      .correction = 1.46875  // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=4
+        .n_eff = 4,
+        .p_hat = 0.75,
+        .correction = 1.46875 // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=4
     });
 
     BOOST_TEST_IS_CLOSE(ci.first, expected.first, eps);
@@ -100,9 +100,9 @@ void run_tests() {
 
     const auto wilson = utility::wilson_interval<double>();
     const auto expected = wilson.wilson_solve({
-      .n_eff = 1,
-      .p_hat = 0,
-      .correction = 10  // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=1
+        .n_eff = 1,
+        .p_hat = 0,
+        .correction = 10 // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=1
     });
 
     BOOST_TEST_IS_CLOSE(ci.first, expected.first, eps);
@@ -116,12 +116,12 @@ void run_tests() {
     BOOST_TEST_EQ(f.value(), 1);
     BOOST_TEST_EQ(f.variance(), 0);
     const auto ci = f.confidence_interval();
-    
+
     const auto wilson = utility::wilson_interval<double>();
     const auto expected = wilson.wilson_solve({
-      .n_eff = 1,
-      .p_hat = 1,
-      .correction = 10  // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=1
+        .n_eff = 1,
+        .p_hat = 1,
+        .correction = 10 // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=1
     });
 
     BOOST_TEST_IS_CLOSE(ci.first, expected.first, eps);

From cd3333ef2f33c1534700c818216296e8b6f77109 Mon Sep 17 00:00:00 2001
From: Ryan <ryan.elandt@gmail.com>
Date: Mon, 5 Jun 2023 09:10:36 -0400
Subject: [PATCH 4/8] Remove designated initializer

---
 .../accumulators/weighted_fraction.hpp        |  5 ++--
 .../histogram/utility/wilson_interval.hpp     | 26 +++++-----------
 test/accumulators_weighted_fraction_test.cpp  | 30 +++++++++----------
 3 files changed, 24 insertions(+), 37 deletions(-)

diff --git a/include/boost/histogram/accumulators/weighted_fraction.hpp b/include/boost/histogram/accumulators/weighted_fraction.hpp
index 09370aa7..fcd5aeb2 100644
--- a/include/boost/histogram/accumulators/weighted_fraction.hpp
+++ b/include/boost/histogram/accumulators/weighted_fraction.hpp
@@ -144,10 +144,9 @@ class weighted_fraction {
   /// Return standard interval with 68.3 % confidence level (Wilson score interval).
   interval_type confidence_interval() const noexcept {
     const real_type n_eff = count();
+    const real_type p_hat = value();
     const real_type correction = score_type::third_order_correction(n_eff);
-
-    return score_type().wilson_solve(
-        {.n_eff = n_eff, .p_hat = value(), .correction = correction});
+    return score_type().wilson_solve_for_neff_phat_correction(n_eff, p_hat, correction);
   }
 
   /// Return the fraction.
diff --git a/include/boost/histogram/utility/wilson_interval.hpp b/include/boost/histogram/utility/wilson_interval.hpp
index 350350b4..336359ff 100644
--- a/include/boost/histogram/utility/wilson_interval.hpp
+++ b/include/boost/histogram/utility/wilson_interval.hpp
@@ -16,17 +16,6 @@ namespace boost {
 namespace histogram {
 namespace utility {
 
-namespace internal {
-
-template <class T>
-struct wilson_interval_problem {
-  T n_eff;
-  T p_hat;
-  T correction;
-};
-
-} // namespace internal
-
 /**
   Wilson interval.
 
@@ -105,8 +94,9 @@ class wilson_interval : public binomial_proportion_interval<ValueType> {
 
     @param p The problem to solve.
   */
-  interval_type wilson_solve(
-      const internal::wilson_interval_problem<value_type>& p) const noexcept {
+  interval_type wilson_solve_for_neff_phat_correction(
+      const value_type& n_eff, const value_type& p_hat,
+      const value_type& correction) const noexcept {
     // Equation 41 from this paper: https://arxiv.org/abs/2110.00294
     //      (p̂ - p)² = p (1 - p) (z² f(n) / n)
     // Multiply by n to avoid floating point error when n = 0.
@@ -121,13 +111,11 @@ class wilson_interval : public binomial_proportion_interval<ValueType> {
     //    b = -2np̂ - z²f(n)
     //    c = np̂²
 
-    const value_type& n = p.n_eff;
-    const value_type& p_hat = p.p_hat;
-    const value_type zz_correction = (z_ * z_) * p.correction;
+    const value_type zz_correction = (z_ * z_) * correction;
 
-    const value_type a = n + zz_correction;
-    const value_type b = -2 * n * p_hat - zz_correction;
-    const value_type c = n * (p_hat * p_hat);
+    const value_type a = n_eff + zz_correction;
+    const value_type b = -2 * n_eff * p_hat - zz_correction;
+    const value_type c = n_eff * (p_hat * p_hat);
 
     return quadratic_roots(a, b, c);
   }
diff --git a/test/accumulators_weighted_fraction_test.cpp b/test/accumulators_weighted_fraction_test.cpp
index 0431e0e7..02991f1a 100644
--- a/test/accumulators_weighted_fraction_test.cpp
+++ b/test/accumulators_weighted_fraction_test.cpp
@@ -80,11 +80,11 @@ void run_tests() {
 
     // const auto expected = utility::wilson_interval<double>()(3, 1);
     const auto wilson = utility::wilson_interval<double>();
-    const auto expected = wilson.wilson_solve({
-        .n_eff = 4,
-        .p_hat = 0.75,
-        .correction = 1.46875 // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=4
-    });
+    const auto expected = wilson.wilson_solve_for_neff_phat_correction(
+        4,      // n_eff = 4
+        0.75,   // p_hat = 0.75
+        1.46875 // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=4
+    );
 
     BOOST_TEST_IS_CLOSE(ci.first, expected.first, eps);
     BOOST_TEST_IS_CLOSE(ci.second, expected.second, eps);
@@ -99,11 +99,11 @@ void run_tests() {
     const auto ci = f.confidence_interval();
 
     const auto wilson = utility::wilson_interval<double>();
-    const auto expected = wilson.wilson_solve({
-        .n_eff = 1,
-        .p_hat = 0,
-        .correction = 10 // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=1
-    });
+    const auto expected = wilson.wilson_solve_for_neff_phat_correction(
+        1, // n_eff = 1
+        0, // p_hat = 0
+        10 // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=1
+    );
 
     BOOST_TEST_IS_CLOSE(ci.first, expected.first, eps);
     BOOST_TEST_IS_CLOSE(ci.second, expected.second, eps);
@@ -118,11 +118,11 @@ void run_tests() {
     const auto ci = f.confidence_interval();
 
     const auto wilson = utility::wilson_interval<double>();
-    const auto expected = wilson.wilson_solve({
-        .n_eff = 1,
-        .p_hat = 1,
-        .correction = 10 // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=1
-    });
+    const auto expected = wilson.wilson_solve_for_neff_phat_correction(
+        1, // n_eff = 1
+        1, // p_hat = 1
+        10 // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=1
+    );
 
     BOOST_TEST_IS_CLOSE(ci.first, expected.first, eps);
     BOOST_TEST_IS_CLOSE(ci.second, expected.second, eps);

From 4a798ffbc55f9145f794966f73c87353b493765d Mon Sep 17 00:00:00 2001
From: Ryan <ryan.elandt@gmail.com>
Date: Mon, 5 Jun 2023 12:29:29 -0400
Subject: [PATCH 5/8] Added serialization functions and tests

---
 .../accumulators/weighted_fraction.hpp        | 13 +++++++++
 include/boost/histogram/fwd.hpp               |  7 +++++
 test/accumulators_serialization_test.cpp      | 29 +++++++++++++++++++
 ...ialization_test_sum_of_weights_squared.xml | 17 +++++++++++
 ...s_serialization_test_weighted_fraction.xml | 24 +++++++++++++++
 5 files changed, 90 insertions(+)
 create mode 100644 test/accumulators_serialization_test_sum_of_weights_squared.xml
 create mode 100644 test/accumulators_serialization_test_weighted_fraction.xml

diff --git a/include/boost/histogram/accumulators/weighted_fraction.hpp b/include/boost/histogram/accumulators/weighted_fraction.hpp
index fcd5aeb2..5fc1e495 100644
--- a/include/boost/histogram/accumulators/weighted_fraction.hpp
+++ b/include/boost/histogram/accumulators/weighted_fraction.hpp
@@ -10,7 +10,9 @@
 #include <boost/core/nvp.hpp>
 #include <boost/histogram/accumulators/fraction.hpp>
 #include <boost/histogram/detail/square.hpp>
+#include <boost/histogram/fwd.hpp> // for weighted_fraction<>
 #include <boost/histogram/weight.hpp>
+#include <type_traits> // for std::common_type
 
 namespace boost {
 namespace histogram {
@@ -65,6 +67,11 @@ class sum_of_weights_squared {
   // Return sum of weights squared.
   const_reference value() const noexcept { return sum_of_weights_squared_; }
 
+  template <class Archive>
+  void serialize(Archive& ar, unsigned /* version */) {
+    ar& make_nvp("sum_of_weights_squared", sum_of_weights_squared_);
+  }
+
 private:
   ValueType sum_of_weights_squared_{};
 };
@@ -155,6 +162,12 @@ class weighted_fraction {
   /// Return the sum of weights squared.
   const value_type& sum_w2() const noexcept { return sum_w2_.value(); }
 
+  template <class Archive>
+  void serialize(Archive& ar, unsigned /* version */) {
+    ar& make_nvp("fraction", f_);
+    ar& make_nvp("sum_of_weights_squared", sum_w2_);
+  }
+
 private:
   fraction_type f_;
   internal::sum_of_weights_squared<ValueType> sum_w2_;
diff --git a/include/boost/histogram/fwd.hpp b/include/boost/histogram/fwd.hpp
index 068c8c7f..edc88e55 100644
--- a/include/boost/histogram/fwd.hpp
+++ b/include/boost/histogram/fwd.hpp
@@ -97,6 +97,13 @@ class count;
 template <class ValueType = double>
 class fraction;
 
+namespace internal {
+
+template <class ValueType = double>
+class sum_of_weights_squared;
+
+} // namespace internal
+
 template <class ValueType = double>
 class weighted_fraction;
 
diff --git a/test/accumulators_serialization_test.cpp b/test/accumulators_serialization_test.cpp
index b09132a8..c42797e0 100644
--- a/test/accumulators_serialization_test.cpp
+++ b/test/accumulators_serialization_test.cpp
@@ -102,5 +102,34 @@ int main(int argc, char** argv) {
     BOOST_TEST(a == b);
   }
 
+  // sum_of_weights_squared
+  {
+    const auto filename =
+        join(argv[1], "accumulators_serialization_test_sum_of_weights_squared.xml");
+    accumulators::internal::sum_of_weights_squared<> a;
+    ++a;
+    print_xml(filename, a);
+
+    accumulators::internal::sum_of_weights_squared<> b;
+    BOOST_TEST_NOT(a == b);
+    load_xml(filename, b);
+    BOOST_TEST(a == b);
+  }
+
+  // weighted_fraction
+  {
+    const auto filename =
+        join(argv[1], "accumulators_serialization_test_weighted_fraction.xml");
+    accumulators::weighted_fraction<> a;
+    a(true);
+    a(weight(6), false);
+    print_xml(filename, a);
+
+    accumulators::weighted_fraction<> b;
+    BOOST_TEST_NOT(a == b);
+    load_xml(filename, b);
+    BOOST_TEST(a == b);
+  }
+
   return boost::report_errors();
 }
diff --git a/test/accumulators_serialization_test_sum_of_weights_squared.xml b/test/accumulators_serialization_test_sum_of_weights_squared.xml
new file mode 100644
index 00000000..dd0c854b
--- /dev/null
+++ b/test/accumulators_serialization_test_sum_of_weights_squared.xml
@@ -0,0 +1,17 @@
+<!--
+	Copyright 2022 Hans Dembinski
+
+	Distributed under the Boost Software License, Version 1.0.
+	(See accompanying file LICENSE_1_0.txt
+	or copy at http://www.boost.org/LICENSE_1_0.txt)
+-->
+
+<?xml version="1.0" encoding="UTF-8" standalone="yes" ?>
+<!DOCTYPE boost_serialization>
+<boost_serialization signature="serialization::archive" version="19">
+<item class_id="0" tracking_level="0" version="0">
+	<sum_of_weights_squared>1.00000000000000000e+00</sum_of_weights_squared>
+</item>
+</boost_serialization>
+
+
diff --git a/test/accumulators_serialization_test_weighted_fraction.xml b/test/accumulators_serialization_test_weighted_fraction.xml
new file mode 100644
index 00000000..4a5ff2f4
--- /dev/null
+++ b/test/accumulators_serialization_test_weighted_fraction.xml
@@ -0,0 +1,24 @@
+<!--
+	Copyright 2022 Hans Dembinski
+
+	Distributed under the Boost Software License, Version 1.0.
+	(See accompanying file LICENSE_1_0.txt
+	or copy at http://www.boost.org/LICENSE_1_0.txt)
+-->
+
+<?xml version="1.0" encoding="UTF-8" standalone="yes" ?>
+<!DOCTYPE boost_serialization>
+<boost_serialization signature="serialization::archive" version="19">
+<item class_id="0" tracking_level="0" version="0">
+	<fraction class_id="1" tracking_level="0" version="0">
+		<successes>1.00000000000000000e+00</successes>
+		<failures>6.00000000000000000e+00</failures>
+	</fraction>
+	<sum_of_weights_squared class_id="2" tracking_level="0" version="0">
+		<sum_of_weights_squared>3.70000000000000000e+01</sum_of_weights_squared>
+	</sum_of_weights_squared>
+</item>
+</boost_serialization>
+
+
+

From 867e5106e15031281f7d4dcb8a6348ba9a59ea2e Mon Sep 17 00:00:00 2001
From: Ryan <ryan.elandt@gmail.com>
Date: Tue, 6 Jun 2023 11:35:02 -0400
Subject: [PATCH 6/8] Make scoretype not publicly visible

---
 include/boost/histogram/accumulators/fraction.hpp  | 14 ++++++++++----
 .../histogram/accumulators/weighted_fraction.hpp   | 13 +++++++++----
 .../boost/histogram/utility/wilson_interval.hpp    |  2 +-
 test/accumulators_weighted_fraction_test.cpp       |  6 +++---
 4 files changed, 23 insertions(+), 12 deletions(-)

diff --git a/include/boost/histogram/accumulators/fraction.hpp b/include/boost/histogram/accumulators/fraction.hpp
index 67ee3cb7..2246cdc1 100644
--- a/include/boost/histogram/accumulators/fraction.hpp
+++ b/include/boost/histogram/accumulators/fraction.hpp
@@ -9,6 +9,7 @@
 
 #include <boost/core/nvp.hpp>
 #include <boost/histogram/fwd.hpp> // for fraction<>
+#include <boost/histogram/utility/binomial_proportion_interval.hpp>
 #include <boost/histogram/utility/wilson_interval.hpp>
 #include <boost/histogram/weight.hpp>
 #include <type_traits> // for std::common_type
@@ -37,8 +38,8 @@ class fraction {
   using const_reference = const value_type&;
   using real_type = typename std::conditional<std::is_floating_point<value_type>::value,
                                               value_type, double>::type;
-  using score_type = typename utility::wilson_interval<real_type>;
-  using interval_type = typename score_type::interval_type;
+  using interval_type =
+      typename utility::binomial_proportion_interval<real_type>::interval_type;
 
   fraction() noexcept = default;
 
@@ -99,8 +100,13 @@ class fraction {
 
   /// Return standard interval with 68.3 % confidence level (Wilson score interval).
   interval_type confidence_interval() const noexcept {
-    return score_type()(static_cast<real_type>(successes()),
-                        static_cast<real_type>(failures()));
+    return confidence_interval(utility::wilson_interval<real_type>());
+  }
+
+  /// Return interval for the given binomial proportion interval computer.
+  interval_type confidence_interval(
+      const utility::binomial_proportion_interval<real_type>& b) const noexcept {
+    return b(static_cast<real_type>(successes()), static_cast<real_type>(failures()));
   }
 
   bool operator==(const fraction& rhs) const noexcept {
diff --git a/include/boost/histogram/accumulators/weighted_fraction.hpp b/include/boost/histogram/accumulators/weighted_fraction.hpp
index 5fc1e495..d35dc934 100644
--- a/include/boost/histogram/accumulators/weighted_fraction.hpp
+++ b/include/boost/histogram/accumulators/weighted_fraction.hpp
@@ -86,8 +86,7 @@ class weighted_fraction {
   using const_reference = const value_type&;
   using fraction_type = fraction<ValueType>;
   using real_type = typename fraction_type::real_type;
-  using score_type = typename fraction_type::score_type;
-  using interval_type = typename score_type::interval_type;
+  using interval_type = typename fraction_type::interval_type;
 
   weighted_fraction() noexcept = default;
 
@@ -150,10 +149,16 @@ class weighted_fraction {
 
   /// Return standard interval with 68.3 % confidence level (Wilson score interval).
   interval_type confidence_interval() const noexcept {
+    return confidence_interval(utility::wilson_interval<real_type>());
+  }
+
+  /// Return the Wilson score interval.
+  interval_type confidence_interval(
+      const utility::wilson_interval<real_type>& w) const noexcept {
     const real_type n_eff = count();
     const real_type p_hat = value();
-    const real_type correction = score_type::third_order_correction(n_eff);
-    return score_type().wilson_solve_for_neff_phat_correction(n_eff, p_hat, correction);
+    const real_type correction = w.third_order_correction(n_eff);
+    return w.solve_for_neff_phat_correction(n_eff, p_hat, correction);
   }
 
   /// Return the fraction.
diff --git a/include/boost/histogram/utility/wilson_interval.hpp b/include/boost/histogram/utility/wilson_interval.hpp
index 336359ff..ae62d9b3 100644
--- a/include/boost/histogram/utility/wilson_interval.hpp
+++ b/include/boost/histogram/utility/wilson_interval.hpp
@@ -94,7 +94,7 @@ class wilson_interval : public binomial_proportion_interval<ValueType> {
 
     @param p The problem to solve.
   */
-  interval_type wilson_solve_for_neff_phat_correction(
+  interval_type solve_for_neff_phat_correction(
       const value_type& n_eff, const value_type& p_hat,
       const value_type& correction) const noexcept {
     // Equation 41 from this paper: https://arxiv.org/abs/2110.00294
diff --git a/test/accumulators_weighted_fraction_test.cpp b/test/accumulators_weighted_fraction_test.cpp
index 02991f1a..55681ec0 100644
--- a/test/accumulators_weighted_fraction_test.cpp
+++ b/test/accumulators_weighted_fraction_test.cpp
@@ -80,7 +80,7 @@ void run_tests() {
 
     // const auto expected = utility::wilson_interval<double>()(3, 1);
     const auto wilson = utility::wilson_interval<double>();
-    const auto expected = wilson.wilson_solve_for_neff_phat_correction(
+    const auto expected = wilson.solve_for_neff_phat_correction(
         4,      // n_eff = 4
         0.75,   // p_hat = 0.75
         1.46875 // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=4
@@ -99,7 +99,7 @@ void run_tests() {
     const auto ci = f.confidence_interval();
 
     const auto wilson = utility::wilson_interval<double>();
-    const auto expected = wilson.wilson_solve_for_neff_phat_correction(
+    const auto expected = wilson.solve_for_neff_phat_correction(
         1, // n_eff = 1
         0, // p_hat = 0
         10 // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=1
@@ -118,7 +118,7 @@ void run_tests() {
     const auto ci = f.confidence_interval();
 
     const auto wilson = utility::wilson_interval<double>();
-    const auto expected = wilson.wilson_solve_for_neff_phat_correction(
+    const auto expected = wilson.solve_for_neff_phat_correction(
         1, // n_eff = 1
         1, // p_hat = 1
         10 // f(n) = (n³ + n² + 2n + 6) / n³ evaluated at n=1

From f2cc6801774456e5ebaede292bf3dc3cc58afc5c Mon Sep 17 00:00:00 2001
From: Ryan <ryan.elandt@gmail.com>
Date: Tue, 20 Jun 2023 11:43:32 -0400
Subject: [PATCH 7/8] made weighted_fraction friend of fraction

---
 .../boost/histogram/accumulators/fraction.hpp | 34 ++++++++++---------
 1 file changed, 18 insertions(+), 16 deletions(-)

diff --git a/include/boost/histogram/accumulators/fraction.hpp b/include/boost/histogram/accumulators/fraction.hpp
index 2246cdc1..bcb967a3 100644
--- a/include/boost/histogram/accumulators/fraction.hpp
+++ b/include/boost/histogram/accumulators/fraction.hpp
@@ -88,27 +88,11 @@ class fraction {
     return variance_for_p_and_n_eff(value(), count());
   }
 
-  /// Calculate the variance for a given success fraction and effective number of samples.
-  template <class T>
-  static real_type variance_for_p_and_n_eff(const real_type& p, const T& n_eff) noexcept {
-    // We want to compute Var(p) for p = X / n with Var(X) = n p (1 - p)
-    // For Var(X) see
-    // https://en.wikipedia.org/wiki/Binomial_distribution#Expected_value_and_variance
-    // Error propagation: Var(p) = p'(X)^2 Var(X) = p (1 - p) / n
-    return p * (1 - p) / n_eff;
-  }
-
   /// Return standard interval with 68.3 % confidence level (Wilson score interval).
   interval_type confidence_interval() const noexcept {
     return confidence_interval(utility::wilson_interval<real_type>());
   }
 
-  /// Return interval for the given binomial proportion interval computer.
-  interval_type confidence_interval(
-      const utility::binomial_proportion_interval<real_type>& b) const noexcept {
-    return b(static_cast<real_type>(successes()), static_cast<real_type>(failures()));
-  }
-
   bool operator==(const fraction& rhs) const noexcept {
     return succ_ == rhs.succ_ && fail_ == rhs.fail_;
   }
@@ -122,6 +106,24 @@ class fraction {
   }
 
 private:
+  friend class weighted_fraction<value_type>;
+
+  // Calculate the variance for a given success fraction and effective number of samples.
+  template <class T>
+  static real_type variance_for_p_and_n_eff(const real_type& p, const T& n_eff) noexcept {
+    // We want to compute Var(p) for p = X / n with Var(X) = n p (1 - p)
+    // For Var(X) see
+    // https://en.wikipedia.org/wiki/Binomial_distribution#Expected_value_and_variance
+    // Error propagation: Var(p) = p'(X)^2 Var(X) = p (1 - p) / n
+    return p * (1 - p) / n_eff;
+  }
+
+  // Return interval for the given binomial proportion interval computer.
+  interval_type confidence_interval(
+      const utility::binomial_proportion_interval<real_type>& b) const noexcept {
+    return b(static_cast<real_type>(successes()), static_cast<real_type>(failures()));
+  }
+
   value_type succ_{};
   value_type fail_{};
 };

From 4b68641a6355b016c296e72d9166d471ddf945a7 Mon Sep 17 00:00:00 2001
From: Ryan <ryan.elandt@gmail.com>
Date: Tue, 20 Jun 2023 11:44:22 -0400
Subject: [PATCH 8/8] moved weighted fraction to experimental folder

---
 include/boost/histogram/accumulators.hpp                        | 2 +-
 .../{accumulators => experimental}/weighted_fraction.hpp        | 0
 test/accumulators_weighted_fraction_test.cpp                    | 2 +-
 3 files changed, 2 insertions(+), 2 deletions(-)
 rename include/boost/histogram/{accumulators => experimental}/weighted_fraction.hpp (100%)

diff --git a/include/boost/histogram/accumulators.hpp b/include/boost/histogram/accumulators.hpp
index d930ad19..2f0e62b4 100644
--- a/include/boost/histogram/accumulators.hpp
+++ b/include/boost/histogram/accumulators.hpp
@@ -19,10 +19,10 @@
 
 #include <boost/histogram/accumulators/count.hpp>
 #include <boost/histogram/accumulators/fraction.hpp>
-#include <boost/histogram/accumulators/weighted_fraction.hpp>
 #include <boost/histogram/accumulators/mean.hpp>
 #include <boost/histogram/accumulators/sum.hpp>
 #include <boost/histogram/accumulators/weighted_mean.hpp>
 #include <boost/histogram/accumulators/weighted_sum.hpp>
+#include <boost/histogram/experimental/weighted_fraction.hpp>
 
 #endif
diff --git a/include/boost/histogram/accumulators/weighted_fraction.hpp b/include/boost/histogram/experimental/weighted_fraction.hpp
similarity index 100%
rename from include/boost/histogram/accumulators/weighted_fraction.hpp
rename to include/boost/histogram/experimental/weighted_fraction.hpp
diff --git a/test/accumulators_weighted_fraction_test.cpp b/test/accumulators_weighted_fraction_test.cpp
index 55681ec0..28d1606d 100644
--- a/test/accumulators_weighted_fraction_test.cpp
+++ b/test/accumulators_weighted_fraction_test.cpp
@@ -6,7 +6,7 @@
 
 #include <boost/core/lightweight_test.hpp>
 #include <boost/histogram/accumulators/ostream.hpp>
-#include <boost/histogram/accumulators/weighted_fraction.hpp>
+#include <boost/histogram/experimental/weighted_fraction.hpp>
 #include <boost/histogram/utility/wilson_interval.hpp>
 #include <cmath>
 #include <limits>