/src/duckdb/third_party/tdigest/t_digest.hpp

Source
/*
 * Licensed to Derrick R. Burns under one or more
 * contributor license agreements.  See the NOTICES 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.
 */

#pragma once

#include <algorithm>
#include <cfloat>
#include <cmath>
#include <queue>
#include <utility>
#include <vector>

#ifdef min
#undef min
#endif

#ifdef max
#undef max
#endif


namespace duckdb_tdigest {

using Value = double;
using Weight = double;
using Index = size_t;

const size_t kHighWater = 40000;
const double pi = 3.14159265358979323846;
class Centroid {
public:
  Centroid() : Centroid(0.0, 0.0) {
  }

  Centroid(Value mean, Weight weight) : mean_(mean), weight_(weight) {
  }

  inline Value mean() const noexcept {
    return mean_;
  }

  inline Weight weight() const noexcept {
    return weight_;
  }

  inline void add(const Centroid &c) {
    //    CHECK_GT(c.weight_, 0);
    if (weight_ != 0.0) {
      weight_ += c.weight_;
      mean_ += c.weight_ * (c.mean_ - mean_) / weight_;
    } else {
      weight_ = c.weight_;
      mean_ = c.mean_;
    }
  }

private:
  Value mean_ = 0;
  Weight weight_ = 0;
};

struct CentroidList {
  explicit CentroidList(const std::vector<Centroid> &s) : iter(s.cbegin()), end(s.cend()) {
  }
  std::vector<Centroid>::const_iterator iter;
  std::vector<Centroid>::const_iterator end;

  bool advance() {
    return ++iter != end;
  }
};

class CentroidListComparator {
public:
  CentroidListComparator() {
  }

  bool operator()(const CentroidList &left, const CentroidList &right) const {
    return left.iter->mean() > right.iter->mean();
  }
};

using CentroidListQueue = std::priority_queue<CentroidList, std::vector<CentroidList>, CentroidListComparator>;

struct CentroidComparator {
  bool operator()(const Centroid &a, const Centroid &b) const {
    return a.mean() < b.mean();
  }
};

class TDigest {
  class TDigestComparator {
  public:
    TDigestComparator() {
    }

    bool operator()(const TDigest *left, const TDigest *right) const {
      return left->totalSize() > right->totalSize();
    }
  };

  using TDigestQueue = std::priority_queue<const TDigest *, std::vector<const TDigest *>, TDigestComparator>;

public:
  TDigest() : TDigest(1000) {
  }

  explicit TDigest(Value compression) : TDigest(compression, 0) {
  }

  TDigest(Value compression, Index bufferSize) : TDigest(compression, bufferSize, 0) {
  }

  TDigest(Value compression, Index unmergedSize, Index mergedSize)
      : compression_(compression), maxProcessed_(processedSize(mergedSize, compression)),
        maxUnprocessed_(unprocessedSize(unmergedSize, compression)) {
    processed_.reserve(maxProcessed_);
    unprocessed_.reserve(maxUnprocessed_ + 1);
  }

  TDigest(std::vector<Centroid> &&processed, std::vector<Centroid> &&unprocessed, Value compression,
          Index unmergedSize, Index mergedSize)
      : TDigest(compression, unmergedSize, mergedSize) {
    processed_ = std::move(processed);
    unprocessed_ = std::move(unprocessed);

    processedWeight_ = weight(processed_);
    unprocessedWeight_ = weight(unprocessed_);
    if (!processed_.empty()) {
      min_ = std::min(min_, processed_[0].mean());
      max_ = std::max(max_, (processed_.cend() - 1)->mean());
    }
    updateCumulative();
  }

  static Weight weight(std::vector<Centroid> &centroids) noexcept {
    Weight w = 0.0;
    for (auto centroid : centroids) {
      w += centroid.weight();
    }
    return w;
  }

  TDigest &operator=(TDigest &&o) {
    compression_ = o.compression_;
    maxProcessed_ = o.maxProcessed_;
    maxUnprocessed_ = o.maxUnprocessed_;
    processedWeight_ = o.processedWeight_;
    unprocessedWeight_ = o.unprocessedWeight_;
    processed_ = std::move(o.processed_);
    unprocessed_ = std::move(o.unprocessed_);
    cumulative_ = std::move(o.cumulative_);
    min_ = o.min_;
    max_ = o.max_;
    return *this;
  }

  TDigest(TDigest &&o)
      : TDigest(std::move(o.processed_), std::move(o.unprocessed_), o.compression_, o.maxUnprocessed_,
                o.maxProcessed_) {
  }

  static inline Index processedSize(Index size, Value compression) noexcept {
    return (size == 0) ? static_cast<Index>(2 * std::ceil(compression)) : size;
  }

  static inline Index unprocessedSize(Index size, Value compression) noexcept {
    return (size == 0) ? static_cast<Index>(8 * std::ceil(compression)) : size;
  }

  // merge in another t-digest
  inline void merge(const TDigest *other) {
    std::vector<const TDigest *> others {other};
    add(others.cbegin(), others.cend());
  }

  const std::vector<Centroid> &processed() const {
    return processed_;
  }

  const std::vector<Centroid> &unprocessed() const {
    return unprocessed_;
  }

  Index maxUnprocessed() const {
    return maxUnprocessed_;
  }

  Index maxProcessed() const {
    return maxProcessed_;
  }

  Value min() const {
    return min_;
  }

  Value max() const {
    return max_;
  }

  // override the tracked minimum/maximum - used when reconstructing a serialized digest
  void setMinMax(Value min, Value max) {
    min_ = min;
    max_ = max;
  }

