T tdigest::get_quantile()

in tdigest/include/tdigest_impl.hpp [151:197]

• 42 lines of code
• 16 McCabe index (conditional complexity)

T tdigest<T, A>::get_quantile(double rank) const {
  if (is_empty()) throw std::runtime_error("operation is undefined for an empty sketch");
  if ((rank < 0.0) || (rank > 1.0)) {
    throw std::invalid_argument("Normalized rank cannot be less than 0 or greater than 1");
  }
  const_cast<tdigest*>(this)->compress(); // side effect
  if (centroids_.size() == 1) return centroids_.front().get_mean();

  // at least 2 centroids
  const double weight = rank * centroids_weight_;
  if (weight < 1) return min_;
  if (weight > centroids_weight_ - 1.0) return max_;
  const double first_weight = centroids_.front().get_weight();
  if (first_weight > 1 && weight < first_weight / 2.0) {
    return min_ + (weight - 1.0) / (first_weight / 2.0 - 1.0) * (centroids_.front().get_mean() - min_);
  }
  const double last_weight = centroids_.back().get_weight();
  if (last_weight > 1 && centroids_weight_ - weight <= last_weight / 2.0) {
    return max_ + (centroids_weight_ - weight - 1.0) / (last_weight / 2.0 - 1.0) * (max_ - centroids_.back().get_mean());
  }

  // interpolate between extremes
  double weight_so_far = first_weight / 2.0;
  for (size_t i = 0; i < centroids_.size() - 1; ++i) {
    const double dw = (centroids_[i].get_weight() + centroids_[i + 1].get_weight()) / 2.0;
    if (weight_so_far + dw > weight) {
      // the target weight is between centroids i and i+1
      double left_weight = 0;
      if (centroids_[i].get_weight() == 1) {
        if (weight - weight_so_far < 0.5) return centroids_[i].get_mean();
        left_weight = 0.5;
      }
      double right_weight = 0;
      if (centroids_[i + 1].get_weight() == 1) {
        if (weight_so_far + dw - weight <= 0.5) return centroids_[i + 1].get_mean();
        right_weight = 0.5;
      }
      const double w1 = weight - weight_so_far - left_weight;
      const double w2 = weight_so_far + dw - weight - right_weight;
      return weighted_average(centroids_[i].get_mean(), w1, centroids_[i + 1].get_mean(), w2);
    }
    weight_so_far += dw;
  }
  const double w1 = weight - centroids_weight_ - centroids_.back().get_weight() / 2.0;
  const double w2 = centroids_.back().get_weight() / 2.0 - w1;
  return weighted_average(centroids_.back().get_weight(), w1, max_, w2);
}