estimated_histogram.hh

/*
 * This file is open source software, licensed to you under the terms
 * of the Apache License, Version 2.0 (the "License").  See the NOTICE file
 * distributed with this work for additional information regarding copyright
 * ownership.  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.
 */
/*
 * Copyright (C) 2023 ScyllaDB.
 */
 
#pragma once
 
#include <cmath>
#include <algorithm>
#include <vector>
#include <chrono>
#include <seastar/core/metrics_types.hh>
#include <seastar/core/print.hh>
#include <seastar/core/bitops.hh>
#include <limits>
#include <array>
 
namespace seastar::metrics::internal {
template<uint64_t Min, uint64_t Max, size_t Precision>
class approximate_exponential_histogram {
    static constexpr unsigned NUM_EXP_RANGES = log2floor(Max/Min);
    static constexpr size_t NUM_BUCKETS = NUM_EXP_RANGES * Precision + 1;
    static constexpr unsigned PRECISION_BITS = log2floor(Precision);
    static constexpr unsigned BASESHIFT = log2floor(Min);
    static constexpr uint64_t LOWER_BITS_MASK = Precision - 1;
    static constexpr size_t MIN_ID = log2ceil(Min) * Precision + 1; // id0 (x,1]
    std::array<uint64_t, NUM_BUCKETS> _buckets;
public:
    approximate_exponential_histogram() {
        clear();
    }
 
    uint64_t get_bucket_lower_limit(uint16_t bucket_id) const {
        if (bucket_id == NUM_BUCKETS - 1) {
            return Max;
        }
        int16_t exp_rang = (bucket_id >> PRECISION_BITS);
        return (Min << exp_rang) +  ((bucket_id & LOWER_BITS_MASK) << (exp_rang + BASESHIFT - PRECISION_BITS));
    }
 
    uint64_t get_bucket_upper_limit(uint16_t bucket_id) const {
        if (bucket_id == NUM_BUCKETS - 1) {
            return std::numeric_limits<uint64_t>::max();
        }
        return get_bucket_lower_limit(bucket_id + 1);
    }
 
    uint16_t find_bucket_index(uint64_t val) const {
        if (val >= Max) {
            return NUM_BUCKETS - 1;
        }
        if (val <= Min) {
            return 0;
        }
        uint16_t range = log2floor(val);
        val >>= range - PRECISION_BITS; // leave the top most N+1 bits where N is the resolution.
        return ((range - BASESHIFT) << PRECISION_BITS) + (val & LOWER_BITS_MASK);
    }
 
    void clear() {
        std::fill(_buckets.begin(), _buckets.end(), 0);
    }
 
    void add(uint64_t n) {
        _buckets.at(find_bucket_index(n))++;
    }
 
    uint64_t min() const {
        for (size_t i = 0; i < NUM_BUCKETS; i ++) {
            if (_buckets[i] > 0) {
                return get_bucket_lower_limit(i);
            }
        }
        return 0;
    }
 
    uint64_t max() const {
        for (int i = NUM_BUCKETS - 1; i >= 0; i--) {
            if (_buckets[i] > 0) {
                return get_bucket_upper_limit(i);
            }
        }
        return 0;
    }
 
    approximate_exponential_histogram& merge(const approximate_exponential_histogram& b) {
        for (size_t i = 0; i < NUM_BUCKETS; i++) {
            _buckets[i] += b.get(i);
        }
        return *this;
    }
 
    template<uint64_t A, uint64_t B, size_t C>
    friend approximate_exponential_histogram<A, B, C> merge(approximate_exponential_histogram<A, B, C> a, const approximate_exponential_histogram<A, B, C>& b);
 
    /*
     * \brief returns the count in the given bucket
     */
    uint64_t get(size_t bucket) const {
        return _buckets[bucket];
    }
 
    uint64_t quantile(float quantile) const {
        if (quantile < 0 || quantile > 1.0) {
            throw std::runtime_error("Invalid quantile value " + std::to_string(quantile) + ". Value should be between 0 and 1");
        }
        auto c = count();
 
        if (!c) {
            return 0; // no data
        }
 
        auto pcount = uint64_t(std::floor(c * quantile));
        uint64_t elements = 0;
        for (size_t i = 0; i < NUM_BUCKETS - 2; i++) {
            if (_buckets[i]) {
                elements += _buckets[i];
                if (elements >= pcount) {
                    return get_bucket_lower_limit(i);
                }
            }
        }
        return Max; // overflowed value is in the requested quantile
    }
 
    uint64_t mean() const {
        uint64_t elements = 0;
        double sum = 0;
        for (size_t i = 0; i < NUM_BUCKETS - 1; i++) {
            elements += _buckets[i];
            sum += _buckets[i] * get_bucket_lower_limit(i);
        }
        return (elements) ?  sum / elements : 0;
    }
 
    size_t size() const {
        return NUM_BUCKETS;
    }
 
    uint64_t count() const {
        uint64_t sum = 0L;
        for (size_t i = 0; i < NUM_BUCKETS; i++) {
            sum += _buckets[i];
        }
        return sum;
    }
 
    approximate_exponential_histogram& operator*=(double v) {
        for (size_t i = 0; i < NUM_BUCKETS; i++) {
            _buckets[i] *= v;
        }
        return *this;
    }
 
    uint64_t& operator[](size_t b) noexcept {
        return _buckets[b];
    }
    int32_t get_schema() const noexcept {
        return PRECISION_BITS;
    }
 
    int32_t min_as_native_histogram_id() const noexcept {
        return MIN_ID;
    }
 
    seastar::metrics::histogram to_metrics_histogram() const noexcept {
        seastar::metrics::histogram res;
        res.buckets.resize(size() - 1);
        uint64_t cummulative_count = 0;
        res.sample_sum = 0;
        res.native_histogram = seastar::metrics::native_histogram_info{get_schema(), min_as_native_histogram_id()};
 
        for (size_t i = 0; i < NUM_BUCKETS - 1; i++) {
            auto& v = res.buckets[i];
            v.upper_bound = get_bucket_lower_limit(i + 1);
            cummulative_count += get(i);
            v.count = cummulative_count;
            res.sample_sum += get(i) * v.upper_bound;
        }
        // The count serves as the infinite bucket
        res.sample_count = cummulative_count + get(NUM_BUCKETS - 1);
        res.sample_sum += get(NUM_BUCKETS - 1) * get_bucket_lower_limit(NUM_BUCKETS - 1);
        return res;
    }
};
 
template<uint64_t Min, uint64_t Max, size_t NumBuckets>
inline approximate_exponential_histogram<Min, Max, NumBuckets> base_estimated_histogram_merge(approximate_exponential_histogram<Min, Max, NumBuckets> a, const approximate_exponential_histogram<Min, Max, NumBuckets>& b) {
    return a.merge(b);
}
 
class time_estimated_histogram : public approximate_exponential_histogram<512, 33554432, 4> {
public:
    time_estimated_histogram& merge(const time_estimated_histogram& b) {
        approximate_exponential_histogram<512, 33554432, 4>::merge(b);
        return *this;
    }
 
    void add_micro(uint64_t n) {
        approximate_exponential_histogram<512, 33554432, 4>::add(n);
    }
 
    template<typename T>
    void add(const T& latency) {
        add_micro(std::chrono::duration_cast<std::chrono::microseconds>(latency).count());
    }
};
 
inline time_estimated_histogram time_estimated_histogram_merge(time_estimated_histogram a, const time_estimated_histogram& b) {
    return a.merge(b);
}
}