  inline void add(std::vector<const TDigest *> digests) {
    add(digests.cbegin(), digests.cend());
  }

  // merge in a vector of tdigests in the most efficient manner possible
  // in constant space
  // works for any value of kHighWater
  void add(std::vector<const TDigest *>::const_iterator iter, std::vector<const TDigest *>::const_iterator end) {
    if (iter != end) {
      auto size = std::distance(iter, end);
      TDigestQueue pq(TDigestComparator {});
      for (; iter != end; iter++) {
        pq.push((*iter));
      }
      std::vector<const TDigest *> batch;
      batch.reserve(size_t(size));

      size_t totalSize = 0;
      while (!pq.empty()) {
        auto td = pq.top();
        batch.push_back(td);
        pq.pop();
        totalSize += td->totalSize();
        if (totalSize >= kHighWater || pq.empty()) {
          mergeProcessed(batch);
          mergeUnprocessed(batch);
          processIfNecessary();
          batch.clear();
          totalSize = 0;
        }
      }
      updateCumulative();
    }
  }

  Weight processedWeight() const {
    return processedWeight_;
  }

  Weight unprocessedWeight() const {
    return unprocessedWeight_;
  }

  bool haveUnprocessed() const {
    return unprocessed_.size() > 0;
  }

  size_t totalSize() const {
    return processed_.size() + unprocessed_.size();
  }

  long totalWeight() const {
    return static_cast<long>(processedWeight_ + unprocessedWeight_);
  }

  // return the cdf on the t-digest
  Value cdf(Value x) {
    if (haveUnprocessed() || isDirty()) {
      process();
    }
    return cdfProcessed(x);
  }

  bool isDirty() {
    return processed_.size() > maxProcessed_ || unprocessed_.size() > maxUnprocessed_;
  }

  // return the cdf on the processed values
  Value cdfProcessed(Value x) const {
    if (processed_.empty()) {
      // no data to examin_e

      return 0.0;
    } else if (processed_.size() == 1) {
      // exactly one centroid, should have max_==min_
      auto width = max_ - min_;
      if (x < min_) {
        return 0.0;
      } else if (x > max_) {
        return 1.0;
      } else if (x - min_ <= width) {
        // min_ and max_ are too close together to do any viable interpolation
        return 0.5;
      } else {
        // interpolate if somehow we have weight > 0 and max_ != min_
        return (x - min_) / (max_ - min_);
      }
    } else {
      auto n = processed_.size();
      if (x <= min_) {
        return 0;
      }

      if (x >= max_) {
        return 1;
      }

      // check for the left tail
      if (x <= mean(0)) {

        // note that this is different than mean(0) > min_ ... this guarantees interpolation works
        if (mean(0) - min_ > 0) {
          return (x - min_) / (mean(0) - min_) * weight(0) / processedWeight_ / 2.0;
        } else {
          return 0;
        }
      }

      // and the right tail
      if (x >= mean(n - 1)) {
        if (max_ - mean(n - 1) > 0) {
          return 1.0 - (max_ - x) / (max_ - mean(n - 1)) * weight(n - 1) / processedWeight_ / 2.0;
        } else {
          return 1;
        }
      }

      CentroidComparator cc;
      auto iter = std::upper_bound(processed_.cbegin(), processed_.cend(), Centroid(x, 0), cc);

      auto i = size_t(std::distance(processed_.cbegin(), iter));
      auto z1 = x - (iter - 1)->mean();
      auto z2 = (iter)->mean() - x;
      return weightedAverage(cumulative_[i - 1], z2, cumulative_[i], z1) / processedWeight_;
    }
  }

  // this returns a quantile on the t-digest
  Value quantile(Value q) {
    if (haveUnprocessed() || isDirty()) {
      process();
    }
    return quantileProcessed(q);
  }

  // this returns a quantile on the currently processed values without changing the t-digest
  // the value will not represent the unprocessed values
  Value quantileProcessed(Value q) const {
    if (q < 0 || q > 1) {
      return NAN;
    }

    if (processed_.size() == 0) {
      // no sorted means no data, no way to get a quantile
      return NAN;
    } else if (processed_.size() == 1) {
      // with one data point, all quantiles lead to Rome

      return mean(0);
    }

    // we know that there are at least two sorted now
    auto n = processed_.size();

    // if values were stored in a sorted array, index would be the offset we are Weighterested in
    const auto index = q * processedWeight_;

    // at the boundaries, we return min_ or max_
    if (index <= weight(0) / 2.0) {
      return min_ + 2.0 * index / weight(0) * (mean(0) - min_);
    }

    auto iter = std::lower_bound(cumulative_.cbegin(), cumulative_.cend(), index);

    if (iter + 1 != cumulative_.cend()) {
      auto i = size_t(std::distance(cumulative_.cbegin(), iter));
      auto z1 = index - *(iter - 1);
      auto z2 = *(iter)-index;
      // LOG(INFO) << "z2 " << z2 << " index " << index << " z1 " << z1;
      return weightedAverage(mean(i - 1), z2, mean(i), z1);
    }
    auto z1 = index - processedWeight_ - weight(n - 1) / 2.0;
    auto z2 = weight(n - 1) / 2 - z1;
    return weightedAverage(mean(n - 1), z1, max_, z2);
  }

  Value compression() const {
    return compression_;
  }

  void add(Value x) {
    add(x, 1);
  }

  inline void compress() {
    process();
  }

  // add a single centroid to the unprocessed vector, processing previously unprocessed sorted if our limit has
  // been reached.
  inline bool add(Value x, Weight w) {
    if (std::isnan(x)) {
      return false;
    }
    unprocessed_.push_back(Centroid(x, w));
    unprocessedWeight_ += w;
    processIfNecessary();
    return true;
  }

