Re: [PR] Sampling: Exact and Bounded, Probability Proportional to Size [datasketches-cpp]

[via GitHub]

AlexanderSaydakov commented on code in PR #396:
URL: https://github.com/apache/datasketches-cpp/pull/396#discussion_r1364366028


##########
sampling/include/ebpps_sketch.hpp:
##########
@@ -0,0 +1,281 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+#ifndef _EBPPS_SKETCH_HPP_
+#define _EBPPS_SKETCH_HPP_
+
+#include "common_defs.hpp"
+#include "ebpps_sample.hpp"
+#include "optional.hpp"
+#include "serde.hpp"
+
+#include <random>
+#include <string>
+#include <vector>
+
+namespace datasketches {
+
+template<typename A> using AllocU8 = typename 
std::allocator_traits<A>::template rebind_alloc<uint8_t>;
+template<typename A> using vector_u8 = std::vector<uint8_t, AllocU8<A>>;
+
+/// EBPPS sketch constants
+namespace ebpps_constants {
+  /// maximum value of parameter K
+  const uint32_t MAX_K = ((uint32_t) 1 << 31) - 2;
+}
+
+/**
+ * An implementation of an Exact and Bounded Sampling Proportional to Size 
sketch.
+ * 
+ * From: "Exact PPS Sampling with Bounded Sample Size",
+ * B. Hentschel, P. J. Haas, Y. Tian. Information Processing Letters, 2023.
+ * 
+ * This sketch samples data from a stream of items propotional to the weight 
of each item.
+ * The sample guarantees the presence of an item in the result is proportional 
to that item's
+ * portion of the total weight seen by the sketch, and returns a sample no 
larger than size k.
+ * 
+ * The sample may be smaller than k and the resulting size of the sample 
potentially includes
+ * a probabilistic component, meaning the resulting sample size is not always 
constant.
+ *
+ * @author Jon Malkin
+ */
+template<
+  typename T,
+  typename A = std::allocator<T>
+>
+class ebpps_sketch {
+
+  public:
+    static const uint32_t MAX_K = ebpps_constants::MAX_K;
+
+    /**
+     * Constructor
+     * @param k sketch size
+     * @param allocator instance of an allocator
+     */
+    explicit ebpps_sketch(uint32_t k, const A& allocator = A());
+
+    /**
+     * Updates this sketch with the given data item with the given weight.
+     * This method takes an lvalue.
+     * @param item an item from a stream of items
+     * @param weight the weight of the item
+     */
+    void update(const T& item, double weight = 1.0);
+
+    /**
+     * Updates this sketch with the given data item with the given weight.
+     * This method takes an rvalue.
+     * @param item an item from a stream of items
+     * @param weight the weight of the item
+     */
+    void update(T&& item, double weight = 1.0);
+
+    /**
+     * Merges the provided sketch into the current one.
+     * This method takes an lvalue.
+     * @param sketch the sketch to merge into the current object
+     */
+    void merge(const ebpps_sketch<T, A>& sketch);
+
+    /**
+     * Merges the provided sketch into the current one.
+     * This method takes an rvalue.
+     * @param sketch the sketch to merge into the current object
+     */
+    void merge(ebpps_sketch<T, A>&& sketch);
+
+    using result_type = typename ebpps_sample<T,A>::result_type;
+
+    /**
+     * @brief Returns a copy of the current sample, as a std::vector
+     */
+    result_type get_result() const;
+
+    /**
+     * Returns the configured maximum sample size.
+     * @return configured maximum sample size
+     */
+    inline uint32_t get_k() const;
+
+    /**
+     * Returns the number of items processed by the sketch, regardless of
+     * item weight.
+     * @return count of items processed by the sketch
+     */
+    inline uint64_t get_n() const;
+
+    /**
+     * Returns the cumulative weight of items processed by the sketch.
+     * @return cumulative weight of items seen
+     */
+    inline double get_cumulative_weight() const;
+
+    /**
+     * Returns the number of samples currently in the sketch
+     * @return stream length
+     */
+    inline double get_c() const;
+    
+    /**
+     * Returns true if the sketch is empty.
+     * @return empty flag
+     */
+    inline bool is_empty() const;
+    
+    /**
+     * Returns an instance of the allocator for this sketch.
+     * @return allocator
+     */
+    A get_allocator() const;
+
+    /**
+     * Resets the sketch to its default, empty state.
+     */
+    void reset();
+
+    /**
+     * Computes size needed to serialize the current state of the sketch.
+     * @param sd instance of a SerDe
+     * @return size in bytes needed to serialize this sketch
+     */
+    template<typename SerDe = serde<T>>
+    inline size_t get_serialized_size_bytes(const SerDe& sd = SerDe()) const;
+
+    // This is a convenience alias for users
+    // The type returned by the following serialize method
+    using vector_bytes = vector_u8<A>;
+
+    /**
+     * This method serializes the sketch as a vector of bytes.
+     * An optional header can be reserved in front of the sketch.
+     * It is a blank space of a given size.
+     * This header is used in Datasketches PostgreSQL extension.
+     * @param header_size_bytes space to reserve in front of the sketch
+     * @param sd instance of a SerDe
+     */
+    template<typename SerDe = serde<T>>
+    vector_bytes serialize(unsigned header_size_bytes = 0, const SerDe& sd = 
SerDe()) const;
+
+    /**
+     * This method serializes the sketch into a given stream in a binary form
+     * @param os output stream
+     * @param sd instance of a SerDe
+     */
+    template<typename SerDe = serde<T>>
+    void serialize(std::ostream& os, const SerDe& sd = SerDe()) const;
+
+    /**
+     * This method deserializes a sketch from a given array of bytes.
+     * @param bytes pointer to the array of bytes
+     * @param size the size of the array
+     * @param sd instance of a SerDe
+     * @param allocator instance of an allocator
+     * @return an instance of a sketch
+     */
+    template<typename SerDe = serde<T>>
+    static ebpps_sketch deserialize(const void* bytes, size_t size, const 
SerDe& sd = SerDe(), const A& allocator = A());
+
+    /**
+     * This method deserializes a sketch from a given stream.
+     * @param is input stream
+     * @param sd instance of a SerDe
+     * @param allocator instance of an allocator
+     * @return an instance of a sketch
+     */
+    template<typename SerDe = serde<T>>
+    static ebpps_sketch deserialize(std::istream& is, const SerDe& sd = 
SerDe(), const A& allocator = A());
+
+    /**
+     * Prints a summary of the sketch.
+     * @param detail if true, prints values of items
+     * @return the summary as a string
+     */
+    string<A> to_string() const;
+
+    /**
+     * Prints the raw sketch items to a string.
+     * Only works for type T with a defined operator<<() and

Review Comment:
   perhaps worth clarifying: std::ostream& operator<<(std::ostream&, const T&)



##########
sampling/include/ebpps_sketch.hpp:
##########
@@ -0,0 +1,281 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+#ifndef _EBPPS_SKETCH_HPP_
+#define _EBPPS_SKETCH_HPP_
+
+#include "common_defs.hpp"
+#include "ebpps_sample.hpp"
+#include "optional.hpp"
+#include "serde.hpp"
+
+#include <random>
+#include <string>
+#include <vector>
+
+namespace datasketches {
+
+template<typename A> using AllocU8 = typename 
std::allocator_traits<A>::template rebind_alloc<uint8_t>;
+template<typename A> using vector_u8 = std::vector<uint8_t, AllocU8<A>>;
+

Review Comment:
   I would suggest getting rid of these global types.



##########
sampling/include/ebpps_sketch.hpp:
##########
@@ -0,0 +1,281 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+#ifndef _EBPPS_SKETCH_HPP_
+#define _EBPPS_SKETCH_HPP_
+
+#include "common_defs.hpp"
+#include "ebpps_sample.hpp"
+#include "optional.hpp"
+#include "serde.hpp"
+
+#include <random>
+#include <string>
+#include <vector>
+
+namespace datasketches {
+
+template<typename A> using AllocU8 = typename 
std::allocator_traits<A>::template rebind_alloc<uint8_t>;
+template<typename A> using vector_u8 = std::vector<uint8_t, AllocU8<A>>;
+
+/// EBPPS sketch constants
+namespace ebpps_constants {
+  /// maximum value of parameter K
+  const uint32_t MAX_K = ((uint32_t) 1 << 31) - 2;
+}
+
+/**
+ * An implementation of an Exact and Bounded Sampling Proportional to Size 
sketch.
+ * 
+ * From: "Exact PPS Sampling with Bounded Sample Size",
+ * B. Hentschel, P. J. Haas, Y. Tian. Information Processing Letters, 2023.
+ * 
+ * This sketch samples data from a stream of items propotional to the weight 
of each item.
+ * The sample guarantees the presence of an item in the result is proportional 
to that item's
+ * portion of the total weight seen by the sketch, and returns a sample no 
larger than size k.
+ * 
+ * The sample may be smaller than k and the resulting size of the sample 
potentially includes
+ * a probabilistic component, meaning the resulting sample size is not always 
constant.
+ *
+ * @author Jon Malkin
+ */
+template<
+  typename T,
+  typename A = std::allocator<T>
+>
+class ebpps_sketch {
+
+  public:
+    static const uint32_t MAX_K = ebpps_constants::MAX_K;
+
+    /**
+     * Constructor
+     * @param k sketch size
+     * @param allocator instance of an allocator
+     */
+    explicit ebpps_sketch(uint32_t k, const A& allocator = A());
+
+    /**
+     * Updates this sketch with the given data item with the given weight.
+     * This method takes an lvalue.
+     * @param item an item from a stream of items
+     * @param weight the weight of the item
+     */
+    void update(const T& item, double weight = 1.0);
+
+    /**
+     * Updates this sketch with the given data item with the given weight.
+     * This method takes an rvalue.
+     * @param item an item from a stream of items
+     * @param weight the weight of the item
+     */
+    void update(T&& item, double weight = 1.0);
+
+    /**
+     * Merges the provided sketch into the current one.
+     * This method takes an lvalue.
+     * @param sketch the sketch to merge into the current object
+     */
+    void merge(const ebpps_sketch<T, A>& sketch);
+
+    /**
+     * Merges the provided sketch into the current one.
+     * This method takes an rvalue.
+     * @param sketch the sketch to merge into the current object
+     */
+    void merge(ebpps_sketch<T, A>&& sketch);
+
+    using result_type = typename ebpps_sample<T,A>::result_type;
+
+    /**
+     * @brief Returns a copy of the current sample, as a std::vector
+     */
+    result_type get_result() const;
+
+    /**
+     * Returns the configured maximum sample size.
+     * @return configured maximum sample size
+     */
+    inline uint32_t get_k() const;
+
+    /**
+     * Returns the number of items processed by the sketch, regardless of
+     * item weight.
+     * @return count of items processed by the sketch
+     */
+    inline uint64_t get_n() const;
+
+    /**
+     * Returns the cumulative weight of items processed by the sketch.
+     * @return cumulative weight of items seen
+     */
+    inline double get_cumulative_weight() const;
+
+    /**
+     * Returns the number of samples currently in the sketch
+     * @return stream length

Review Comment:
   this must be a mistake



##########
sampling/include/ebpps_sketch_impl.hpp:
##########
@@ -0,0 +1,533 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+#ifndef _EBPPS_SKETCH_IMPL_HPP_
+#define _EBPPS_SKETCH_IMPL_HPP_
+
+#include <memory>
+#include <sstream>
+#include <cmath>
+#include <random>
+#include <algorithm>
+#include <stdexcept>
+#include <utility>
+
+#include "ebpps_sketch.hpp"
+
+namespace datasketches {
+
+template<typename T, typename A>
+ebpps_sketch<T, A>::ebpps_sketch(uint32_t k, const A& allocator) :
+  allocator_(allocator),
+  k_(k),
+  n_(0),
+  cumulative_wt_(0.0),
+  wt_max_(0.0),
+  rho_(1.0),
+  sample_(check_k(k), allocator)
+  {}
+
+template<typename T, typename A>
+ebpps_sketch<T,A>::ebpps_sketch(uint32_t k, uint64_t n, double cumulative_wt,
+                                double wt_max, double rho,
+                                ebpps_sample<T,A>&& sample, const A& 
allocator) :
+  allocator_(allocator),
+  k_(k),
+  n_(n),
+  cumulative_wt_(cumulative_wt),
+  wt_max_(wt_max),
+  rho_(rho),
+  sample_(sample)
+  {}
+
+template<typename T, typename A>
+uint32_t ebpps_sketch<T, A>::get_k() const {
+  return k_;
+}
+
+template<typename T, typename A>
+uint64_t ebpps_sketch<T, A>::get_n() const {
+  return n_;
+}
+
+template<typename T, typename A>
+double ebpps_sketch<T, A>::get_c() const {
+  return sample_.get_c();
+}
+
+template<typename T, typename A>
+double ebpps_sketch<T, A>::get_cumulative_weight() const {
+  return cumulative_wt_;
+}
+
+template<typename T, typename A>
+bool ebpps_sketch<T, A>::is_empty() const {
+  return n_ == 0;
+}
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::reset() {
+  n_ = 0;
+  cumulative_wt_ = 0.0;
+  wt_max_ = 0.0;
+  rho_ = 1.0;
+  sample_.reset();
+}
+
+template<typename T, typename A>
+string<A> ebpps_sketch<T, A>::to_string() const {
+  // Using a temporary stream for implementation here does not comply with 
AllocatorAwareContainer requirements.
+  // The stream does not support passing an allocator instance, and 
alternatives are complicated.
+  std::ostringstream os;
+  os << "### EBPPS Sketch SUMMARY:" << std::endl;
+  os << "   k            : " << k_ << std::endl;
+  os << "   n            : " << n_ << std::endl;
+  os << "   cum. weight  : " << cumulative_wt_ << std::endl;
+  os << "   wt_mac       : " << wt_max_ << std::endl;
+  os << "   rho          : " << rho_ << std::endl;
+  os << "   C            : " << sample_.get_c() << std::endl;
+  os << "### END SKETCH SUMMARY" << std::endl;
+  return string<A>(os.str().c_str(), allocator_);
+}
+
+template<typename T, typename A>
+string<A> ebpps_sketch<T, A>::items_to_string() const {
+  // Using a temporary stream for implementation here does not comply with 
AllocatorAwareContainer requirements.
+  // The stream does not support passing an allocator instance, and 
alternatives are complicated.
+  std::ostringstream os;
+  os << "### Sketch Items" << std::endl;
+  os << sample_.to_string(); // assumes std::endl at end
+  return string<A>(os.str().c_str(), allocator_);
+}
+
+template<typename T, typename A>
+A ebpps_sketch<T, A>::get_allocator() const {
+  return allocator_;
+}
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::update(const T& item, double weight) {
+  return internal_update(item, weight);
+}
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::update(T&& item, double weight) {
+  return internal_update(std::move(item), weight);
+}
+
+template<typename T, typename A>
+template<typename FwdItem>
+void ebpps_sketch<T, A>::internal_update(FwdItem&& item, double weight) {
+  if (weight < 0.0 || std::isnan(weight) || std::isinf(weight)) {
+    throw std::invalid_argument("Item weights must be nonnegative and finite. 
Found: "
+                                + std::to_string(weight));
+  } else if (weight == 0.0) {
+    return;
+  }
+
+  double new_cum_wt = cumulative_wt_ + weight;
+  double new_wt_max = std::max(wt_max_, weight);
+  double new_rho = std::min(1.0 / new_wt_max, k_ / new_cum_wt);
+
+  if (cumulative_wt_ > 0.0)
+    sample_.downsample(new_rho / rho_);
+  
+  ebpps_sample<T,A> tmp(conditional_forward<FwdItem>(item), new_rho * weight, 
allocator_);
+
+  sample_.merge(tmp);
+
+  cumulative_wt_ = new_cum_wt;
+  wt_max_ = new_wt_max;
+  rho_ = new_rho;
+  ++n_;
+}
+
+template<typename T, typename A>
+auto ebpps_sketch<T,A>::get_result() const -> result_type {
+  return sample_.get_sample();
+}
+
+/* Merging
+ * There is a trivial merge algorithm that involves downsampling each sketch A 
and B
+ * as A.cum_wt / (A.cum_wt + B.cum_wt) and B.cum_wt / (A.cum_wt + B.cum_wt),
+ * respectively. That merge does preserve first-order probabilities, 
specifically
+ * the probability proportional to size property, and like all other known 
merge
+ * algorithms distorts second-order probabilities (co-occurrences). There are
+ * pathological cases, most obvious with k=2 and A.cum_wt == B.cum_wt where 
that
+ * approach will always take exactly 1 item from A and 1 from B, meaning the
+ * co-occurrence rate for two items from either sketch is guaranteed to be 0.0.
+ * 
+ * With EBPPS, once an item is accepted into the sketch we no longer need to
+ * track the item's weight: All accepted items are treated equally. As a 
result, we
+ * can take inspiration from the reservoir sampling merge in the 
datasketches-java
+ * library. We need to merge the smaller sketch into the larger one, swapping 
as
+ * needed to ensure that, at which point we simply call update() with the items
+ * in the smaller sketch as long as we adjust the weight appropriately.
+ * Merging smaller into larger is essential to ensure that no item has a
+ * contribution to C > 1.0.
+ */
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::merge(ebpps_sketch<T, A>&& sk) {
+  if (sk.get_cumulative_weight() == 0.0) return;
+  else if (sk.get_cumulative_weight() > get_cumulative_weight()) {
+    // need to swap this with sk to merge smaller into larger
+    std::swap(*this, sk);
+  }
+
+  internal_merge(sk);
+}
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::merge(const ebpps_sketch<T, A>& sk) {
+  if (sk.get_cumulative_weight() == 0.0) return;
+  else if (sk.get_cumulative_weight() > get_cumulative_weight()) {
+    // need to swap this with sk to merge, so make a copy, swap,
+    // and use that to merge
+    ebpps_sketch sk_copy(sk);
+    swap(*this, sk_copy);
+    internal_merge(sk_copy);
+  } else {
+    internal_merge(sk);
+  }
+}
+
+template<typename T, typename A>
+template<typename O>
+void ebpps_sketch<T, A>::internal_merge(O&& sk) {
+  // assumes that sk.cumulative_wt_ <= cumulative_wt_,
+  // which must be checked before calling this
+  
+  const ebpps_sample<T,A>& other_sample = sk.sample_;
+
+  double final_cum_wt = cumulative_wt_ + sk.cumulative_wt_;

Review Comment:
   const?



##########
sampling/include/ebpps_sample_impl.hpp:
##########
@@ -0,0 +1,535 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+#ifndef _EBPPS_SAMPLE_IMPL_HPP_
+#define _EBPPS_SAMPLE_IMPL_HPP_
+
+#include "common_defs.hpp"
+#include "conditional_forward.hpp"
+#include "ebpps_sample.hpp"
+#include "serde.hpp"
+
+#include <cmath>
+#include <string>
+#include <sstream>
+
+namespace datasketches {
+
+template<typename T, typename A>
+ebpps_sample<T,A>::ebpps_sample(uint32_t reserved_size, const A& allocator) :
+  allocator_(allocator),
+  c_(0.0),
+  partial_item_(),
+  data_(allocator)
+  {
+    data_.reserve(reserved_size);
+  }
+
+template<typename T, typename A>
+template<typename TT>
+ebpps_sample<T,A>::ebpps_sample(TT&& item, double theta, const A& allocator) :
+  allocator_(allocator),
+  c_(theta),
+  partial_item_(),
+  data_(allocator)
+  {
+    if (theta == 1.0) {
+      data_.reserve(1);
+      data_.emplace_back(std::forward<TT>(item));
+    } else {
+      partial_item_.emplace(std::forward<TT>(item));
+    }
+  }
+
+template<typename T, typename A>
+ebpps_sample<T,A>::ebpps_sample(std::vector<T, A>&& data, optional<T>&& 
partial_item, double c, const A& allocator) :
+  allocator_(allocator),
+  c_(c),
+  partial_item_(partial_item),
+  data_(data, allocator)
+  {}
+
+template<typename T, typename A>
+auto ebpps_sample<T,A>::get_sample() const -> result_type {
+  double unused;
+  double c_frac = std::modf(c_, &unused);
+  bool include_partial = next_double() < c_frac;
+  uint32_t result_size = static_cast<uint32_t>(data_.size()) + 
(include_partial ? 1 : 0);
+
+  result_type result;
+  result.reserve(result_size);
+  std::copy(data_.begin(), data_.end(), std::back_inserter(result));
+  if (include_partial)
+    result.emplace_back(static_cast<const T&>(*partial_item_));
+  
+  return result;
+}
+
+template<typename T, typename A>
+void ebpps_sample<T,A>::downsample(double theta) {
+  if (theta >= 1.0) return;
+
+  double new_c = theta * c_;
+  double new_c_int;
+  double new_c_frac = std::modf(new_c, &new_c_int);
+  double c_int;
+  double c_frac = std::modf(c_, &c_int);
+
+  if (new_c_int == 0.0) {
+    // no full items retained
+    if (next_double() > (c_frac / c_)) {
+      swap_with_partial();
+    }
+    data_.clear();
+  } else if (new_c_int == c_int) {
+    // no items deleted
+    if (next_double() > (1 - theta * c_frac)/(1 - new_c_frac)) {
+      swap_with_partial();
+    }
+  } else {
+    if (next_double() < theta * c_frac) {
+      // subsample data in random order; last item is partial
+      // create sample size new_c_int then swap_with_partial()
+      subsample(static_cast<uint32_t>(new_c_int));
+      swap_with_partial();
+    } else {
+      // create sample size new_c_int + 1 then move_one_to_partial)
+      subsample(static_cast<uint32_t>(new_c_int) + 1);
+      move_one_to_partial();
+    }
+  }
+
+  if (new_c == new_c_int)
+    partial_item_.reset();
+
+  c_ = new_c;
+}
+
+template<typename T, typename A>
+template<typename FwdSample>
+void ebpps_sample<T,A>::merge(FwdSample&& other) {
+  double c_int;
+  double c_frac = std::modf(c_, &c_int);
+
+  double unused;
+  double other_c_frac = std::modf(other.c_, &unused);
+
+  // update c_ here but do NOT recompute fractional part yet
+  c_ += other.c_;
+
+  for (uint32_t i = 0; i < other.data_.size(); ++i)
+    data_.emplace_back(conditional_forward<FwdSample>(other.data_[i]));
+
+  // This modifies the original algorithm slightly due to numeric
+  // precision issues. Specifically, the test if c_frac + other_c_frac == 1.0
+  // happens before tests for < 1.0 or > 1.0 and can also be triggered
+  // if c_ == floor(c_) (the updated value of c_, not the input).
+  //
+  // We can still run into issues where c_frac + other_c_frac == epsilon
+  // and the first case would have ideally triggered. As a result, we must
+  // check if the partial item exists before adding to the data_ vector.
+
+  if (c_frac == 0.0 && other_c_frac == 0.0) {
+    partial_item_.reset();
+  } else if (c_frac + other_c_frac == 1.0 || c_ == std::floor(c_)) {
+    if (next_double() <= c_frac) {
+      if (partial_item_)
+        data_.emplace_back(std::move(*partial_item_));
+    } else {
+      if (other.partial_item_)
+        
data_.emplace_back(conditional_forward<FwdSample>(*other.partial_item_));
+    }
+    partial_item_.reset();
+  } else if (c_frac + other_c_frac < 1.0) {
+    if (next_double() > c_frac / (c_frac + other_c_frac)) {
+      set_partial(conditional_forward<FwdSample>(*other.partial_item_));
+    }
+  } else { // c_frac + other_c_frac > 1
+    if (next_double() <= (1 - c_frac) / ((1 - c_frac) + (1 - other_c_frac))) {
+      data_.emplace_back(conditional_forward<FwdSample>(*other.partial_item_));
+    } else {
+      data_.emplace_back(std::move(*partial_item_));
+      partial_item_.reset();
+      set_partial(conditional_forward<FwdSample>(*other.partial_item_));
+    }
+  }
+}
+
+template<typename T, typename A>
+string<A> ebpps_sample<T ,A>::to_string() const {
+  std::ostringstream oss;
+  oss << "   sample:" << std::endl;
+  uint32_t idx = 0;
+  for (const T& item : data_)
+    oss << "\t" << idx++ << ":\t" << item << std::endl;
+  oss << "   partial: " << (bool(partial_item_) ? (*partial_item_) : "NULL") 
<< std::endl;
+
+  return oss.str();
+}
+
+template<typename T, typename A>
+void ebpps_sample<T,A>::subsample(uint32_t num_samples) {
+  // we can perform a Fisher-Yates style shuffle, stopping after
+  // num_samples points since subsequent swaps would only be
+  // between items after num_samples. This is valid since a
+  // point from anywhere in the initial array would be eligible
+  // to end up in the final subsample.
+
+  if (num_samples == data_.size()) return;
+
+  auto erase_start = data_.begin();
+  uint32_t data_len = static_cast<uint32_t>(data_.size());
+  for (uint32_t i = 0; i < num_samples; ++i, ++erase_start) {
+    uint32_t j = i + random_idx(data_len - i);
+    std::swap(data_[i], data_[j]);
+  }
+
+  data_.erase(erase_start, data_.end());
+}
+
+template<typename T, typename A>
+template<typename FwdItem>
+void ebpps_sample<T,A>::set_partial(FwdItem&& item) {
+  if (partial_item_)
+    *partial_item_ = conditional_forward<FwdItem>(item);
+  else
+    partial_item_.emplace(conditional_forward<FwdItem>(item));
+}
+
+template<typename T, typename A>
+void ebpps_sample<T,A>::move_one_to_partial() {
+  size_t idx = random_idx(static_cast<uint32_t>(data_.size()));
+  // swap selected item to end so we can delete it easily
+  size_t last_idx = data_.size() - 1;
+  if (idx != last_idx) {
+    std::swap(data_[idx], data_[last_idx]);
+  }
+
+  set_partial(std::move(data_[last_idx]));
+
+  data_.pop_back();
+}
+
+template<typename T, typename A>
+void ebpps_sample<T,A>::swap_with_partial() {
+  if (partial_item_) {
+    size_t idx = random_idx(static_cast<uint32_t>(data_.size()));
+    std::swap(data_[idx], *partial_item_);
+  } else {
+    move_one_to_partial();
+  }
+}
+
+template<typename T, typename A>
+void ebpps_sample<T,A>::reset() {
+  c_ = 0.0;
+  partial_item_.reset();
+  data_.clear();
+}
+
+template<typename T, typename A>
+double ebpps_sample<T,A>::get_c() const {
+  return c_;
+}
+
+template<typename T, typename A>
+auto ebpps_sample<T,A>::get_full_items() const -> result_type {
+  return result_type(data_);
+}
+
+template<typename T, typename A>
+bool ebpps_sample<T,A>::has_partial_item() const {
+  return bool(partial_item_);
+}
+
+template<typename T, typename A>
+T ebpps_sample<T,A>::get_partial_item() const {
+  if (!partial_item_) throw std::runtime_error("Call to get_partial_item() 
when no partial item exists");
+  return *partial_item_;
+}
+
+template<typename T, typename A>
+uint32_t ebpps_sample<T,A>::random_idx(uint32_t max) {
+  std::uniform_int_distribution<uint32_t> dist(0, max - 1);

Review Comment:
   should this be static to initialize once?



##########
sampling/include/ebpps_sketch.hpp:
##########
@@ -0,0 +1,281 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+#ifndef _EBPPS_SKETCH_HPP_
+#define _EBPPS_SKETCH_HPP_
+
+#include "common_defs.hpp"
+#include "ebpps_sample.hpp"
+#include "optional.hpp"
+#include "serde.hpp"
+
+#include <random>
+#include <string>
+#include <vector>
+
+namespace datasketches {
+
+template<typename A> using AllocU8 = typename 
std::allocator_traits<A>::template rebind_alloc<uint8_t>;
+template<typename A> using vector_u8 = std::vector<uint8_t, AllocU8<A>>;
+
+/// EBPPS sketch constants
+namespace ebpps_constants {
+  /// maximum value of parameter K
+  const uint32_t MAX_K = ((uint32_t) 1 << 31) - 2;
+}
+
+/**
+ * An implementation of an Exact and Bounded Sampling Proportional to Size 
sketch.
+ * 
+ * From: "Exact PPS Sampling with Bounded Sample Size",
+ * B. Hentschel, P. J. Haas, Y. Tian. Information Processing Letters, 2023.
+ * 
+ * This sketch samples data from a stream of items propotional to the weight 
of each item.
+ * The sample guarantees the presence of an item in the result is proportional 
to that item's
+ * portion of the total weight seen by the sketch, and returns a sample no 
larger than size k.
+ * 
+ * The sample may be smaller than k and the resulting size of the sample 
potentially includes
+ * a probabilistic component, meaning the resulting sample size is not always 
constant.
+ *
+ * @author Jon Malkin
+ */
+template<
+  typename T,
+  typename A = std::allocator<T>
+>
+class ebpps_sketch {
+
+  public:
+    static const uint32_t MAX_K = ebpps_constants::MAX_K;
+
+    /**
+     * Constructor
+     * @param k sketch size
+     * @param allocator instance of an allocator
+     */
+    explicit ebpps_sketch(uint32_t k, const A& allocator = A());
+
+    /**
+     * Updates this sketch with the given data item with the given weight.
+     * This method takes an lvalue.
+     * @param item an item from a stream of items
+     * @param weight the weight of the item
+     */
+    void update(const T& item, double weight = 1.0);
+
+    /**
+     * Updates this sketch with the given data item with the given weight.
+     * This method takes an rvalue.
+     * @param item an item from a stream of items
+     * @param weight the weight of the item
+     */
+    void update(T&& item, double weight = 1.0);
+
+    /**
+     * Merges the provided sketch into the current one.
+     * This method takes an lvalue.
+     * @param sketch the sketch to merge into the current object
+     */
+    void merge(const ebpps_sketch<T, A>& sketch);
+
+    /**
+     * Merges the provided sketch into the current one.
+     * This method takes an rvalue.
+     * @param sketch the sketch to merge into the current object
+     */
+    void merge(ebpps_sketch<T, A>&& sketch);
+
+    using result_type = typename ebpps_sample<T,A>::result_type;
+
+    /**
+     * @brief Returns a copy of the current sample, as a std::vector
+     */
+    result_type get_result() const;
+
+    /**
+     * Returns the configured maximum sample size.
+     * @return configured maximum sample size
+     */
+    inline uint32_t get_k() const;
+
+    /**
+     * Returns the number of items processed by the sketch, regardless of
+     * item weight.
+     * @return count of items processed by the sketch
+     */
+    inline uint64_t get_n() const;
+
+    /**
+     * Returns the cumulative weight of items processed by the sketch.
+     * @return cumulative weight of items seen
+     */
+    inline double get_cumulative_weight() const;
+
+    /**
+     * Returns the number of samples currently in the sketch
+     * @return stream length
+     */
+    inline double get_c() const;
+    
+    /**
+     * Returns true if the sketch is empty.
+     * @return empty flag
+     */
+    inline bool is_empty() const;
+    
+    /**
+     * Returns an instance of the allocator for this sketch.
+     * @return allocator
+     */
+    A get_allocator() const;
+
+    /**
+     * Resets the sketch to its default, empty state.
+     */
+    void reset();
+
+    /**
+     * Computes size needed to serialize the current state of the sketch.
+     * @param sd instance of a SerDe
+     * @return size in bytes needed to serialize this sketch
+     */
+    template<typename SerDe = serde<T>>
+    inline size_t get_serialized_size_bytes(const SerDe& sd = SerDe()) const;
+
+    // This is a convenience alias for users
+    // The type returned by the following serialize method
+    using vector_bytes = vector_u8<A>;

Review Comment:
   this can be
   using vector_bytes = std::vector<uint8_t, typename 
std::allocator_traits<Allocator>::template rebind_alloc<uint8_t>>;
   so that the vector_u8 would be unnecessary



##########
sampling/include/ebpps_sketch.hpp:
##########
@@ -0,0 +1,281 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+#ifndef _EBPPS_SKETCH_HPP_
+#define _EBPPS_SKETCH_HPP_
+
+#include "common_defs.hpp"
+#include "ebpps_sample.hpp"
+#include "optional.hpp"
+#include "serde.hpp"
+
+#include <random>
+#include <string>
+#include <vector>
+
+namespace datasketches {
+
+template<typename A> using AllocU8 = typename 
std::allocator_traits<A>::template rebind_alloc<uint8_t>;
+template<typename A> using vector_u8 = std::vector<uint8_t, AllocU8<A>>;
+
+/// EBPPS sketch constants
+namespace ebpps_constants {
+  /// maximum value of parameter K
+  const uint32_t MAX_K = ((uint32_t) 1 << 31) - 2;
+}
+
+/**
+ * An implementation of an Exact and Bounded Sampling Proportional to Size 
sketch.
+ * 
+ * From: "Exact PPS Sampling with Bounded Sample Size",
+ * B. Hentschel, P. J. Haas, Y. Tian. Information Processing Letters, 2023.
+ * 
+ * This sketch samples data from a stream of items propotional to the weight 
of each item.
+ * The sample guarantees the presence of an item in the result is proportional 
to that item's
+ * portion of the total weight seen by the sketch, and returns a sample no 
larger than size k.
+ * 
+ * The sample may be smaller than k and the resulting size of the sample 
potentially includes
+ * a probabilistic component, meaning the resulting sample size is not always 
constant.
+ *
+ * @author Jon Malkin
+ */
+template<
+  typename T,
+  typename A = std::allocator<T>
+>
+class ebpps_sketch {
+
+  public:
+    static const uint32_t MAX_K = ebpps_constants::MAX_K;
+
+    /**
+     * Constructor
+     * @param k sketch size
+     * @param allocator instance of an allocator
+     */
+    explicit ebpps_sketch(uint32_t k, const A& allocator = A());
+
+    /**
+     * Updates this sketch with the given data item with the given weight.
+     * This method takes an lvalue.
+     * @param item an item from a stream of items
+     * @param weight the weight of the item
+     */
+    void update(const T& item, double weight = 1.0);
+
+    /**
+     * Updates this sketch with the given data item with the given weight.
+     * This method takes an rvalue.
+     * @param item an item from a stream of items
+     * @param weight the weight of the item
+     */
+    void update(T&& item, double weight = 1.0);
+
+    /**
+     * Merges the provided sketch into the current one.
+     * This method takes an lvalue.
+     * @param sketch the sketch to merge into the current object
+     */
+    void merge(const ebpps_sketch<T, A>& sketch);
+
+    /**
+     * Merges the provided sketch into the current one.
+     * This method takes an rvalue.
+     * @param sketch the sketch to merge into the current object
+     */
+    void merge(ebpps_sketch<T, A>&& sketch);
+
+    using result_type = typename ebpps_sample<T,A>::result_type;
+
+    /**
+     * @brief Returns a copy of the current sample, as a std::vector
+     */
+    result_type get_result() const;
+
+    /**
+     * Returns the configured maximum sample size.
+     * @return configured maximum sample size
+     */
+    inline uint32_t get_k() const;
+
+    /**
+     * Returns the number of items processed by the sketch, regardless of
+     * item weight.
+     * @return count of items processed by the sketch
+     */
+    inline uint64_t get_n() const;
+
+    /**
+     * Returns the cumulative weight of items processed by the sketch.
+     * @return cumulative weight of items seen
+     */
+    inline double get_cumulative_weight() const;
+
+    /**
+     * Returns the number of samples currently in the sketch
+     * @return stream length
+     */
+    inline double get_c() const;
+    
+    /**
+     * Returns true if the sketch is empty.
+     * @return empty flag
+     */
+    inline bool is_empty() const;
+    
+    /**
+     * Returns an instance of the allocator for this sketch.
+     * @return allocator
+     */
+    A get_allocator() const;
+
+    /**
+     * Resets the sketch to its default, empty state.
+     */
+    void reset();
+
+    /**
+     * Computes size needed to serialize the current state of the sketch.
+     * @param sd instance of a SerDe
+     * @return size in bytes needed to serialize this sketch
+     */
+    template<typename SerDe = serde<T>>
+    inline size_t get_serialized_size_bytes(const SerDe& sd = SerDe()) const;
+
+    // This is a convenience alias for users
+    // The type returned by the following serialize method
+    using vector_bytes = vector_u8<A>;
+
+    /**
+     * This method serializes the sketch as a vector of bytes.
+     * An optional header can be reserved in front of the sketch.
+     * It is a blank space of a given size.
+     * This header is used in Datasketches PostgreSQL extension.
+     * @param header_size_bytes space to reserve in front of the sketch
+     * @param sd instance of a SerDe
+     */
+    template<typename SerDe = serde<T>>
+    vector_bytes serialize(unsigned header_size_bytes = 0, const SerDe& sd = 
SerDe()) const;
+
+    /**
+     * This method serializes the sketch into a given stream in a binary form
+     * @param os output stream
+     * @param sd instance of a SerDe
+     */
+    template<typename SerDe = serde<T>>
+    void serialize(std::ostream& os, const SerDe& sd = SerDe()) const;
+
+    /**
+     * This method deserializes a sketch from a given array of bytes.
+     * @param bytes pointer to the array of bytes
+     * @param size the size of the array
+     * @param sd instance of a SerDe
+     * @param allocator instance of an allocator
+     * @return an instance of a sketch
+     */
+    template<typename SerDe = serde<T>>
+    static ebpps_sketch deserialize(const void* bytes, size_t size, const 
SerDe& sd = SerDe(), const A& allocator = A());
+
+    /**
+     * This method deserializes a sketch from a given stream.
+     * @param is input stream
+     * @param sd instance of a SerDe
+     * @param allocator instance of an allocator
+     * @return an instance of a sketch
+     */
+    template<typename SerDe = serde<T>>
+    static ebpps_sketch deserialize(std::istream& is, const SerDe& sd = 
SerDe(), const A& allocator = A());
+
+    /**
+     * Prints a summary of the sketch.
+     * @param detail if true, prints values of items
+     * @return the summary as a string
+     */
+    string<A> to_string() const;
+
+    /**
+     * Prints the raw sketch items to a string.
+     * Only works for type T with a defined operator<<() and
+     * kept separate from to_string() to allow compilation even if
+     * T does not have such an operator defined.
+     * @return a string with the sketch items
+     */
+    string<A> items_to_string() const;
+
+    /**
+     * Iterator pointing to the first item in the sketch.
+     * If the sketch is empty, the returned iterator must not be dereferenced 
or incremented.
+     * @return iterator pointing to the first item in the sketch
+     */
+    typename ebpps_sample<T,A>::const_iterator begin() const;
+
+    /**
+     * Iterator pointing to the past-the-end item in the sketch.
+     * The past-the-end item is the hypothetical item that would follow the 
last item.
+     * It does not point to any item, and must not be dereferenced or 
incremented.
+     * @return iterator pointing to the past-the-end item in the sketch
+     */
+    typename ebpps_sample<T,A>::const_iterator end() const;
+
+  private:
+    typedef typename std::allocator_traits<A>::template rebind_alloc<double> 
AllocDouble;

Review Comment:
   I thought we prefer "using"



##########
sampling/include/ebpps_sketch.hpp:
##########
@@ -0,0 +1,281 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+#ifndef _EBPPS_SKETCH_HPP_
+#define _EBPPS_SKETCH_HPP_
+
+#include "common_defs.hpp"
+#include "ebpps_sample.hpp"
+#include "optional.hpp"
+#include "serde.hpp"
+
+#include <random>
+#include <string>
+#include <vector>
+
+namespace datasketches {
+
+template<typename A> using AllocU8 = typename 
std::allocator_traits<A>::template rebind_alloc<uint8_t>;
+template<typename A> using vector_u8 = std::vector<uint8_t, AllocU8<A>>;
+
+/// EBPPS sketch constants
+namespace ebpps_constants {
+  /// maximum value of parameter K
+  const uint32_t MAX_K = ((uint32_t) 1 << 31) - 2;
+}
+
+/**
+ * An implementation of an Exact and Bounded Sampling Proportional to Size 
sketch.
+ * 
+ * From: "Exact PPS Sampling with Bounded Sample Size",
+ * B. Hentschel, P. J. Haas, Y. Tian. Information Processing Letters, 2023.
+ * 
+ * This sketch samples data from a stream of items propotional to the weight 
of each item.
+ * The sample guarantees the presence of an item in the result is proportional 
to that item's
+ * portion of the total weight seen by the sketch, and returns a sample no 
larger than size k.
+ * 
+ * The sample may be smaller than k and the resulting size of the sample 
potentially includes
+ * a probabilistic component, meaning the resulting sample size is not always 
constant.
+ *
+ * @author Jon Malkin
+ */
+template<
+  typename T,
+  typename A = std::allocator<T>
+>
+class ebpps_sketch {
+
+  public:
+    static const uint32_t MAX_K = ebpps_constants::MAX_K;
+
+    /**
+     * Constructor
+     * @param k sketch size
+     * @param allocator instance of an allocator
+     */
+    explicit ebpps_sketch(uint32_t k, const A& allocator = A());
+
+    /**
+     * Updates this sketch with the given data item with the given weight.
+     * This method takes an lvalue.
+     * @param item an item from a stream of items
+     * @param weight the weight of the item
+     */
+    void update(const T& item, double weight = 1.0);
+
+    /**
+     * Updates this sketch with the given data item with the given weight.
+     * This method takes an rvalue.
+     * @param item an item from a stream of items
+     * @param weight the weight of the item
+     */
+    void update(T&& item, double weight = 1.0);
+
+    /**
+     * Merges the provided sketch into the current one.
+     * This method takes an lvalue.
+     * @param sketch the sketch to merge into the current object
+     */
+    void merge(const ebpps_sketch<T, A>& sketch);
+
+    /**
+     * Merges the provided sketch into the current one.
+     * This method takes an rvalue.
+     * @param sketch the sketch to merge into the current object
+     */
+    void merge(ebpps_sketch<T, A>&& sketch);
+
+    using result_type = typename ebpps_sample<T,A>::result_type;
+
+    /**
+     * @brief Returns a copy of the current sample, as a std::vector
+     */
+    result_type get_result() const;
+
+    /**
+     * Returns the configured maximum sample size.
+     * @return configured maximum sample size
+     */
+    inline uint32_t get_k() const;
+
+    /**
+     * Returns the number of items processed by the sketch, regardless of
+     * item weight.
+     * @return count of items processed by the sketch
+     */
+    inline uint64_t get_n() const;
+
+    /**
+     * Returns the cumulative weight of items processed by the sketch.
+     * @return cumulative weight of items seen
+     */
+    inline double get_cumulative_weight() const;
+
+    /**
+     * Returns the number of samples currently in the sketch
+     * @return stream length
+     */
+    inline double get_c() const;
+    
+    /**
+     * Returns true if the sketch is empty.
+     * @return empty flag
+     */
+    inline bool is_empty() const;
+    
+    /**
+     * Returns an instance of the allocator for this sketch.
+     * @return allocator
+     */
+    A get_allocator() const;
+
+    /**
+     * Resets the sketch to its default, empty state.
+     */
+    void reset();
+
+    /**
+     * Computes size needed to serialize the current state of the sketch.
+     * @param sd instance of a SerDe
+     * @return size in bytes needed to serialize this sketch
+     */
+    template<typename SerDe = serde<T>>
+    inline size_t get_serialized_size_bytes(const SerDe& sd = SerDe()) const;
+
+    // This is a convenience alias for users
+    // The type returned by the following serialize method
+    using vector_bytes = vector_u8<A>;
+
+    /**
+     * This method serializes the sketch as a vector of bytes.
+     * An optional header can be reserved in front of the sketch.
+     * It is a blank space of a given size.
+     * This header is used in Datasketches PostgreSQL extension.
+     * @param header_size_bytes space to reserve in front of the sketch
+     * @param sd instance of a SerDe
+     */
+    template<typename SerDe = serde<T>>
+    vector_bytes serialize(unsigned header_size_bytes = 0, const SerDe& sd = 
SerDe()) const;
+
+    /**
+     * This method serializes the sketch into a given stream in a binary form
+     * @param os output stream
+     * @param sd instance of a SerDe
+     */
+    template<typename SerDe = serde<T>>
+    void serialize(std::ostream& os, const SerDe& sd = SerDe()) const;
+
+    /**
+     * This method deserializes a sketch from a given array of bytes.
+     * @param bytes pointer to the array of bytes
+     * @param size the size of the array
+     * @param sd instance of a SerDe
+     * @param allocator instance of an allocator
+     * @return an instance of a sketch
+     */
+    template<typename SerDe = serde<T>>
+    static ebpps_sketch deserialize(const void* bytes, size_t size, const 
SerDe& sd = SerDe(), const A& allocator = A());
+
+    /**
+     * This method deserializes a sketch from a given stream.
+     * @param is input stream
+     * @param sd instance of a SerDe
+     * @param allocator instance of an allocator
+     * @return an instance of a sketch
+     */
+    template<typename SerDe = serde<T>>
+    static ebpps_sketch deserialize(std::istream& is, const SerDe& sd = 
SerDe(), const A& allocator = A());
+
+    /**
+     * Prints a summary of the sketch.
+     * @param detail if true, prints values of items
+     * @return the summary as a string
+     */
+    string<A> to_string() const;
+
+    /**
+     * Prints the raw sketch items to a string.
+     * Only works for type T with a defined operator<<() and
+     * kept separate from to_string() to allow compilation even if
+     * T does not have such an operator defined.
+     * @return a string with the sketch items
+     */
+    string<A> items_to_string() const;
+
+    /**
+     * Iterator pointing to the first item in the sketch.
+     * If the sketch is empty, the returned iterator must not be dereferenced 
or incremented.
+     * @return iterator pointing to the first item in the sketch
+     */
+    typename ebpps_sample<T,A>::const_iterator begin() const;
+
+    /**
+     * Iterator pointing to the past-the-end item in the sketch.
+     * The past-the-end item is the hypothetical item that would follow the 
last item.
+     * It does not point to any item, and must not be dereferenced or 
incremented.
+     * @return iterator pointing to the past-the-end item in the sketch
+     */
+    typename ebpps_sample<T,A>::const_iterator end() const;
+
+  private:
+    typedef typename std::allocator_traits<A>::template rebind_alloc<double> 
AllocDouble;
+
+    static const uint8_t PREAMBLE_LONGS_EMPTY  = 1;
+    static const uint8_t PREAMBLE_LONGS_FULL   = 5; // C is part of sample_
+    static const uint8_t SER_VER = 1;
+    static const uint8_t FAMILY_ID  = 19;
+    static const uint8_t EMPTY_FLAG_MASK  = 4;
+    static const uint8_t HAS_PARTIAL_ITEM_MASK = 8;
+
+    A allocator_;
+    uint32_t k_;                    // max size of sketch, in items
+    uint64_t n_;                    // total number of items processed by the 
sketch
+
+    double cumulative_wt_;          // total weight of items processed by the 
sketch
+    double wt_max_;                 // maximum weight seen so far
+    double rho_;                    // latest scaling parameter for 
downsampling
+
+    ebpps_sample<T,A> sample_;      // Object holding the current state of the 
sample
+
+    // handles merge after ensuring other.cumulative_wt_ <= 
this->cumulative_wt_
+    // so we can send items in individually
+    template<typename O>
+    void internal_merge(O&& other);
+
+    ebpps_sketch(uint32_t k, uint64_t n, double cumulative_wt, double wt_max, 
double rho,
+                 ebpps_sample<T,A>&& sample, const A& allocator = A());
+
+    string<A> items_to_string(bool print_gap) const;
+
+    // internal-use-only update

Review Comment:
   unnecessary comment. the method has "internal" in its name and it is private.



##########
sampling/include/ebpps_sample_impl.hpp:
##########
@@ -0,0 +1,535 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+#ifndef _EBPPS_SAMPLE_IMPL_HPP_
+#define _EBPPS_SAMPLE_IMPL_HPP_
+
+#include "common_defs.hpp"
+#include "conditional_forward.hpp"
+#include "ebpps_sample.hpp"
+#include "serde.hpp"
+
+#include <cmath>
+#include <string>
+#include <sstream>
+
+namespace datasketches {
+
+template<typename T, typename A>
+ebpps_sample<T,A>::ebpps_sample(uint32_t reserved_size, const A& allocator) :
+  allocator_(allocator),
+  c_(0.0),
+  partial_item_(),
+  data_(allocator)
+  {
+    data_.reserve(reserved_size);
+  }
+
+template<typename T, typename A>
+template<typename TT>
+ebpps_sample<T,A>::ebpps_sample(TT&& item, double theta, const A& allocator) :
+  allocator_(allocator),
+  c_(theta),
+  partial_item_(),
+  data_(allocator)
+  {
+    if (theta == 1.0) {
+      data_.reserve(1);
+      data_.emplace_back(std::forward<TT>(item));
+    } else {
+      partial_item_.emplace(std::forward<TT>(item));
+    }
+  }
+
+template<typename T, typename A>
+ebpps_sample<T,A>::ebpps_sample(std::vector<T, A>&& data, optional<T>&& 
partial_item, double c, const A& allocator) :
+  allocator_(allocator),
+  c_(c),
+  partial_item_(partial_item),
+  data_(data, allocator)
+  {}
+
+template<typename T, typename A>
+auto ebpps_sample<T,A>::get_sample() const -> result_type {
+  double unused;
+  double c_frac = std::modf(c_, &unused);
+  bool include_partial = next_double() < c_frac;
+  uint32_t result_size = static_cast<uint32_t>(data_.size()) + 
(include_partial ? 1 : 0);
+
+  result_type result;
+  result.reserve(result_size);
+  std::copy(data_.begin(), data_.end(), std::back_inserter(result));
+  if (include_partial)
+    result.emplace_back(static_cast<const T&>(*partial_item_));
+  
+  return result;
+}
+
+template<typename T, typename A>
+void ebpps_sample<T,A>::downsample(double theta) {
+  if (theta >= 1.0) return;
+
+  double new_c = theta * c_;

Review Comment:
   const?



##########
sampling/include/ebpps_sample_impl.hpp:
##########
@@ -0,0 +1,535 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+#ifndef _EBPPS_SAMPLE_IMPL_HPP_
+#define _EBPPS_SAMPLE_IMPL_HPP_
+
+#include "common_defs.hpp"
+#include "conditional_forward.hpp"
+#include "ebpps_sample.hpp"
+#include "serde.hpp"
+
+#include <cmath>
+#include <string>
+#include <sstream>
+
+namespace datasketches {
+
+template<typename T, typename A>
+ebpps_sample<T,A>::ebpps_sample(uint32_t reserved_size, const A& allocator) :
+  allocator_(allocator),
+  c_(0.0),
+  partial_item_(),
+  data_(allocator)
+  {
+    data_.reserve(reserved_size);
+  }
+
+template<typename T, typename A>
+template<typename TT>
+ebpps_sample<T,A>::ebpps_sample(TT&& item, double theta, const A& allocator) :
+  allocator_(allocator),
+  c_(theta),
+  partial_item_(),
+  data_(allocator)
+  {
+    if (theta == 1.0) {
+      data_.reserve(1);
+      data_.emplace_back(std::forward<TT>(item));
+    } else {
+      partial_item_.emplace(std::forward<TT>(item));
+    }
+  }
+
+template<typename T, typename A>
+ebpps_sample<T,A>::ebpps_sample(std::vector<T, A>&& data, optional<T>&& 
partial_item, double c, const A& allocator) :
+  allocator_(allocator),
+  c_(c),
+  partial_item_(partial_item),
+  data_(data, allocator)
+  {}
+
+template<typename T, typename A>
+auto ebpps_sample<T,A>::get_sample() const -> result_type {
+  double unused;
+  double c_frac = std::modf(c_, &unused);
+  bool include_partial = next_double() < c_frac;
+  uint32_t result_size = static_cast<uint32_t>(data_.size()) + 
(include_partial ? 1 : 0);
+
+  result_type result;
+  result.reserve(result_size);
+  std::copy(data_.begin(), data_.end(), std::back_inserter(result));
+  if (include_partial)
+    result.emplace_back(static_cast<const T&>(*partial_item_));
+  
+  return result;
+}
+
+template<typename T, typename A>
+void ebpps_sample<T,A>::downsample(double theta) {
+  if (theta >= 1.0) return;
+
+  double new_c = theta * c_;
+  double new_c_int;
+  double new_c_frac = std::modf(new_c, &new_c_int);
+  double c_int;
+  double c_frac = std::modf(c_, &c_int);
+
+  if (new_c_int == 0.0) {
+    // no full items retained
+    if (next_double() > (c_frac / c_)) {
+      swap_with_partial();
+    }
+    data_.clear();
+  } else if (new_c_int == c_int) {
+    // no items deleted
+    if (next_double() > (1 - theta * c_frac)/(1 - new_c_frac)) {
+      swap_with_partial();
+    }
+  } else {
+    if (next_double() < theta * c_frac) {
+      // subsample data in random order; last item is partial
+      // create sample size new_c_int then swap_with_partial()
+      subsample(static_cast<uint32_t>(new_c_int));
+      swap_with_partial();
+    } else {
+      // create sample size new_c_int + 1 then move_one_to_partial)
+      subsample(static_cast<uint32_t>(new_c_int) + 1);
+      move_one_to_partial();
+    }
+  }
+
+  if (new_c == new_c_int)
+    partial_item_.reset();
+
+  c_ = new_c;
+}
+
+template<typename T, typename A>
+template<typename FwdSample>
+void ebpps_sample<T,A>::merge(FwdSample&& other) {
+  double c_int;
+  double c_frac = std::modf(c_, &c_int);
+
+  double unused;
+  double other_c_frac = std::modf(other.c_, &unused);

Review Comment:
   const?



##########
sampling/include/ebpps_sketch.hpp:
##########
@@ -0,0 +1,281 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+#ifndef _EBPPS_SKETCH_HPP_
+#define _EBPPS_SKETCH_HPP_
+
+#include "common_defs.hpp"
+#include "ebpps_sample.hpp"
+#include "optional.hpp"
+#include "serde.hpp"
+
+#include <random>
+#include <string>
+#include <vector>
+
+namespace datasketches {
+
+template<typename A> using AllocU8 = typename 
std::allocator_traits<A>::template rebind_alloc<uint8_t>;
+template<typename A> using vector_u8 = std::vector<uint8_t, AllocU8<A>>;
+
+/// EBPPS sketch constants
+namespace ebpps_constants {
+  /// maximum value of parameter K
+  const uint32_t MAX_K = ((uint32_t) 1 << 31) - 2;
+}
+
+/**
+ * An implementation of an Exact and Bounded Sampling Proportional to Size 
sketch.
+ * 
+ * From: "Exact PPS Sampling with Bounded Sample Size",
+ * B. Hentschel, P. J. Haas, Y. Tian. Information Processing Letters, 2023.
+ * 
+ * This sketch samples data from a stream of items propotional to the weight 
of each item.
+ * The sample guarantees the presence of an item in the result is proportional 
to that item's
+ * portion of the total weight seen by the sketch, and returns a sample no 
larger than size k.
+ * 
+ * The sample may be smaller than k and the resulting size of the sample 
potentially includes
+ * a probabilistic component, meaning the resulting sample size is not always 
constant.
+ *
+ * @author Jon Malkin
+ */
+template<
+  typename T,
+  typename A = std::allocator<T>
+>
+class ebpps_sketch {
+
+  public:
+    static const uint32_t MAX_K = ebpps_constants::MAX_K;
+
+    /**
+     * Constructor
+     * @param k sketch size
+     * @param allocator instance of an allocator
+     */
+    explicit ebpps_sketch(uint32_t k, const A& allocator = A());
+
+    /**
+     * Updates this sketch with the given data item with the given weight.
+     * This method takes an lvalue.
+     * @param item an item from a stream of items
+     * @param weight the weight of the item
+     */
+    void update(const T& item, double weight = 1.0);
+
+    /**
+     * Updates this sketch with the given data item with the given weight.
+     * This method takes an rvalue.
+     * @param item an item from a stream of items
+     * @param weight the weight of the item
+     */
+    void update(T&& item, double weight = 1.0);
+
+    /**
+     * Merges the provided sketch into the current one.
+     * This method takes an lvalue.
+     * @param sketch the sketch to merge into the current object
+     */
+    void merge(const ebpps_sketch<T, A>& sketch);
+
+    /**
+     * Merges the provided sketch into the current one.
+     * This method takes an rvalue.
+     * @param sketch the sketch to merge into the current object
+     */
+    void merge(ebpps_sketch<T, A>&& sketch);
+
+    using result_type = typename ebpps_sample<T,A>::result_type;
+
+    /**
+     * @brief Returns a copy of the current sample, as a std::vector
+     */
+    result_type get_result() const;
+
+    /**
+     * Returns the configured maximum sample size.
+     * @return configured maximum sample size
+     */
+    inline uint32_t get_k() const;
+
+    /**
+     * Returns the number of items processed by the sketch, regardless of
+     * item weight.
+     * @return count of items processed by the sketch
+     */
+    inline uint64_t get_n() const;
+
+    /**
+     * Returns the cumulative weight of items processed by the sketch.
+     * @return cumulative weight of items seen
+     */
+    inline double get_cumulative_weight() const;
+
+    /**
+     * Returns the number of samples currently in the sketch

Review Comment:
   I think this needs better explanation. Fractional number of samples is 
confusing.



##########
sampling/include/ebpps_sketch_impl.hpp:
##########
@@ -0,0 +1,533 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+#ifndef _EBPPS_SKETCH_IMPL_HPP_
+#define _EBPPS_SKETCH_IMPL_HPP_
+
+#include <memory>
+#include <sstream>
+#include <cmath>
+#include <random>
+#include <algorithm>
+#include <stdexcept>
+#include <utility>
+
+#include "ebpps_sketch.hpp"
+
+namespace datasketches {
+
+template<typename T, typename A>
+ebpps_sketch<T, A>::ebpps_sketch(uint32_t k, const A& allocator) :
+  allocator_(allocator),
+  k_(k),
+  n_(0),
+  cumulative_wt_(0.0),
+  wt_max_(0.0),
+  rho_(1.0),
+  sample_(check_k(k), allocator)
+  {}
+
+template<typename T, typename A>
+ebpps_sketch<T,A>::ebpps_sketch(uint32_t k, uint64_t n, double cumulative_wt,
+                                double wt_max, double rho,
+                                ebpps_sample<T,A>&& sample, const A& 
allocator) :
+  allocator_(allocator),
+  k_(k),
+  n_(n),
+  cumulative_wt_(cumulative_wt),
+  wt_max_(wt_max),
+  rho_(rho),
+  sample_(sample)
+  {}
+
+template<typename T, typename A>
+uint32_t ebpps_sketch<T, A>::get_k() const {
+  return k_;
+}
+
+template<typename T, typename A>
+uint64_t ebpps_sketch<T, A>::get_n() const {
+  return n_;
+}
+
+template<typename T, typename A>
+double ebpps_sketch<T, A>::get_c() const {
+  return sample_.get_c();
+}
+
+template<typename T, typename A>
+double ebpps_sketch<T, A>::get_cumulative_weight() const {
+  return cumulative_wt_;
+}
+
+template<typename T, typename A>
+bool ebpps_sketch<T, A>::is_empty() const {
+  return n_ == 0;
+}
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::reset() {
+  n_ = 0;
+  cumulative_wt_ = 0.0;
+  wt_max_ = 0.0;
+  rho_ = 1.0;
+  sample_.reset();
+}
+
+template<typename T, typename A>
+string<A> ebpps_sketch<T, A>::to_string() const {
+  // Using a temporary stream for implementation here does not comply with 
AllocatorAwareContainer requirements.
+  // The stream does not support passing an allocator instance, and 
alternatives are complicated.
+  std::ostringstream os;
+  os << "### EBPPS Sketch SUMMARY:" << std::endl;
+  os << "   k            : " << k_ << std::endl;
+  os << "   n            : " << n_ << std::endl;
+  os << "   cum. weight  : " << cumulative_wt_ << std::endl;
+  os << "   wt_mac       : " << wt_max_ << std::endl;
+  os << "   rho          : " << rho_ << std::endl;
+  os << "   C            : " << sample_.get_c() << std::endl;
+  os << "### END SKETCH SUMMARY" << std::endl;
+  return string<A>(os.str().c_str(), allocator_);
+}
+
+template<typename T, typename A>
+string<A> ebpps_sketch<T, A>::items_to_string() const {
+  // Using a temporary stream for implementation here does not comply with 
AllocatorAwareContainer requirements.
+  // The stream does not support passing an allocator instance, and 
alternatives are complicated.
+  std::ostringstream os;
+  os << "### Sketch Items" << std::endl;
+  os << sample_.to_string(); // assumes std::endl at end
+  return string<A>(os.str().c_str(), allocator_);
+}
+
+template<typename T, typename A>
+A ebpps_sketch<T, A>::get_allocator() const {
+  return allocator_;
+}
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::update(const T& item, double weight) {
+  return internal_update(item, weight);
+}
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::update(T&& item, double weight) {
+  return internal_update(std::move(item), weight);
+}
+
+template<typename T, typename A>
+template<typename FwdItem>
+void ebpps_sketch<T, A>::internal_update(FwdItem&& item, double weight) {
+  if (weight < 0.0 || std::isnan(weight) || std::isinf(weight)) {
+    throw std::invalid_argument("Item weights must be nonnegative and finite. 
Found: "
+                                + std::to_string(weight));
+  } else if (weight == 0.0) {
+    return;
+  }
+
+  double new_cum_wt = cumulative_wt_ + weight;

Review Comment:
   const?



##########
sampling/include/ebpps_sample_impl.hpp:
##########
@@ -0,0 +1,535 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+#ifndef _EBPPS_SAMPLE_IMPL_HPP_
+#define _EBPPS_SAMPLE_IMPL_HPP_
+
+#include "common_defs.hpp"
+#include "conditional_forward.hpp"
+#include "ebpps_sample.hpp"
+#include "serde.hpp"
+
+#include <cmath>
+#include <string>
+#include <sstream>
+
+namespace datasketches {
+
+template<typename T, typename A>
+ebpps_sample<T,A>::ebpps_sample(uint32_t reserved_size, const A& allocator) :
+  allocator_(allocator),
+  c_(0.0),
+  partial_item_(),
+  data_(allocator)
+  {
+    data_.reserve(reserved_size);
+  }
+
+template<typename T, typename A>
+template<typename TT>
+ebpps_sample<T,A>::ebpps_sample(TT&& item, double theta, const A& allocator) :
+  allocator_(allocator),
+  c_(theta),
+  partial_item_(),
+  data_(allocator)
+  {
+    if (theta == 1.0) {
+      data_.reserve(1);
+      data_.emplace_back(std::forward<TT>(item));
+    } else {
+      partial_item_.emplace(std::forward<TT>(item));
+    }
+  }
+
+template<typename T, typename A>
+ebpps_sample<T,A>::ebpps_sample(std::vector<T, A>&& data, optional<T>&& 
partial_item, double c, const A& allocator) :
+  allocator_(allocator),
+  c_(c),
+  partial_item_(partial_item),
+  data_(data, allocator)
+  {}
+
+template<typename T, typename A>
+auto ebpps_sample<T,A>::get_sample() const -> result_type {
+  double unused;
+  double c_frac = std::modf(c_, &unused);
+  bool include_partial = next_double() < c_frac;
+  uint32_t result_size = static_cast<uint32_t>(data_.size()) + 
(include_partial ? 1 : 0);
+
+  result_type result;
+  result.reserve(result_size);
+  std::copy(data_.begin(), data_.end(), std::back_inserter(result));
+  if (include_partial)
+    result.emplace_back(static_cast<const T&>(*partial_item_));
+  
+  return result;
+}
+
+template<typename T, typename A>
+void ebpps_sample<T,A>::downsample(double theta) {
+  if (theta >= 1.0) return;
+
+  double new_c = theta * c_;
+  double new_c_int;
+  double new_c_frac = std::modf(new_c, &new_c_int);
+  double c_int;
+  double c_frac = std::modf(c_, &c_int);
+
+  if (new_c_int == 0.0) {
+    // no full items retained
+    if (next_double() > (c_frac / c_)) {
+      swap_with_partial();
+    }
+    data_.clear();
+  } else if (new_c_int == c_int) {
+    // no items deleted
+    if (next_double() > (1 - theta * c_frac)/(1 - new_c_frac)) {
+      swap_with_partial();
+    }
+  } else {
+    if (next_double() < theta * c_frac) {
+      // subsample data in random order; last item is partial
+      // create sample size new_c_int then swap_with_partial()
+      subsample(static_cast<uint32_t>(new_c_int));
+      swap_with_partial();
+    } else {
+      // create sample size new_c_int + 1 then move_one_to_partial)
+      subsample(static_cast<uint32_t>(new_c_int) + 1);
+      move_one_to_partial();
+    }
+  }
+
+  if (new_c == new_c_int)
+    partial_item_.reset();
+
+  c_ = new_c;
+}
+
+template<typename T, typename A>
+template<typename FwdSample>
+void ebpps_sample<T,A>::merge(FwdSample&& other) {
+  double c_int;
+  double c_frac = std::modf(c_, &c_int);
+
+  double unused;
+  double other_c_frac = std::modf(other.c_, &unused);
+
+  // update c_ here but do NOT recompute fractional part yet
+  c_ += other.c_;
+
+  for (uint32_t i = 0; i < other.data_.size(); ++i)
+    data_.emplace_back(conditional_forward<FwdSample>(other.data_[i]));
+
+  // This modifies the original algorithm slightly due to numeric
+  // precision issues. Specifically, the test if c_frac + other_c_frac == 1.0
+  // happens before tests for < 1.0 or > 1.0 and can also be triggered
+  // if c_ == floor(c_) (the updated value of c_, not the input).
+  //
+  // We can still run into issues where c_frac + other_c_frac == epsilon
+  // and the first case would have ideally triggered. As a result, we must
+  // check if the partial item exists before adding to the data_ vector.
+
+  if (c_frac == 0.0 && other_c_frac == 0.0) {
+    partial_item_.reset();
+  } else if (c_frac + other_c_frac == 1.0 || c_ == std::floor(c_)) {
+    if (next_double() <= c_frac) {
+      if (partial_item_)
+        data_.emplace_back(std::move(*partial_item_));
+    } else {
+      if (other.partial_item_)
+        
data_.emplace_back(conditional_forward<FwdSample>(*other.partial_item_));
+    }
+    partial_item_.reset();
+  } else if (c_frac + other_c_frac < 1.0) {
+    if (next_double() > c_frac / (c_frac + other_c_frac)) {
+      set_partial(conditional_forward<FwdSample>(*other.partial_item_));
+    }
+  } else { // c_frac + other_c_frac > 1
+    if (next_double() <= (1 - c_frac) / ((1 - c_frac) + (1 - other_c_frac))) {
+      data_.emplace_back(conditional_forward<FwdSample>(*other.partial_item_));
+    } else {
+      data_.emplace_back(std::move(*partial_item_));
+      partial_item_.reset();
+      set_partial(conditional_forward<FwdSample>(*other.partial_item_));
+    }
+  }
+}
+
+template<typename T, typename A>
+string<A> ebpps_sample<T ,A>::to_string() const {
+  std::ostringstream oss;
+  oss << "   sample:" << std::endl;
+  uint32_t idx = 0;
+  for (const T& item : data_)
+    oss << "\t" << idx++ << ":\t" << item << std::endl;
+  oss << "   partial: " << (bool(partial_item_) ? (*partial_item_) : "NULL") 
<< std::endl;
+
+  return oss.str();
+}
+
+template<typename T, typename A>
+void ebpps_sample<T,A>::subsample(uint32_t num_samples) {
+  // we can perform a Fisher-Yates style shuffle, stopping after
+  // num_samples points since subsequent swaps would only be
+  // between items after num_samples. This is valid since a
+  // point from anywhere in the initial array would be eligible
+  // to end up in the final subsample.
+
+  if (num_samples == data_.size()) return;
+
+  auto erase_start = data_.begin();
+  uint32_t data_len = static_cast<uint32_t>(data_.size());
+  for (uint32_t i = 0; i < num_samples; ++i, ++erase_start) {
+    uint32_t j = i + random_idx(data_len - i);
+    std::swap(data_[i], data_[j]);
+  }
+
+  data_.erase(erase_start, data_.end());
+}
+
+template<typename T, typename A>
+template<typename FwdItem>
+void ebpps_sample<T,A>::set_partial(FwdItem&& item) {
+  if (partial_item_)
+    *partial_item_ = conditional_forward<FwdItem>(item);
+  else
+    partial_item_.emplace(conditional_forward<FwdItem>(item));
+}
+
+template<typename T, typename A>
+void ebpps_sample<T,A>::move_one_to_partial() {
+  size_t idx = random_idx(static_cast<uint32_t>(data_.size()));
+  // swap selected item to end so we can delete it easily
+  size_t last_idx = data_.size() - 1;
+  if (idx != last_idx) {
+    std::swap(data_[idx], data_[last_idx]);
+  }
+
+  set_partial(std::move(data_[last_idx]));
+
+  data_.pop_back();
+}
+
+template<typename T, typename A>
+void ebpps_sample<T,A>::swap_with_partial() {
+  if (partial_item_) {
+    size_t idx = random_idx(static_cast<uint32_t>(data_.size()));

Review Comment:
   const



##########
sampling/include/ebpps_sketch_impl.hpp:
##########
@@ -0,0 +1,533 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+#ifndef _EBPPS_SKETCH_IMPL_HPP_
+#define _EBPPS_SKETCH_IMPL_HPP_
+
+#include <memory>
+#include <sstream>
+#include <cmath>
+#include <random>
+#include <algorithm>
+#include <stdexcept>
+#include <utility>
+
+#include "ebpps_sketch.hpp"
+
+namespace datasketches {
+
+template<typename T, typename A>
+ebpps_sketch<T, A>::ebpps_sketch(uint32_t k, const A& allocator) :
+  allocator_(allocator),
+  k_(k),
+  n_(0),
+  cumulative_wt_(0.0),
+  wt_max_(0.0),
+  rho_(1.0),
+  sample_(check_k(k), allocator)
+  {}
+
+template<typename T, typename A>
+ebpps_sketch<T,A>::ebpps_sketch(uint32_t k, uint64_t n, double cumulative_wt,
+                                double wt_max, double rho,
+                                ebpps_sample<T,A>&& sample, const A& 
allocator) :
+  allocator_(allocator),
+  k_(k),
+  n_(n),
+  cumulative_wt_(cumulative_wt),
+  wt_max_(wt_max),
+  rho_(rho),
+  sample_(sample)
+  {}
+
+template<typename T, typename A>
+uint32_t ebpps_sketch<T, A>::get_k() const {
+  return k_;
+}
+
+template<typename T, typename A>
+uint64_t ebpps_sketch<T, A>::get_n() const {
+  return n_;
+}
+
+template<typename T, typename A>
+double ebpps_sketch<T, A>::get_c() const {
+  return sample_.get_c();
+}
+
+template<typename T, typename A>
+double ebpps_sketch<T, A>::get_cumulative_weight() const {
+  return cumulative_wt_;
+}
+
+template<typename T, typename A>
+bool ebpps_sketch<T, A>::is_empty() const {
+  return n_ == 0;
+}
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::reset() {
+  n_ = 0;
+  cumulative_wt_ = 0.0;
+  wt_max_ = 0.0;
+  rho_ = 1.0;
+  sample_.reset();
+}
+
+template<typename T, typename A>
+string<A> ebpps_sketch<T, A>::to_string() const {
+  // Using a temporary stream for implementation here does not comply with 
AllocatorAwareContainer requirements.
+  // The stream does not support passing an allocator instance, and 
alternatives are complicated.
+  std::ostringstream os;
+  os << "### EBPPS Sketch SUMMARY:" << std::endl;
+  os << "   k            : " << k_ << std::endl;
+  os << "   n            : " << n_ << std::endl;
+  os << "   cum. weight  : " << cumulative_wt_ << std::endl;
+  os << "   wt_mac       : " << wt_max_ << std::endl;
+  os << "   rho          : " << rho_ << std::endl;
+  os << "   C            : " << sample_.get_c() << std::endl;
+  os << "### END SKETCH SUMMARY" << std::endl;
+  return string<A>(os.str().c_str(), allocator_);
+}
+
+template<typename T, typename A>
+string<A> ebpps_sketch<T, A>::items_to_string() const {
+  // Using a temporary stream for implementation here does not comply with 
AllocatorAwareContainer requirements.
+  // The stream does not support passing an allocator instance, and 
alternatives are complicated.
+  std::ostringstream os;
+  os << "### Sketch Items" << std::endl;
+  os << sample_.to_string(); // assumes std::endl at end
+  return string<A>(os.str().c_str(), allocator_);
+}
+
+template<typename T, typename A>
+A ebpps_sketch<T, A>::get_allocator() const {
+  return allocator_;
+}
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::update(const T& item, double weight) {
+  return internal_update(item, weight);
+}
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::update(T&& item, double weight) {
+  return internal_update(std::move(item), weight);
+}
+
+template<typename T, typename A>
+template<typename FwdItem>
+void ebpps_sketch<T, A>::internal_update(FwdItem&& item, double weight) {
+  if (weight < 0.0 || std::isnan(weight) || std::isinf(weight)) {
+    throw std::invalid_argument("Item weights must be nonnegative and finite. 
Found: "
+                                + std::to_string(weight));
+  } else if (weight == 0.0) {
+    return;
+  }
+
+  double new_cum_wt = cumulative_wt_ + weight;
+  double new_wt_max = std::max(wt_max_, weight);
+  double new_rho = std::min(1.0 / new_wt_max, k_ / new_cum_wt);
+
+  if (cumulative_wt_ > 0.0)
+    sample_.downsample(new_rho / rho_);
+  
+  ebpps_sample<T,A> tmp(conditional_forward<FwdItem>(item), new_rho * weight, 
allocator_);
+
+  sample_.merge(tmp);
+
+  cumulative_wt_ = new_cum_wt;
+  wt_max_ = new_wt_max;
+  rho_ = new_rho;
+  ++n_;
+}
+
+template<typename T, typename A>
+auto ebpps_sketch<T,A>::get_result() const -> result_type {
+  return sample_.get_sample();
+}
+
+/* Merging
+ * There is a trivial merge algorithm that involves downsampling each sketch A 
and B
+ * as A.cum_wt / (A.cum_wt + B.cum_wt) and B.cum_wt / (A.cum_wt + B.cum_wt),
+ * respectively. That merge does preserve first-order probabilities, 
specifically
+ * the probability proportional to size property, and like all other known 
merge
+ * algorithms distorts second-order probabilities (co-occurrences). There are
+ * pathological cases, most obvious with k=2 and A.cum_wt == B.cum_wt where 
that
+ * approach will always take exactly 1 item from A and 1 from B, meaning the
+ * co-occurrence rate for two items from either sketch is guaranteed to be 0.0.
+ * 
+ * With EBPPS, once an item is accepted into the sketch we no longer need to
+ * track the item's weight: All accepted items are treated equally. As a 
result, we
+ * can take inspiration from the reservoir sampling merge in the 
datasketches-java
+ * library. We need to merge the smaller sketch into the larger one, swapping 
as
+ * needed to ensure that, at which point we simply call update() with the items
+ * in the smaller sketch as long as we adjust the weight appropriately.
+ * Merging smaller into larger is essential to ensure that no item has a
+ * contribution to C > 1.0.
+ */
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::merge(ebpps_sketch<T, A>&& sk) {
+  if (sk.get_cumulative_weight() == 0.0) return;
+  else if (sk.get_cumulative_weight() > get_cumulative_weight()) {
+    // need to swap this with sk to merge smaller into larger
+    std::swap(*this, sk);
+  }
+
+  internal_merge(sk);
+}
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::merge(const ebpps_sketch<T, A>& sk) {
+  if (sk.get_cumulative_weight() == 0.0) return;
+  else if (sk.get_cumulative_weight() > get_cumulative_weight()) {
+    // need to swap this with sk to merge, so make a copy, swap,
+    // and use that to merge
+    ebpps_sketch sk_copy(sk);
+    swap(*this, sk_copy);
+    internal_merge(sk_copy);
+  } else {
+    internal_merge(sk);
+  }
+}
+
+template<typename T, typename A>
+template<typename O>
+void ebpps_sketch<T, A>::internal_merge(O&& sk) {
+  // assumes that sk.cumulative_wt_ <= cumulative_wt_,
+  // which must be checked before calling this
+  
+  const ebpps_sample<T,A>& other_sample = sk.sample_;
+
+  double final_cum_wt = cumulative_wt_ + sk.cumulative_wt_;
+  double new_wt_max = std::max(wt_max_, sk.wt_max_);
+  k_ = std::min(k_, sk.k_);
+  uint64_t new_n = n_ + sk.n_;
+
+  // Insert sk's items with the cumulative weight
+  // split between the input items. We repeat the same process
+  // for full items and the partial item, scaling the input
+  // weight appropriately.
+  // We handle all C input items, meaning we always process
+  // the partial item using a scaled down weight.
+  // Handling the partial item by probabilistically including
+  // it as a full item would be correct on average but would
+  // introduce bias for any specific merge operation.
+  double avg_wt = sk.get_cumulative_weight() / sk.get_c();
+  auto items = other_sample.get_full_items();
+  for (size_t i = 0; i < items.size(); ++i) {
+    // new_wt_max is pre-computed
+    double new_cum_wt = cumulative_wt_ + avg_wt;
+    double new_rho = std::min(1.0 / new_wt_max, k_ / new_cum_wt);
+
+    if (cumulative_wt_ > 0.0)
+      sample_.downsample(new_rho / rho_);
+  
+    ebpps_sample<T,A> tmp(conditional_forward<O>(items[i]), new_rho * avg_wt, 
allocator_);
+
+    sample_.merge(tmp);
+
+    cumulative_wt_ = new_cum_wt;
+    rho_ = new_rho;
+  }
+
+  // insert partial item with weight scaled by the fractional part of C
+  if (other_sample.has_partial_item()) {
+    double unused;
+    double other_c_frac = std::modf(other_sample.get_c(), &unused);
+
+    double new_cum_wt = cumulative_wt_ + (other_c_frac * avg_wt);
+    double new_rho = std::min(1.0 / new_wt_max, k_ / new_cum_wt);
+
+    if (cumulative_wt_ > 0.0)
+      sample_.downsample(new_rho / rho_);
+  
+    ebpps_sample<T,A> 
tmp(conditional_forward<O>(other_sample.get_partial_item()), new_rho * 
other_c_frac * avg_wt, allocator_);
+
+    sample_.merge(tmp);
+
+    cumulative_wt_ = new_cum_wt;
+    rho_ = new_rho;
+  }
+
+  // avoid numeric issues by setting cumulative weight to the
+  // pre-computed value
+  cumulative_wt_ = final_cum_wt;
+  n_ = new_n;
+}
+
+/*
+ * An empty sketch requires 8 bytes.
+ *
+ * <pre>
+ * Long || Start Byte Adr:
+ * Adr:
+ *      ||       0        |    1   |    2   |    3   |    4   |    5   |    6  
 |    7   |
+ *  0   || Preamble_Longs | SerVer | FamID  |  Flags |---------Max Res. Size 
(K)---------|
+ * </pre>
+ *
+ * A non-empty sketch requires 40 bytes of preamble. C looks like part of
+ * the preamble but is serialized as part of the internal sample state.
+ *
+ * The count of items seen is not used but preserved as the value seems like a 
useful
+ * count to track.
+ * 
+ * <pre>
+ * Long || Start Byte Adr:
+ * Adr:
+ *      ||       0        |    1   |    2   |    3   |    4   |    5   |    6  
 |    7   |
+ *  0   || Preamble_Longs | SerVer | FamID  |  Flags |---------Max Res. Size 
(K)---------|
+ *
+ *      ||       8        |    9   |   10   |   11   |   12   |   13   |   14  
 |   15   |
+ *  1   ||---------------------------Items Seen Count 
(N)--------------------------------|
+ *
+ *      ||      16        |   17   |   18   |   19   |   20   |   21   |   22  
 |   23   |
+ *  2   ||----------------------------Cumulative 
Weight----------------------------------|
+ *
+ *      ||      24        |   25   |   26   |   27   |   28   |   29   |   30  
 |   31   |
+ *  3   ||-----------------------------Max Item 
Weight-----------------------------------|
+ *
+ *      ||      32        |   33   |   34   |   35   |   36   |   37   |   38  
 |   39   |
+ *  4   
||----------------------------------Rho------------------------------------------|
+ *
+ *      ||      40        |   41   |   42   |   43   |   44   |   45   |   46  
 |   47   |
+ *  5   
||-----------------------------------C-------------------------------------------|
+ *
+ *      ||      40+                      |
+ *  6+  ||  {Items Array}                |
+ *      ||  {Optional Item (if needed)}  |
+ * </pre>
+ */
+
+template<typename T, typename A>
+template<typename SerDe>
+size_t ebpps_sketch<T, A>::get_serialized_size_bytes(const SerDe& sd) const {
+  if (is_empty()) { return PREAMBLE_LONGS_EMPTY << 3; }
+  return (PREAMBLE_LONGS_FULL << 3) + sample_.get_serialized_size_bytes(sd);
+}
+
+template<typename T, typename A>
+template<typename SerDe>
+auto ebpps_sketch<T,A>::serialize(unsigned header_size_bytes, const SerDe& sd) 
const -> vector_bytes {
+  bool empty = is_empty();
+  uint8_t prelongs = (empty ? PREAMBLE_LONGS_EMPTY : PREAMBLE_LONGS_FULL);
+
+  const size_t size = header_size_bytes + (prelongs << 3) + 
sample_.get_serialized_size_bytes(sd);
+  vector_bytes bytes(size, 0, allocator_);
+  uint8_t* ptr = bytes.data() + header_size_bytes;
+  uint8_t* end_ptr = ptr + size;
+
+  uint8_t flags = 0;
+  if (empty) {
+    flags |= EMPTY_FLAG_MASK;
+  } else {
+    flags |= sample_.has_partial_item() ? HAS_PARTIAL_ITEM_MASK : 0;
+  }
+
+  // first prelong
+  uint8_t ser_ver = SER_VER;
+  uint8_t family = FAMILY_ID;
+  ptr += copy_to_mem(prelongs, ptr);
+  ptr += copy_to_mem(ser_ver, ptr);
+  ptr += copy_to_mem(family, ptr);
+  ptr += copy_to_mem(flags, ptr);
+  ptr += copy_to_mem(k_, ptr);
+
+  if (!empty) {
+    // remaining preamble
+    ptr += copy_to_mem(n_, ptr);
+    ptr += copy_to_mem(cumulative_wt_, ptr);
+    ptr += copy_to_mem(wt_max_, ptr);
+    ptr += copy_to_mem(rho_, ptr);
+    ptr += sample_.serialize(ptr, end_ptr, sd);
+  }
+
+  return bytes;
+}
+
+template<typename T, typename A>
+template<typename SerDe>
+void ebpps_sketch<T,A>::serialize(std::ostream& os, const SerDe& sd) const {
+  const bool empty = is_empty();

Review Comment:
   and this



##########
sampling/include/ebpps_sketch_impl.hpp:
##########
@@ -0,0 +1,533 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+#ifndef _EBPPS_SKETCH_IMPL_HPP_
+#define _EBPPS_SKETCH_IMPL_HPP_
+
+#include <memory>
+#include <sstream>
+#include <cmath>
+#include <random>
+#include <algorithm>
+#include <stdexcept>
+#include <utility>
+
+#include "ebpps_sketch.hpp"
+
+namespace datasketches {
+
+template<typename T, typename A>
+ebpps_sketch<T, A>::ebpps_sketch(uint32_t k, const A& allocator) :
+  allocator_(allocator),
+  k_(k),
+  n_(0),
+  cumulative_wt_(0.0),
+  wt_max_(0.0),
+  rho_(1.0),
+  sample_(check_k(k), allocator)
+  {}
+
+template<typename T, typename A>
+ebpps_sketch<T,A>::ebpps_sketch(uint32_t k, uint64_t n, double cumulative_wt,
+                                double wt_max, double rho,
+                                ebpps_sample<T,A>&& sample, const A& 
allocator) :
+  allocator_(allocator),
+  k_(k),
+  n_(n),
+  cumulative_wt_(cumulative_wt),
+  wt_max_(wt_max),
+  rho_(rho),
+  sample_(sample)
+  {}
+
+template<typename T, typename A>
+uint32_t ebpps_sketch<T, A>::get_k() const {
+  return k_;
+}
+
+template<typename T, typename A>
+uint64_t ebpps_sketch<T, A>::get_n() const {
+  return n_;
+}
+
+template<typename T, typename A>
+double ebpps_sketch<T, A>::get_c() const {
+  return sample_.get_c();
+}
+
+template<typename T, typename A>
+double ebpps_sketch<T, A>::get_cumulative_weight() const {
+  return cumulative_wt_;
+}
+
+template<typename T, typename A>
+bool ebpps_sketch<T, A>::is_empty() const {
+  return n_ == 0;
+}
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::reset() {
+  n_ = 0;
+  cumulative_wt_ = 0.0;
+  wt_max_ = 0.0;
+  rho_ = 1.0;
+  sample_.reset();
+}
+
+template<typename T, typename A>
+string<A> ebpps_sketch<T, A>::to_string() const {
+  // Using a temporary stream for implementation here does not comply with 
AllocatorAwareContainer requirements.
+  // The stream does not support passing an allocator instance, and 
alternatives are complicated.
+  std::ostringstream os;
+  os << "### EBPPS Sketch SUMMARY:" << std::endl;
+  os << "   k            : " << k_ << std::endl;
+  os << "   n            : " << n_ << std::endl;
+  os << "   cum. weight  : " << cumulative_wt_ << std::endl;
+  os << "   wt_mac       : " << wt_max_ << std::endl;
+  os << "   rho          : " << rho_ << std::endl;
+  os << "   C            : " << sample_.get_c() << std::endl;
+  os << "### END SKETCH SUMMARY" << std::endl;
+  return string<A>(os.str().c_str(), allocator_);
+}
+
+template<typename T, typename A>
+string<A> ebpps_sketch<T, A>::items_to_string() const {
+  // Using a temporary stream for implementation here does not comply with 
AllocatorAwareContainer requirements.
+  // The stream does not support passing an allocator instance, and 
alternatives are complicated.
+  std::ostringstream os;
+  os << "### Sketch Items" << std::endl;
+  os << sample_.to_string(); // assumes std::endl at end
+  return string<A>(os.str().c_str(), allocator_);
+}
+
+template<typename T, typename A>
+A ebpps_sketch<T, A>::get_allocator() const {
+  return allocator_;
+}
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::update(const T& item, double weight) {
+  return internal_update(item, weight);
+}
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::update(T&& item, double weight) {
+  return internal_update(std::move(item), weight);
+}
+
+template<typename T, typename A>
+template<typename FwdItem>
+void ebpps_sketch<T, A>::internal_update(FwdItem&& item, double weight) {
+  if (weight < 0.0 || std::isnan(weight) || std::isinf(weight)) {
+    throw std::invalid_argument("Item weights must be nonnegative and finite. 
Found: "
+                                + std::to_string(weight));
+  } else if (weight == 0.0) {
+    return;
+  }
+
+  double new_cum_wt = cumulative_wt_ + weight;
+  double new_wt_max = std::max(wt_max_, weight);
+  double new_rho = std::min(1.0 / new_wt_max, k_ / new_cum_wt);
+
+  if (cumulative_wt_ > 0.0)
+    sample_.downsample(new_rho / rho_);
+  
+  ebpps_sample<T,A> tmp(conditional_forward<FwdItem>(item), new_rho * weight, 
allocator_);
+
+  sample_.merge(tmp);
+
+  cumulative_wt_ = new_cum_wt;
+  wt_max_ = new_wt_max;
+  rho_ = new_rho;
+  ++n_;
+}
+
+template<typename T, typename A>
+auto ebpps_sketch<T,A>::get_result() const -> result_type {
+  return sample_.get_sample();
+}
+
+/* Merging
+ * There is a trivial merge algorithm that involves downsampling each sketch A 
and B
+ * as A.cum_wt / (A.cum_wt + B.cum_wt) and B.cum_wt / (A.cum_wt + B.cum_wt),
+ * respectively. That merge does preserve first-order probabilities, 
specifically
+ * the probability proportional to size property, and like all other known 
merge
+ * algorithms distorts second-order probabilities (co-occurrences). There are
+ * pathological cases, most obvious with k=2 and A.cum_wt == B.cum_wt where 
that
+ * approach will always take exactly 1 item from A and 1 from B, meaning the
+ * co-occurrence rate for two items from either sketch is guaranteed to be 0.0.
+ * 
+ * With EBPPS, once an item is accepted into the sketch we no longer need to
+ * track the item's weight: All accepted items are treated equally. As a 
result, we
+ * can take inspiration from the reservoir sampling merge in the 
datasketches-java
+ * library. We need to merge the smaller sketch into the larger one, swapping 
as
+ * needed to ensure that, at which point we simply call update() with the items
+ * in the smaller sketch as long as we adjust the weight appropriately.
+ * Merging smaller into larger is essential to ensure that no item has a
+ * contribution to C > 1.0.
+ */
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::merge(ebpps_sketch<T, A>&& sk) {
+  if (sk.get_cumulative_weight() == 0.0) return;
+  else if (sk.get_cumulative_weight() > get_cumulative_weight()) {
+    // need to swap this with sk to merge smaller into larger
+    std::swap(*this, sk);
+  }
+
+  internal_merge(sk);
+}
+
+template<typename T, typename A>
+void ebpps_sketch<T, A>::merge(const ebpps_sketch<T, A>& sk) {
+  if (sk.get_cumulative_weight() == 0.0) return;
+  else if (sk.get_cumulative_weight() > get_cumulative_weight()) {
+    // need to swap this with sk to merge, so make a copy, swap,
+    // and use that to merge
+    ebpps_sketch sk_copy(sk);
+    swap(*this, sk_copy);
+    internal_merge(sk_copy);
+  } else {
+    internal_merge(sk);
+  }
+}
+
+template<typename T, typename A>
+template<typename O>
+void ebpps_sketch<T, A>::internal_merge(O&& sk) {
+  // assumes that sk.cumulative_wt_ <= cumulative_wt_,
+  // which must be checked before calling this
+  
+  const ebpps_sample<T,A>& other_sample = sk.sample_;
+
+  double final_cum_wt = cumulative_wt_ + sk.cumulative_wt_;
+  double new_wt_max = std::max(wt_max_, sk.wt_max_);
+  k_ = std::min(k_, sk.k_);
+  uint64_t new_n = n_ + sk.n_;
+
+  // Insert sk's items with the cumulative weight
+  // split between the input items. We repeat the same process
+  // for full items and the partial item, scaling the input
+  // weight appropriately.
+  // We handle all C input items, meaning we always process
+  // the partial item using a scaled down weight.
+  // Handling the partial item by probabilistically including
+  // it as a full item would be correct on average but would
+  // introduce bias for any specific merge operation.
+  double avg_wt = sk.get_cumulative_weight() / sk.get_c();
+  auto items = other_sample.get_full_items();
+  for (size_t i = 0; i < items.size(); ++i) {
+    // new_wt_max is pre-computed
+    double new_cum_wt = cumulative_wt_ + avg_wt;
+    double new_rho = std::min(1.0 / new_wt_max, k_ / new_cum_wt);
+
+    if (cumulative_wt_ > 0.0)
+      sample_.downsample(new_rho / rho_);
+  
+    ebpps_sample<T,A> tmp(conditional_forward<O>(items[i]), new_rho * avg_wt, 
allocator_);
+
+    sample_.merge(tmp);
+
+    cumulative_wt_ = new_cum_wt;
+    rho_ = new_rho;
+  }
+
+  // insert partial item with weight scaled by the fractional part of C
+  if (other_sample.has_partial_item()) {
+    double unused;
+    double other_c_frac = std::modf(other_sample.get_c(), &unused);
+
+    double new_cum_wt = cumulative_wt_ + (other_c_frac * avg_wt);
+    double new_rho = std::min(1.0 / new_wt_max, k_ / new_cum_wt);
+
+    if (cumulative_wt_ > 0.0)
+      sample_.downsample(new_rho / rho_);
+  
+    ebpps_sample<T,A> 
tmp(conditional_forward<O>(other_sample.get_partial_item()), new_rho * 
other_c_frac * avg_wt, allocator_);
+
+    sample_.merge(tmp);
+
+    cumulative_wt_ = new_cum_wt;
+    rho_ = new_rho;
+  }
+
+  // avoid numeric issues by setting cumulative weight to the
+  // pre-computed value
+  cumulative_wt_ = final_cum_wt;
+  n_ = new_n;
+}
+
+/*
+ * An empty sketch requires 8 bytes.
+ *
+ * <pre>
+ * Long || Start Byte Adr:
+ * Adr:
+ *      ||       0        |    1   |    2   |    3   |    4   |    5   |    6  
 |    7   |
+ *  0   || Preamble_Longs | SerVer | FamID  |  Flags |---------Max Res. Size 
(K)---------|
+ * </pre>
+ *
+ * A non-empty sketch requires 40 bytes of preamble. C looks like part of
+ * the preamble but is serialized as part of the internal sample state.
+ *
+ * The count of items seen is not used but preserved as the value seems like a 
useful
+ * count to track.
+ * 
+ * <pre>
+ * Long || Start Byte Adr:
+ * Adr:
+ *      ||       0        |    1   |    2   |    3   |    4   |    5   |    6  
 |    7   |
+ *  0   || Preamble_Longs | SerVer | FamID  |  Flags |---------Max Res. Size 
(K)---------|
+ *
+ *      ||       8        |    9   |   10   |   11   |   12   |   13   |   14  
 |   15   |
+ *  1   ||---------------------------Items Seen Count 
(N)--------------------------------|
+ *
+ *      ||      16        |   17   |   18   |   19   |   20   |   21   |   22  
 |   23   |
+ *  2   ||----------------------------Cumulative 
Weight----------------------------------|
+ *
+ *      ||      24        |   25   |   26   |   27   |   28   |   29   |   30  
 |   31   |
+ *  3   ||-----------------------------Max Item 
Weight-----------------------------------|
+ *
+ *      ||      32        |   33   |   34   |   35   |   36   |   37   |   38  
 |   39   |
+ *  4   
||----------------------------------Rho------------------------------------------|
+ *
+ *      ||      40        |   41   |   42   |   43   |   44   |   45   |   46  
 |   47   |
+ *  5   
||-----------------------------------C-------------------------------------------|
+ *
+ *      ||      40+                      |
+ *  6+  ||  {Items Array}                |
+ *      ||  {Optional Item (if needed)}  |
+ * </pre>
+ */
+
+template<typename T, typename A>
+template<typename SerDe>
+size_t ebpps_sketch<T, A>::get_serialized_size_bytes(const SerDe& sd) const {
+  if (is_empty()) { return PREAMBLE_LONGS_EMPTY << 3; }
+  return (PREAMBLE_LONGS_FULL << 3) + sample_.get_serialized_size_bytes(sd);
+}
+
+template<typename T, typename A>
+template<typename SerDe>
+auto ebpps_sketch<T,A>::serialize(unsigned header_size_bytes, const SerDe& sd) 
const -> vector_bytes {
+  bool empty = is_empty();

Review Comment:
   do we really need this copy?



-- 
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.

To unsubscribe, e-mail: commits-unsubscr...@datasketches.apache.org

For queries about this service, please contact Infrastructure at:
us...@infra.apache.org


---------------------------------------------------------------------
To unsubscribe, e-mail: commits-unsubscr...@datasketches.apache.org
For additional commands, e-mail: commits-h...@datasketches.apache.org