#pragma once #include #include #include #include #include #include #include #include #include #include #include #include #include namespace DB { namespace ErrorCodes { extern const int LOGICAL_ERROR; } } /// Implementation of Reservoir Sampling algorithm. Incrementally selects from the added objects a random subset of the `sample_count` size. /// Can approximately get quantiles. /// The `quantile` call takes O(sample_count log sample_count), if after the previous call `quantile` there was at least one call to insert. Otherwise, O(1). /// That is, it makes sense to first add, then get quantiles without adding. namespace DB { namespace ErrorCodes { extern const int MEMORY_LIMIT_EXCEEDED; } } namespace detail { const size_t DEFAULT_SAMPLE_COUNT = 8192; const auto MAX_SKIP_DEGREE = sizeof(UInt32) * 8; } /// What if there is not a single value - throw an exception, or return 0 or NaN in the case of double? enum class ReservoirSamplerDeterministicOnEmpty { THROW, RETURN_NAN_OR_ZERO, }; template class ReservoirSamplerDeterministic { bool good(const UInt32 hash) { return hash == ((hash >> skip_degree) << skip_degree); } public: ReservoirSamplerDeterministic(const size_t sample_count_ = DEFAULT_SAMPLE_COUNT) : sample_count{sample_count_} { } void clear() { samples.clear(); sorted = false; total_values = 0; } void insert(const T & v, const UInt64 determinator) { if (isNaN(v)) return; const UInt32 hash = intHash64(determinator); if (!good(hash)) return; insertImpl(v, hash); sorted = false; ++total_values; } size_t size() const { return total_values; } T quantileNearest(double level) { if (samples.empty()) return onEmpty(); sortIfNeeded(); double index = level * (samples.size() - 1); size_t int_index = static_cast(index + 0.5); int_index = std::max(0LU, std::min(samples.size() - 1, int_index)); return samples[int_index].first; } /** If T is not a numeric type, using this method causes a compilation error, * but use of error class does not cause. SFINAE. * Not SFINAE. Functions members of type templates are simply not checked until they are used. */ double quantileInterpolated(double level) { if (samples.empty()) return onEmpty(); sortIfNeeded(); const double index = std::max(0., std::min(samples.size() - 1., level * (samples.size() - 1))); /// To get a value from a fractional index, we linearly interpolate between adjacent values. size_t left_index = static_cast(index); size_t right_index = left_index + 1; if (right_index == samples.size()) return static_cast(samples[left_index].first); const double left_coef = right_index - index; const double right_coef = index - left_index; return static_cast(samples[left_index].first) * left_coef + static_cast(samples[right_index].first) * right_coef; } void merge(const ReservoirSamplerDeterministic & b) { if (sample_count != b.sample_count) throw Poco::Exception("Cannot merge ReservoirSamplerDeterministic's with different sample_count"); sorted = false; if (b.skip_degree > skip_degree) { skip_degree = b.skip_degree; thinOut(); } for (const auto & sample : b.samples) if (good(sample.second)) insertImpl(sample.first, sample.second); total_values += b.total_values; } void read(DB::ReadBuffer & buf) { DB::readIntBinary(sample_count, buf); DB::readIntBinary(total_values, buf); samples.resize(std::min(total_values, sample_count)); for (size_t i = 0; i < samples.size(); ++i) DB::readPODBinary(samples[i], buf); sorted = false; } void write(DB::WriteBuffer & buf) const { DB::writeIntBinary(sample_count, buf); DB::writeIntBinary(total_values, buf); for (size_t i = 0; i < std::min(sample_count, total_values); ++i) DB::writePODBinary(samples[i], buf); } private: /// We allocate some memory on the stack to avoid allocations when there are many objects with a small number of elements. using Element = std::pair; using Array = DB::PODArray; size_t sample_count; size_t total_values{}; bool sorted{}; Array samples; UInt8 skip_degree{}; void insertImpl(const T & v, const UInt32 hash) { /// @todo why + 1? I don't quite recall while (samples.size() + 1 >= sample_count) { if (++skip_degree > detail::MAX_SKIP_DEGREE) throw DB::Exception{"skip_degree exceeds maximum value", DB::ErrorCodes::MEMORY_LIMIT_EXCEEDED}; thinOut(); } samples.emplace_back(v, hash); } void thinOut() { auto size = samples.size(); for (size_t i = 0; i < size;) { if (!good(samples[i].second)) { /// swap current element with the last one std::swap(samples[size - 1], samples[i]); --size; } else ++i; } if (size != samples.size()) { samples.resize(size); sorted = false; } } void sortIfNeeded() { if (sorted) return; sorted = true; std::sort(samples.begin(), samples.end(), [] (const std::pair & lhs, const std::pair & rhs) { return lhs.first < rhs.first; }); } template ResultType onEmpty() const { if (OnEmpty == ReservoirSamplerDeterministicOnEmpty::THROW) throw DB::Exception(DB::ErrorCodes::LOGICAL_ERROR, "Quantile of empty ReservoirSamplerDeterministic"); else return NanLikeValueConstructor>::getValue(); } };