  inline void add(std::vector<Centroid>::const_iterator iter, std::vector<Centroid>::const_iterator end) {
    while (iter != end) {
      const size_t diff = size_t(std::distance(iter, end));
      const size_t room = maxUnprocessed_ - unprocessed_.size();
      auto mid = iter + int64_t(std::min(diff, room));
      while (iter != mid) {
        unprocessed_.push_back(*(iter++));
      }
      if (unprocessed_.size() >= maxUnprocessed_) {
        process();
      }
    }
  }

private:
  Value compression_;

  Value min_ = std::numeric_limits<Value>::max();

  Value max_ = std::numeric_limits<Value>::min();

  Index maxProcessed_;

  Index maxUnprocessed_;

  Value processedWeight_ = 0.0;

  Value unprocessedWeight_ = 0.0;

  std::vector<Centroid> processed_;

  std::vector<Centroid> unprocessed_;

  std::vector<Weight> cumulative_;

  // return mean of i-th centroid
  inline Value mean(size_t i) const noexcept {
    return processed_[i].mean();
  }

  // return weight of i-th centroid
  inline Weight weight(size_t i) const noexcept {
    return processed_[i].weight();
  }

  // append all unprocessed centroids into current unprocessed vector
  void mergeUnprocessed(const std::vector<const TDigest *> &tdigests) {
    if (tdigests.size() == 0) {
      return;
    }

    size_t total = unprocessed_.size();
    for (auto &td : tdigests) {
      total += td->unprocessed_.size();
    }

    unprocessed_.reserve(total);
    for (auto &td : tdigests) {
      unprocessed_.insert(unprocessed_.end(), td->unprocessed_.cbegin(), td->unprocessed_.cend());
      unprocessedWeight_ += td->unprocessedWeight_;
    }
  }

  // merge all processed centroids together into a single sorted vector
  void mergeProcessed(const std::vector<const TDigest *> &tdigests) {
    if (tdigests.size() == 0) {
      return;
    }

    size_t total = 0;
    CentroidListQueue pq(CentroidListComparator {});
    for (auto &td : tdigests) {
      auto &sorted = td->processed_;
      auto size = sorted.size();
      if (size > 0) {
        pq.push(CentroidList(sorted));
        total += size;
        processedWeight_ += td->processedWeight_;
      }
    }
    if (total == 0) {
      return;
    }

    if (processed_.size() > 0) {
      pq.push(CentroidList(processed_));
      total += processed_.size();
    }

    std::vector<Centroid> sorted;
    sorted.reserve(total);

    while (!pq.empty()) {
      auto best = pq.top();
      pq.pop();
      sorted.push_back(*(best.iter));
      if (best.advance()) {
        pq.push(best);
      }
    }
    processed_ = std::move(sorted);
    if (processed_.size() > 0) {
      min_ = std::min(min_, processed_[0].mean());
      max_ = std::max(max_, (processed_.cend() - 1)->mean());
    }
  }

  inline void processIfNecessary() {
    if (isDirty()) {
      process();
    }
  }

  void updateCumulative() {
    const auto n = processed_.size();
    cumulative_.clear();
    cumulative_.reserve(n + 1);
    auto previous = 0.0;
    for (Index i = 0; i < n; i++) {
      auto current = weight(i);
      auto halfCurrent = current / 2.0;
      cumulative_.push_back(previous + halfCurrent);
      previous = previous + current;
    }
    cumulative_.push_back(previous);
  }

  // merges unprocessed_ centroids and processed_ centroids together and processes them
  // when complete, unprocessed_ will be empty and processed_ will have at most maxProcessed_ centroids
  inline void process() {
    CentroidComparator cc;
    std::sort(unprocessed_.begin(), unprocessed_.end(), cc);
    auto count = unprocessed_.size();
    unprocessed_.insert(unprocessed_.end(), processed_.cbegin(), processed_.cend());
    std::inplace_merge(unprocessed_.begin(), unprocessed_.begin() + int64_t(count), unprocessed_.end(), cc);

    processedWeight_ += unprocessedWeight_;
    unprocessedWeight_ = 0;
    processed_.clear();

    processed_.push_back(unprocessed_[0]);
    Weight wSoFar = unprocessed_[0].weight();
    Weight wLimit = processedWeight_ * integratedQ(1.0);

    auto end = unprocessed_.end();
    for (auto iter = unprocessed_.cbegin() + 1; iter < end; iter++) {
      auto &centroid = *iter;
      Weight projectedW = wSoFar + centroid.weight();
      if (projectedW <= wLimit) {
        wSoFar = projectedW;
        (processed_.end() - 1)->add(centroid);
      } else {
        auto k1 = integratedLocation(wSoFar / processedWeight_);
        wLimit = processedWeight_ * integratedQ(k1 + 1.0);
        wSoFar += centroid.weight();
        processed_.emplace_back(centroid);
      }
    }
    unprocessed_.clear();
    min_ = std::min(min_, processed_[0].mean());
    max_ = std::max(max_, (processed_.cend() - 1)->mean());
    updateCumulative();
  }

  inline size_t checkWeights() {
    return checkWeights(processed_, processedWeight_);
  }

  size_t checkWeights(const std::vector<Centroid> &sorted, Value total) {
    size_t badWeight = 0;
    auto k1 = 0.0;
    auto q = 0.0;
    for (auto iter = sorted.cbegin(); iter != sorted.cend(); iter++) {
      auto w = iter->weight();
      auto dq = w / total;
      auto k2 = integratedLocation(q + dq);
      if (k2 - k1 > 1 && w != 1) {
        badWeight++;
      }
      if (k2 - k1 > 1.5 && w != 1) {
        badWeight++;
      }
      q += dq;
      k1 = k2;
    }

    return badWeight;
  }

  /**
   * Converts a quantile into a centroid scale value.  The centroid scale is nomin_ally
   * the number k of the centroid that a quantile point q should belong to.  Due to
   * round-offs, however, we can't align things perfectly without splitting points
   * and sorted.  We don't want to do that, so we have to allow for offsets.
   * In the end, the criterion is that any quantile range that spans a centroid
   * scale range more than one should be split across more than one centroid if
   * possible.  This won't be possible if the quantile range refers to a single point
   * or an already existing centroid.
   * <p/>
   * This mapping is steep near q=0 or q=1 so each centroid there will correspond to
   * less q range.  Near q=0.5, the mapping is flatter so that sorted there will
   * represent a larger chunk of quantiles.
   *
   * @param q The quantile scale value to be mapped.
   * @return The centroid scale value corresponding to q.
   */
  inline Value integratedLocation(Value q) const {
    return compression_ * (std::asin(2.0 * q - 1.0) + pi / 2) / pi;
  }

  inline Value integratedQ(Value k) const {
    return (std::sin(std::min(k, compression_) * pi / compression_ - pi / 2) + 1) / 2;
  }

  /**
   * Same as {@link #weightedAverageSorted(Value, Value, Value, Value)} but flips
   * the order of the variables if <code>x2</code> is greater than
   * <code>x1</code>.
   */
  static Value weightedAverage(Value x1, Value w1, Value x2, Value w2) {
    return (x1 <= x2) ? weightedAverageSorted(x1, w1, x2, w2) : weightedAverageSorted(x2, w2, x1, w1);
  }

  /**
   * Compute the weighted average between <code>x1</code> with a weight of
   * <code>w1</code> and <code>x2</code> with a weight of <code>w2</code>.
   * This expects <code>x1</code> to be less than or equal to <code>x2</code>
   * and is guaranteed to return a number between <code>x1</code> and
   * <code>x2</code>.
   */
  static Value weightedAverageSorted(Value x1, Value w1, Value x2, Value w2) {
    const Value x = (x1 * w1 + x2 * w2) / (w1 + w2);
    return std::max(x1, std::min(x, x2));
  }

  static Value interpolate(Value x, Value x0, Value x1) {
    return (x - x0) / (x1 - x0);
  }

  /**
   * Computes an interpolated value of a quantile that is between two sorted.
   *
   * Index is the quantile desired multiplied by the total number of samples - 1.
   *
   * @param index              Denormalized quantile desired
   * @param previousIndex      The denormalized quantile corresponding to the center of the previous centroid.
   * @param nextIndex          The denormalized quantile corresponding to the center of the following centroid.
   * @param previousMean       The mean of the previous centroid.
   * @param nextMean           The mean of the following centroid.
   * @return  The interpolated mean.
   */
  static Value quantile(Value index, Value previousIndex, Value nextIndex, Value previousMean, Value nextMean) {
    const auto delta = nextIndex - previousIndex;
    const auto previousWeight = (nextIndex - index) / delta;
    const auto nextWeight = (index - previousIndex) / delta;
    return previousMean * previousWeight + nextMean * nextWeight;
  }
};

} // namespace tdigest

Coverage Report

Created: 2026-06-30 06:38

Line	Count	Source
1		/*
2		* Licensed to Derrick R. Burns under one or more
3		* contributor license agreements. See the NOTICES file distributed with
4		* this work for additional information regarding copyright ownership.
5		* The ASF licenses this file to You under the Apache License, Version 2.0
6		* (the "License"); you may not use this file except in compliance with
7		* the License. You may obtain a copy of the License at
8		*
9		* http://www.apache.org/licenses/LICENSE-2.0
10		*
11		* Unless required by applicable law or agreed to in writing, software
12		* distributed under the License is distributed on an "AS IS" BASIS,
13		* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14		* See the License for the specific language governing permissions and
15		* limitations under the License.
16		*/
17
18		#pragma once
19
20		#include <algorithm>
21		#include <cfloat>
22		#include <cmath>
23		#include <queue>
24		#include <utility>
25		#include <vector>
26
27		#ifdef min
28		#undef min
29		#endif
30
31		#ifdef max
32		#undef max
33		#endif
34
35
36		namespace duckdb_tdigest {
37
38		using Value = double;
39		using Weight = double;
40		using Index = size_t;
41
42		const size_t kHighWater = 40000;
43		const double pi = 3.14159265358979323846;
44		class Centroid {
45		public:
46	0	Centroid() : Centroid(0.0, 0.0) {
47	0	}
48
49	0	Centroid(Value mean, Weight weight) : mean_(mean), weight_(weight) {
50	0	}
51
52	0	inline Value mean() const noexcept {
53	0	return mean_;
54	0	}
55
56	0	inline Weight weight() const noexcept {
57	0	return weight_;
58	0	}
59
60	0	inline void add(const Centroid &c) {
61		// CHECK_GT(c.weight_, 0);
62	0	if (weight_ != 0.0) {
63	0	weight_ += c.weight_;
64	0	mean_ += c.weight_ * (c.mean_ - mean_) / weight_;
65	0	} else {
66	0	weight_ = c.weight_;
67	0	mean_ = c.mean_;
68	0	}
69	0	}
70
71		private:
72		Value mean_ = 0;
73		Weight weight_ = 0;
74		};
75
76		struct CentroidList {
77	0	explicit CentroidList(const std::vector<Centroid> &s) : iter(s.cbegin()), end(s.cend()) {
78	0	}
79		std::vector<Centroid>::const_iterator iter;
80		std::vector<Centroid>::const_iterator end;
81
82	0	bool advance() {
83	0	return ++iter != end;
84	0	}
85		};
86
87		class CentroidListComparator {
88		public:
89	0	CentroidListComparator() {
90	0	}
91
92	0	bool operator()(const CentroidList &left, const CentroidList &right) const {
93	0	return left.iter->mean() > right.iter->mean();
94	0	}
95		};
96
97		using CentroidListQueue = std::priority_queue<CentroidList, std::vector<CentroidList>, CentroidListComparator>;
98
99		struct CentroidComparator {
100	0	bool operator()(const Centroid &a, const Centroid &b) const {
101	0	return a.mean() < b.mean();
102	0	}
103		};
104
105		class TDigest {
106		class TDigestComparator {
107		public:
108	0	TDigestComparator() {
109	0	}
110
111	0	bool operator()(const TDigest left, const TDigest right) const {
112	0	return left->totalSize() > right->totalSize();
113	0	}
114		};
115
116		using TDigestQueue = std::priority_queue<const TDigest , std::vector<const TDigest >, TDigestComparator>;
117
118		public:
119	0	TDigest() : TDigest(1000) {
120	0	}
121
122	0	explicit TDigest(Value compression) : TDigest(compression, 0) {
123	0	}
124
125	0	TDigest(Value compression, Index bufferSize) : TDigest(compression, bufferSize, 0) {
126	0	}
127
128		TDigest(Value compression, Index unmergedSize, Index mergedSize)
129	0	: compression_(compression), maxProcessed_(processedSize(mergedSize, compression)),
130	0	maxUnprocessed_(unprocessedSize(unmergedSize, compression)) {
131	0	processed_.reserve(maxProcessed_);
132	0	unprocessed_.reserve(maxUnprocessed_ + 1);
133	0	}
134
135		TDigest(std::vector<Centroid> &&processed, std::vector<Centroid> &&unprocessed, Value compression,
136		Index unmergedSize, Index mergedSize)
137	0	: TDigest(compression, unmergedSize, mergedSize) {
138	0	processed_ = std::move(processed);
139	0	unprocessed_ = std::move(unprocessed);
140
141	0	processedWeight_ = weight(processed_);
142	0	unprocessedWeight_ = weight(unprocessed_);
143	0	if (!processed_.empty()) {
144	0	min_ = std::min(min_, processed_[0].mean());
145	0	max_ = std::max(max_, (processed_.cend() - 1)->mean());
146	0	}
147	0	updateCumulative();
148	0	}
149
150	0	static Weight weight(std::vector<Centroid> &centroids) noexcept {
151	0	Weight w = 0.0;
152	0	for (auto centroid : centroids) {
153	0	w += centroid.weight();
154	0	}
155	0	return w;
156	0	}
157
158	0	TDigest &operator=(TDigest &&o) {
159	0	compression_ = o.compression_;
160	0	maxProcessed_ = o.maxProcessed_;
161	0	maxUnprocessed_ = o.maxUnprocessed_;
162	0	processedWeight_ = o.processedWeight_;
163	0	unprocessedWeight_ = o.unprocessedWeight_;
164	0	processed_ = std::move(o.processed_);
165	0	unprocessed_ = std::move(o.unprocessed_);
166	0	cumulative_ = std::move(o.cumulative_);
167	0	min_ = o.min_;
168	0	max_ = o.max_;
169	0	return *this;
170	0	}
171
172		TDigest(TDigest &&o)
173		: TDigest(std::move(o.processed_), std::move(o.unprocessed_), o.compression_, o.maxUnprocessed_,
174	0	o.maxProcessed_) {
175	0	}
176
177	0	static inline Index processedSize(Index size, Value compression) noexcept {
178	0	return (size == 0) ? static_cast<Index>(2 * std::ceil(compression)) : size;
179	0	}
180
181	0	static inline Index unprocessedSize(Index size, Value compression) noexcept {
182	0	return (size == 0) ? static_cast<Index>(8 * std::ceil(compression)) : size;
183	0	}
184
185		// merge in another t-digest
186	0	inline void merge(const TDigest *other) {
187	0	std::vector<const TDigest *> others {other};
188	0	add(others.cbegin(), others.cend());
189	0	}
190
191	0	const std::vector<Centroid> &processed() const {
192	0	return processed_;
193	0	}
194
195	0	const std::vector<Centroid> &unprocessed() const {
196	0	return unprocessed_;
197	0	}
198
199	0	Index maxUnprocessed() const {
200	0	return maxUnprocessed_;
201	0	}
202
203	0	Index maxProcessed() const {
204	0	return maxProcessed_;
205	0	}
206
207	0	Value min() const {
208	0	return min_;
209	0	}
210
211	0	Value max() const {
212	0	return max_;
213	0	}
214
215		// override the tracked minimum/maximum - used when reconstructing a serialized digest
216	0	void setMinMax(Value min, Value max) {
217	0	min_ = min;
218	0	max_ = max;
219	0	}
220
221	0	inline void add(std::vector<const TDigest *> digests) {
222	0	add(digests.cbegin(), digests.cend());
223	0	}
224
225		// merge in a vector of tdigests in the most efficient manner possible
226		// in constant space
227		// works for any value of kHighWater
228	0	void add(std::vector<const TDigest >::const_iterator iter, std::vector<const TDigest >::const_iterator end) {
229	0	if (iter != end) {
230	0	auto size = std::distance(iter, end);
231	0	TDigestQueue pq(TDigestComparator {});
232	0	for (; iter != end; iter++) {
233	0	pq.push((*iter));
234	0	}
235	0	std::vector<const TDigest *> batch;
236	0	batch.reserve(size_t(size));
237
238	0	size_t totalSize = 0;
239	0	while (!pq.empty()) {
240	0	auto td = pq.top();
241	0	batch.push_back(td);
242	0	pq.pop();
243	0	totalSize += td->totalSize();
244	0	if (totalSize >= kHighWater \|\| pq.empty()) {
245	0	mergeProcessed(batch);
246	0	mergeUnprocessed(batch);
247	0	processIfNecessary();
248	0	batch.clear();
249	0	totalSize = 0;
250	0	}
251	0	}
252	0	updateCumulative();
253	0	}
254	0	}
255
256	0	Weight processedWeight() const {
257	0	return processedWeight_;
258	0	}
259
260	0	Weight unprocessedWeight() const {
261	0	return unprocessedWeight_;
262	0	}
263
264	0	bool haveUnprocessed() const {
265	0	return unprocessed_.size() > 0;
266	0	}
267
268	0	size_t totalSize() const {
269	0	return processed_.size() + unprocessed_.size();
270	0	}
271
272	0	long totalWeight() const {
273	0	return static_cast<long>(processedWeight_ + unprocessedWeight_);
274	0	}
275
276		// return the cdf on the t-digest
277	0	Value cdf(Value x) {
278	0	if (haveUnprocessed() \|\| isDirty()) {
279	0	process();
280	0	}
281	0	return cdfProcessed(x);
282	0	}
283
284	0	bool isDirty() {
285	0	return processed_.size() > maxProcessed_ \|\| unprocessed_.size() > maxUnprocessed_;
286	0	}
287
288		// return the cdf on the processed values
289	0	Value cdfProcessed(Value x) const {
290	0	if (processed_.empty()) {
291	0	// no data to examin_e
292	0
293	0	return 0.0;
294	0	} else if (processed_.size() == 1) {
295	0	// exactly one centroid, should have max_==min_
296	0	auto width = max_ - min_;
297	0	if (x < min_) {
298	0	return 0.0;
299	0	} else if (x > max_) {
300	0	return 1.0;
301	0	} else if (x - min_ <= width) {
302	0	// min_ and max_ are too close together to do any viable interpolation
303	0	return 0.5;
304	0	} else {
305	0	// interpolate if somehow we have weight > 0 and max_ != min_
306	0	return (x - min_) / (max_ - min_);
307	0	}
308	0	} else {
309	0	auto n = processed_.size();
310	0	if (x <= min_) {
311	0	return 0;
312	0	}
313	0
314	0	if (x >= max_) {
315	0	return 1;
316	0	}
317	0
318	0	// check for the left tail
319	0	if (x <= mean(0)) {
320	0
321	0	// note that this is different than mean(0) > min_ ... this guarantees interpolation works
322	0	if (mean(0) - min_ > 0) {
323	0	return (x - min_) / (mean(0) - min_) * weight(0) / processedWeight_ / 2.0;
324	0	} else {
325	0	return 0;
326	0	}
327	0	}
328	0
329	0	// and the right tail
330	0	if (x >= mean(n - 1)) {
331	0	if (max_ - mean(n - 1) > 0) {
332	0	return 1.0 - (max_ - x) / (max_ - mean(n - 1)) * weight(n - 1) / processedWeight_ / 2.0;
333	0	} else {
334	0	return 1;
335	0	}
336	0	}
337	0
338	0	CentroidComparator cc;
339	0	auto iter = std::upper_bound(processed_.cbegin(), processed_.cend(), Centroid(x, 0), cc);
340	0
341	0	auto i = size_t(std::distance(processed_.cbegin(), iter));
342	0	auto z1 = x - (iter - 1)->mean();
343	0	auto z2 = (iter)->mean() - x;
344	0	return weightedAverage(cumulative_[i - 1], z2, cumulative_[i], z1) / processedWeight_;
345	0	}
346	0	}
347
348		// this returns a quantile on the t-digest
349	0	Value quantile(Value q) {
350	0	if (haveUnprocessed() \|\| isDirty()) {
351	0	process();
352	0	}
353	0	return quantileProcessed(q);
354	0	}
355
356		// this returns a quantile on the currently processed values without changing the t-digest
357		// the value will not represent the unprocessed values
358	0	Value quantileProcessed(Value q) const {
359	0	if (q < 0 \|\| q > 1) {
360	0	return NAN;
361	0	}
362
363	0	if (processed_.size() == 0) {
364		// no sorted means no data, no way to get a quantile
365	0	return NAN;
366	0	} else if (processed_.size() == 1) {
367		// with one data point, all quantiles lead to Rome
368
369	0	return mean(0);
370	0	}
371
372		// we know that there are at least two sorted now
373	0	auto n = processed_.size();
374
375		// if values were stored in a sorted array, index would be the offset we are Weighterested in
376	0	const auto index = q * processedWeight_;
377
378		// at the boundaries, we return min_ or max_
379	0	if (index <= weight(0) / 2.0) {
380	0	return min_ + 2.0 * index / weight(0) * (mean(0) - min_);
381	0	}
382
383	0	auto iter = std::lower_bound(cumulative_.cbegin(), cumulative_.cend(), index);
384
385	0	if (iter + 1 != cumulative_.cend()) {
386	0	auto i = size_t(std::distance(cumulative_.cbegin(), iter));
387	0	auto z1 = index - *(iter - 1);
388	0	auto z2 = *(iter)-index;
389		// LOG(INFO) << "z2 " << z2 << " index " << index << " z1 " << z1;
390	0	return weightedAverage(mean(i - 1), z2, mean(i), z1);
391	0	}
392	0	auto z1 = index - processedWeight_ - weight(n - 1) / 2.0;
393	0	auto z2 = weight(n - 1) / 2 - z1;
394	0	return weightedAverage(mean(n - 1), z1, max_, z2);
395	0	}
396
397	0	Value compression() const {
398	0	return compression_;
399	0	}
400
401	0	void add(Value x) {
402	0	add(x, 1);
403	0	}
404
405	0	inline void compress() {
406	0	process();
407	0	}
408
409		// add a single centroid to the unprocessed vector, processing previously unprocessed sorted if our limit has
410		// been reached.
411	0	inline bool add(Value x, Weight w) {
412	0	if (std::isnan(x)) {
413	0	return false;
414	0	}
415	0	unprocessed_.push_back(Centroid(x, w));
416	0	unprocessedWeight_ += w;
417	0	processIfNecessary();
418	0	return true;
419	0	}
420
421	0	inline void add(std::vector<Centroid>::const_iterator iter, std::vector<Centroid>::const_iterator end) {
422	0	while (iter != end) {
423	0	const size_t diff = size_t(std::distance(iter, end));
424	0	const size_t room = maxUnprocessed_ - unprocessed_.size();
425	0	auto mid = iter + int64_t(std::min(diff, room));
426	0	while (iter != mid) {
427	0	unprocessed_.push_back(*(iter++));
428	0	}
429	0	if (unprocessed_.size() >= maxUnprocessed_) {
430	0	process();
431	0	}
432	0	}
433	0	}
434
435		private:
436		Value compression_;
437
438		Value min_ = std::numeric_limits<Value>::max();
439
440		Value max_ = std::numeric_limits<Value>::min();
441
442		Index maxProcessed_;
443
444		Index maxUnprocessed_;
445
446		Value processedWeight_ = 0.0;
447
448		Value unprocessedWeight_ = 0.0;
449
450		std::vector<Centroid> processed_;
451
452		std::vector<Centroid> unprocessed_;
453
454		std::vector<Weight> cumulative_;
455
456		// return mean of i-th centroid
457	0	inline Value mean(size_t i) const noexcept {
458	0	return processed_[i].mean();
459	0	}
460
461		// return weight of i-th centroid
462	0	inline Weight weight(size_t i) const noexcept {
463	0	return processed_[i].weight();
464	0	}
465
466		// append all unprocessed centroids into current unprocessed vector
467	0	void mergeUnprocessed(const std::vector<const TDigest *> &tdigests) {
468	0	if (tdigests.size() == 0) {
469	0	return;
470	0	}
471
472	0	size_t total = unprocessed_.size();
473	0	for (auto &td : tdigests) {
474	0	total += td->unprocessed_.size();
475	0	}
476
477	0	unprocessed_.reserve(total);
478	0	for (auto &td : tdigests) {
479	0	unprocessed_.insert(unprocessed_.end(), td->unprocessed_.cbegin(), td->unprocessed_.cend());
480	0	unprocessedWeight_ += td->unprocessedWeight_;
481	0	}
482	0	}
483
484		// merge all processed centroids together into a single sorted vector
485	0	void mergeProcessed(const std::vector<const TDigest *> &tdigests) {
486	0	if (tdigests.size() == 0) {
487	0	return;
488	0	}
489
490	0	size_t total = 0;
491	0	CentroidListQueue pq(CentroidListComparator {});
492	0	for (auto &td : tdigests) {
493	0	auto &sorted = td->processed_;
494	0	auto size = sorted.size();
495	0	if (size > 0) {
496	0	pq.push(CentroidList(sorted));
497	0	total += size;
498	0	processedWeight_ += td->processedWeight_;
499	0	}
500	0	}
501	0	if (total == 0) {
502	0	return;
503	0	}
504
505	0	if (processed_.size() > 0) {
506	0	pq.push(CentroidList(processed_));
507	0	total += processed_.size();
508	0	}
509
510	0	std::vector<Centroid> sorted;
511	0	sorted.reserve(total);
512
513	0	while (!pq.empty()) {
514	0	auto best = pq.top();
515	0	pq.pop();
516	0	sorted.push_back(*(best.iter));
517	0	if (best.advance()) {
518	0	pq.push(best);
519	0	}
520	0	}
521	0	processed_ = std::move(sorted);
522	0	if (processed_.size() > 0) {
523	0	min_ = std::min(min_, processed_[0].mean());
524	0	max_ = std::max(max_, (processed_.cend() - 1)->mean());
525	0	}
526	0	}
527
528	0	inline void processIfNecessary() {
529	0	if (isDirty()) {
530	0	process();
531	0	}
532	0	}
533
534	0	void updateCumulative() {
535	0	const auto n = processed_.size();
536	0	cumulative_.clear();
537	0	cumulative_.reserve(n + 1);
538	0	auto previous = 0.0;
539	0	for (Index i = 0; i < n; i++) {
540	0	auto current = weight(i);
541	0	auto halfCurrent = current / 2.0;
542	0	cumulative_.push_back(previous + halfCurrent);
543	0	previous = previous + current;
544	0	}
545	0	cumulative_.push_back(previous);
546	0	}
547
548		// merges unprocessed_ centroids and processed_ centroids together and processes them
549		// when complete, unprocessed_ will be empty and processed_ will have at most maxProcessed_ centroids
550	0	inline void process() {
551	0	CentroidComparator cc;
552	0	std::sort(unprocessed_.begin(), unprocessed_.end(), cc);
553	0	auto count = unprocessed_.size();
554	0	unprocessed_.insert(unprocessed_.end(), processed_.cbegin(), processed_.cend());
555	0	std::inplace_merge(unprocessed_.begin(), unprocessed_.begin() + int64_t(count), unprocessed_.end(), cc);
556
557	0	processedWeight_ += unprocessedWeight_;
558	0	unprocessedWeight_ = 0;
559	0	processed_.clear();
560
561	0	processed_.push_back(unprocessed_[0]);
562	0	Weight wSoFar = unprocessed_[0].weight();
563	0	Weight wLimit = processedWeight_ * integratedQ(1.0);
564
565	0	auto end = unprocessed_.end();
566	0	for (auto iter = unprocessed_.cbegin() + 1; iter < end; iter++) {
567	0	auto &centroid = *iter;
568	0	Weight projectedW = wSoFar + centroid.weight();
569	0	if (projectedW <= wLimit) {
570	0	wSoFar = projectedW;
571	0	(processed_.end() - 1)->add(centroid);
572	0	} else {
573	0	auto k1 = integratedLocation(wSoFar / processedWeight_);
574	0	wLimit = processedWeight_ * integratedQ(k1 + 1.0);
575	0	wSoFar += centroid.weight();
576	0	processed_.emplace_back(centroid);
577	0	}
578	0	}
579	0	unprocessed_.clear();
580	0	min_ = std::min(min_, processed_[0].mean());
581	0	max_ = std::max(max_, (processed_.cend() - 1)->mean());
582	0	updateCumulative();
583	0	}
584
585	0	inline size_t checkWeights() {
586	0	return checkWeights(processed_, processedWeight_);
587	0	}
588
589	0	size_t checkWeights(const std::vector<Centroid> &sorted, Value total) {
590	0	size_t badWeight = 0;
591	0	auto k1 = 0.0;
592	0	auto q = 0.0;
593	0	for (auto iter = sorted.cbegin(); iter != sorted.cend(); iter++) {
594	0	auto w = iter->weight();
595	0	auto dq = w / total;
596	0	auto k2 = integratedLocation(q + dq);
597	0	if (k2 - k1 > 1 && w != 1) {
598	0	badWeight++;
599	0	}
600	0	if (k2 - k1 > 1.5 && w != 1) {
601	0	badWeight++;
602	0	}
603	0	q += dq;
604	0	k1 = k2;
605	0	}
606	0
607	0	return badWeight;
608	0	}
609
610		/**
611		* Converts a quantile into a centroid scale value. The centroid scale is nomin_ally
612		* the number k of the centroid that a quantile point q should belong to. Due to
613		* round-offs, however, we can't align things perfectly without splitting points
614		* and sorted. We don't want to do that, so we have to allow for offsets.
615		* In the end, the criterion is that any quantile range that spans a centroid
616		* scale range more than one should be split across more than one centroid if
617		* possible. This won't be possible if the quantile range refers to a single point
618		* or an already existing centroid.
619		* <p/>
620		* This mapping is steep near q=0 or q=1 so each centroid there will correspond to
621		* less q range. Near q=0.5, the mapping is flatter so that sorted there will
622		* represent a larger chunk of quantiles.
623		*
624		* @param q The quantile scale value to be mapped.
625		* @return The centroid scale value corresponding to q.
626		*/
627	0	inline Value integratedLocation(Value q) const {
628	0	return compression_ * (std::asin(2.0 * q - 1.0) + pi / 2) / pi;
629	0	}
630
631	0	inline Value integratedQ(Value k) const {
632	0	return (std::sin(std::min(k, compression_) * pi / compression_ - pi / 2) + 1) / 2;
633	0	}
634
635		/**
636		* Same as {@link #weightedAverageSorted(Value, Value, Value, Value)} but flips
637		* the order of the variables if <code>x2</code> is greater than
638		* <code>x1</code>.
639		*/
640	0	static Value weightedAverage(Value x1, Value w1, Value x2, Value w2) {
641	0	return (x1 <= x2) ? weightedAverageSorted(x1, w1, x2, w2) : weightedAverageSorted(x2, w2, x1, w1);
642	0	}
643
644		/**
645		* Compute the weighted average between <code>x1</code> with a weight of
646		* <code>w1</code> and <code>x2</code> with a weight of <code>w2</code>.
647		* This expects <code>x1</code> to be less than or equal to <code>x2</code>
648		* and is guaranteed to return a number between <code>x1</code> and
649		* <code>x2</code>.
650		*/
651	0	static Value weightedAverageSorted(Value x1, Value w1, Value x2, Value w2) {
652	0	const Value x = (x1 * w1 + x2 * w2) / (w1 + w2);
653	0	return std::max(x1, std::min(x, x2));
654	0	}
655
656	0	static Value interpolate(Value x, Value x0, Value x1) {
657	0	return (x - x0) / (x1 - x0);
658	0	}
659
660		/**
661		* Computes an interpolated value of a quantile that is between two sorted.
662		*
663		* Index is the quantile desired multiplied by the total number of samples - 1.
664		*
665		* @param index Denormalized quantile desired
666		* @param previousIndex The denormalized quantile corresponding to the center of the previous centroid.
667		* @param nextIndex The denormalized quantile corresponding to the center of the following centroid.
668		* @param previousMean The mean of the previous centroid.
669		* @param nextMean The mean of the following centroid.
670		* @return The interpolated mean.
671		*/
672	0	static Value quantile(Value index, Value previousIndex, Value nextIndex, Value previousMean, Value nextMean) {
673	0	const auto delta = nextIndex - previousIndex;
674	0	const auto previousWeight = (nextIndex - index) / delta;
675	0	const auto nextWeight = (index - previousIndex) / delta;
676	0	return previousMean * previousWeight + nextMean * nextWeight;
677	0	}
678		};
679
680		} // namespace tdigest