Coverage Report

Created: 2026-10-09 14:29

next uncovered line (L), next uncovered region (R), next uncovered branch (B)
/root/doris/be/src/util/tdigest.h
Line
Count
Source
1
// Licensed to the Apache Software Foundation (ASF) under one
2
// or more contributor license agreements.  See the NOTICE file
3
// distributed with this work for additional information
4
// regarding copyright ownership.  The ASF licenses this file
5
// to you under the Apache License, Version 2.0 (the
6
// "License"); you may not use this file except in compliance
7
// with 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,
12
// software distributed under the License is distributed on an
13
// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
14
// KIND, either express or implied.  See the License for the
15
// specific language governing permissions and limitations
16
// under the License.
17
18
/*
19
 * Licensed to Derrick R. Burns under one or more
20
 * contributor license agreements.  See the NOTICES file distributed with
21
 * this work for additional information regarding copyright ownership.
22
 * The ASF licenses this file to You under the Apache License, Version 2.0
23
 * (the "License"); you may not use this file except in compliance with
24
 * the License.  You may obtain a copy of the License at
25
 *
26
 *     http://www.apache.org/licenses/LICENSE-2.0
27
 *
28
 * Unless required by applicable law or agreed to in writing, software
29
 * distributed under the License is distributed on an "AS IS" BASIS,
30
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
31
 * See the License for the specific language governing permissions and
32
 * limitations under the License.
33
 */
34
35
// T-Digest :  Percentile and Quantile Estimation of Big Data
36
// A new data structure for accurate on-line accumulation of rank-based statistics
37
// such as quantiles and trimmed means.
38
// See original paper: "Computing extremely accurate quantiles using t-digest"
39
// by Ted Dunning and Otmar Ertl for more details
40
// https://github.com/tdunning/t-digest/blob/07b8f2ca2be8d0a9f04df2feadad5ddc1bb73c88/docs/t-digest-paper/histo.pdf.
41
// https://github.com/derrickburns/tdigest
42
43
#pragma once
44
45
#include <pdqsort.h>
46
47
#include <algorithm>
48
#include <cfloat>
49
#include <cmath>
50
#include <iostream>
51
#include <memory>
52
#include <queue>
53
#include <utility>
54
#include <vector>
55
56
#include "common/factory_creator.h"
57
#include "common/logging.h"
58
59
namespace doris {
60
#include "common/compile_check_begin.h"
61
62
using Value = float;
63
using Weight = float;
64
using Index = size_t;
65
66
const size_t kHighWater = 40000;
67
68
class Centroid {
69
public:
70
423k
    Centroid() : Centroid(0.0, 0.0) {}
71
72
652k
    Centroid(Value mean, Weight weight) : _mean(mean), _weight(weight) {}
73
74
297M
    Value mean() const noexcept { return _mean; }
75
76
439M
    Weight weight() const noexcept { return _weight; }
77
78
27.8k
    Value& mean() noexcept { return _mean; }
79
80
26.2k
    Weight& weight() noexcept { return _weight; }
81
82
1.08M
    void add(const Centroid& c) {
83
1.08M
        DCHECK_GT(c._weight, 0);
84
1.08M
        if (_weight != 0.0) {
85
1.08M
            _weight += c._weight;
86
1.08M
            _mean += c._weight * (c._mean - _mean) / _weight;
87
18.4E
        } else {
88
18.4E
            _weight = c._weight;
89
18.4E
            _mean = c._mean;
90
18.4E
        }
91
1.08M
    }
92
93
private:
94
    Value _mean = 0;
95
    Weight _weight = 0;
96
};
97
98
struct CentroidList {
99
845
    CentroidList(const std::vector<Centroid>& s) : iter(s.cbegin()), end(s.cend()) {}
100
    std::vector<Centroid>::const_iterator iter;
101
    std::vector<Centroid>::const_iterator end;
102
103
1.99M
    bool advance() { return ++iter != end; }
104
};
105
106
class CentroidListComparator {
107
public:
108
    CentroidListComparator() = default;
109
110
1.83M
    bool operator()(const CentroidList& left, const CentroidList& right) const {
111
1.83M
        return left.iter->mean() > right.iter->mean();
112
1.83M
    }
113
};
114
115
using CentroidListQueue =
116
        std::priority_queue<CentroidList, std::vector<CentroidList>, CentroidListComparator>;
117
118
struct CentroidComparator {
119
147M
    bool operator()(const Centroid& a, const Centroid& b) const { return a.mean() < b.mean(); }
120
};
121
122
class TDigest {
123
    ENABLE_FACTORY_CREATOR(TDigest);
124
125
    class TDigestComparator {
126
    public:
127
        TDigestComparator() = default;
128
129
0
        bool operator()(const TDigest* left, const TDigest* right) const {
130
0
            return left->totalSize() > right->totalSize();
131
0
        }
132
    };
133
    using TDigestQueue =
134
            std::priority_queue<const TDigest*, std::vector<const TDigest*>, TDigestComparator>;
135
136
public:
137
0
    TDigest() : TDigest(10000) {}
138
139
1.40k
    explicit TDigest(Value compression) : TDigest(compression, 0) {}
140
141
1.40k
    TDigest(Value compression, Index bufferSize) : TDigest(compression, bufferSize, 0) {}
142
143
    TDigest(Value compression, Index unmergedSize, Index mergedSize)
144
1.40k
            : _compression(compression),
145
1.40k
              _max_processed(processedSize(mergedSize, compression)),
146
1.40k
              _max_unprocessed(unprocessedSize(unmergedSize, compression)) {
147
1.40k
        _processed.reserve(_max_processed);
148
1.40k
        _unprocessed.reserve(_max_unprocessed + 1);
149
1.40k
    }
150
151
    TDigest(std::vector<Centroid>&& processed, std::vector<Centroid>&& unprocessed,
152
            Value compression, Index unmergedSize, Index mergedSize)
153
0
            : TDigest(compression, unmergedSize, mergedSize) {
154
0
        _processed = std::move(processed);
155
0
        _unprocessed = std::move(unprocessed);
156
0
157
0
        _processed_weight = weight(_processed);
158
0
        _unprocessed_weight = weight(_unprocessed);
159
0
        if (_processed.size() > 0) {
160
0
            _min = std::min(_min, _processed[0].mean());
161
0
            _max = std::max(_max, (_processed.cend() - 1)->mean());
162
0
        }
163
0
        updateCumulative();
164
0
    }
165
166
0
    static Weight weight(std::vector<Centroid>& centroids) noexcept {
167
0
        Weight w = 0.0;
168
0
        for (auto centroid : centroids) {
169
0
            w += centroid.weight();
170
0
        }
171
0
        return w;
172
0
    }
173
174
24.0k
    TDigest(const TDigest&) = default;
175
176
0
    TDigest& operator=(TDigest&& o) {
177
0
        _compression = o._compression;
178
0
        _max_processed = o._max_processed;
179
0
        _max_unprocessed = o._max_unprocessed;
180
0
        _processed_weight = o._processed_weight;
181
0
        _unprocessed_weight = o._unprocessed_weight;
182
0
        _processed = std::move(o._processed);
183
0
        _unprocessed = std::move(o._unprocessed);
184
0
        _cumulative = std::move(o._cumulative);
185
0
        _min = o._min;
186
0
        _max = o._max;
187
0
        return *this;
188
0
    }
189
190
    TDigest(TDigest&& o)
191
            : TDigest(std::move(o._processed), std::move(o._unprocessed), o._compression,
192
0
                      o._max_unprocessed, o._max_processed) {}
193
194
1.40k
    static inline Index processedSize(Index size, Value compression) noexcept {
195
1.40k
        return (size == 0) ? static_cast<Index>(2 * std::ceil(compression)) : size;
196
1.40k
    }
197
198
1.40k
    static inline Index unprocessedSize(Index size, Value compression) noexcept {
199
1.40k
        return (size == 0) ? static_cast<Index>(8 * std::ceil(compression)) : size;
200
1.40k
    }
201
202
    // merge in another t-digest
203
718
    void merge(const TDigest* other) {
204
718
        std::vector<const TDigest*> others {other};
205
718
        add(others.cbegin(), others.cend());
206
718
    }
207
208
    const std::vector<Centroid>& processed() const { return _processed; }
209
210
    const std::vector<Centroid>& unprocessed() const { return _unprocessed; }
211
212
0
    Index maxUnprocessed() const { return _max_unprocessed; }
213
214
0
    Index maxProcessed() const { return _max_processed; }
215
216
    void add(std::vector<const TDigest*> digests) { add(digests.cbegin(), digests.cend()); }
217
218
    // merge in a vector of tdigests in the most efficient manner possible
219
    // in constant space
220
    // works for any value of kHighWater
221
    void add(std::vector<const TDigest*>::const_iterator iter,
222
719
             std::vector<const TDigest*>::const_iterator end) {
223
719
        if (iter != end) {
224
719
            auto size = std::distance(iter, end);
225
719
            TDigestQueue pq(TDigestComparator {});
226
1.43k
            for (; iter != end; iter++) {
227
719
                pq.push((*iter));
228
719
            }
229
719
            std::vector<const TDigest*> batch;
230
719
            batch.reserve(size);
231
232
719
            size_t totalSize = 0;
233
1.43k
            while (!pq.empty()) {
234
719
                const auto* td = pq.top();
235
719
                batch.push_back(td);
236
719
                pq.pop();
237
719
                totalSize += td->totalSize();
238
719
                if (totalSize >= kHighWater || pq.empty()) {
239
719
                    mergeProcessed(batch);
240
719
                    mergeUnprocessed(batch);
241
719
                    processIfNecessary();
242
719
                    batch.clear();
243
719
                    totalSize = 0;
244
719
                }
245
719
            }
246
719
            updateCumulative();
247
719
        }
248
719
    }
249
250
0
    Weight processedWeight() const { return _processed_weight; }
251
252
0
    Weight unprocessedWeight() const { return _unprocessed_weight; }
253
254
49.9k
    bool haveUnprocessed() const { return _unprocessed.size() > 0; }
255
256
719
    size_t totalSize() const { return _processed.size() + _unprocessed.size(); }
257
258
    long totalWeight() const { return static_cast<long>(_processed_weight + _unprocessed_weight); }
259
260
    // return the cdf on the t-digest
261
    Value cdf(Value x) {
262
        if (haveUnprocessed() || isDirty()) {
263
            process();
264
        }
265
        return cdfProcessed(x);
266
    }
267
268
230k
    bool isDirty() {
269
230k
        return _processed.size() > _max_processed || _unprocessed.size() > _max_unprocessed;
270
230k
    }
271
272
    // return the cdf on the processed values
273
0
    Value cdfProcessed(Value x) const {
274
0
        VLOG_CRITICAL << "cdf value " << x;
275
0
        VLOG_CRITICAL << "processed size " << _processed.size();
276
0
        if (_processed.size() == 0) {
277
0
            // no data to examine
278
0
            VLOG_CRITICAL << "no processed values";
279
0
280
0
            return 0.0;
281
0
        } else if (_processed.size() == 1) {
282
0
            VLOG_CRITICAL << "one processed value "
283
0
                          << " _min " << _min << " _max " << _max;
284
0
            // exactly one centroid, should have _max==_min
285
0
            auto width = _max - _min;
286
0
            if (x < _min) {
287
0
                return 0.0;
288
0
            } else if (x > _max) {
289
0
                return 1.0;
290
0
            } else if (x - _min <= width) {
291
0
                // _min and _max are too close together to do any viable interpolation
292
0
                return 0.5;
293
0
            } else {
294
0
                // interpolate if somehow we have weight > 0 and _max != _min
295
0
                return (x - _min) / (_max - _min);
296
0
            }
297
0
        } else {
298
0
            auto n = _processed.size();
299
0
            if (x <= _min) {
300
0
                VLOG_CRITICAL << "below _min "
301
0
                              << " _min " << _min << " x " << x;
302
0
                return 0;
303
0
            }
304
0
305
0
            if (x >= _max) {
306
0
                VLOG_CRITICAL << "above _max "
307
0
                              << " _max " << _max << " x " << x;
308
0
                return 1;
309
0
            }
310
0
311
0
            // check for the left tail
312
0
            if (x <= mean(0)) {
313
0
                VLOG_CRITICAL << "left tail "
314
0
                              << " _min " << _min << " mean(0) " << mean(0) << " x " << x;
315
0
316
0
                // note that this is different than mean(0) > _min ... this guarantees interpolation works
317
0
                if (mean(0) - _min > 0) {
318
0
                    return static_cast<Value>((x - _min) / (mean(0) - _min) * weight(0) /
319
0
                                              _processed_weight / 2.0);
320
0
                } else {
321
0
                    return 0;
322
0
                }
323
0
            }
324
0
325
0
            // and the right tail
326
0
            if (x >= mean(n - 1)) {
327
0
                VLOG_CRITICAL << "right tail"
328
0
                              << " _max " << _max << " mean(n - 1) " << mean(n - 1) << " x " << x;
329
0
330
0
                if (_max - mean(n - 1) > 0) {
331
0
                    return static_cast<Value>(1.0 - (_max - x) / (_max - mean(n - 1)) *
332
0
                                                            weight(n - 1) / _processed_weight /
333
0
                                                            2.0);
334
0
                } else {
335
0
                    return 1;
336
0
                }
337
0
            }
338
0
339
0
            CentroidComparator cc;
340
0
            auto iter =
341
0
                    std::upper_bound(_processed.cbegin(), _processed.cend(), Centroid(x, 0), cc);
342
0
343
0
            auto i = std::distance(_processed.cbegin(), iter);
344
0
            auto z1 = x - (iter - 1)->mean();
345
0
            auto z2 = (iter)->mean() - x;
346
0
            DCHECK_LE(0.0, z1);
347
0
            DCHECK_LE(0.0, z2);
348
0
            VLOG_CRITICAL << "middle "
349
0
                          << " z1 " << z1 << " z2 " << z2 << " x " << x;
350
0
351
0
            return weightedAverage(_cumulative[i - 1], z2, _cumulative[i], z1) / _processed_weight;
352
0
        }
353
0
    }
354
355
    // this returns a quantile on the t-digest
356
524
    Value quantile(Value q) {
357
524
        if (haveUnprocessed() || isDirty()) {
358
443
            process();
359
443
        }
360
524
        return quantileProcessed(q);
361
524
    }
362
363
    // this returns a quantile on the currently processed values without changing the t-digest
364
    // the value will not represent the unprocessed values
365
1.43k
    Value quantileProcessed(Value q) const {
366
1.43k
        if (q < 0 || q > 1) {
367
0
            VLOG_CRITICAL << "q should be in [0,1], got " << q;
368
0
            return NAN;
369
0
        }
370
371
1.43k
        if (_processed.size() == 0) {
372
            // no sorted means no data, no way to get a quantile
373
0
            return NAN;
374
1.43k
        } else if (_processed.size() == 1) {
375
            // with one data point, all quantiles lead to Rome
376
377
178
            return mean(0);
378
178
        }
379
380
        // we know that there are at least two sorted now
381
1.25k
        auto n = _processed.size();
382
383
        // if values were stored in a sorted array, index would be the offset we are Weighterested in
384
1.25k
        const auto index = q * _processed_weight;
385
386
        // at the boundaries, we return _min or _max
387
1.25k
        if (index <= weight(0) / 2.0) {
388
47
            DCHECK_GT(weight(0), 0);
389
47
            return static_cast<Value>(_min + 2.0 * index / weight(0) * (mean(0) - _min));
390
47
        }
391
392
1.20k
        auto iter = std::lower_bound(_cumulative.cbegin(), _cumulative.cend(), index);
393
394
1.20k
        if (iter != _cumulative.cend() && iter != _cumulative.cbegin() &&
395
1.20k
            iter + 1 != _cumulative.cend()) {
396
501
            auto i = std::distance(_cumulative.cbegin(), iter);
397
501
            auto z1 = index - *(iter - 1);
398
501
            auto z2 = *(iter)-index;
399
            // VLOG_CRITICAL << "z2 " << z2 << " index " << index << " z1 " << z1;
400
501
            return weightedAverage(mean(i - 1), z2, mean(i), z1);
401
501
        }
402
403
1.20k
        DCHECK_LE(index, _processed_weight);
404
704
        DCHECK_GE(index, _processed_weight - weight(n - 1) / 2.0);
405
406
704
        auto z1 = static_cast<Value>(index - _processed_weight - weight(n - 1) / 2.0);
407
704
        auto z2 = static_cast<Value>(weight(n - 1) / 2 - z1);
408
704
        return weightedAverage(mean(n - 1), z1, _max, z2);
409
1.20k
    }
410
411
0
    Value compression() const { return _compression; }
412
413
219k
    void add(Value x) { add(x, 1); }
414
415
24.1k
    void compress() { process(); }
416
417
    // add a single centroid to the unprocessed vector, processing previously unprocessed sorted if our limit has
418
    // been reached.
419
229k
    bool add(Value x, Weight w) {
420
229k
        if (std::isnan(x)) {
421
0
            return false;
422
0
        }
423
229k
        _unprocessed.emplace_back(x, w);
424
229k
        _unprocessed_weight += w;
425
229k
        processIfNecessary();
426
229k
        return true;
427
229k
    }
428
429
    void add(std::vector<Centroid>::const_iterator iter,
430
0
             std::vector<Centroid>::const_iterator end) {
431
0
        while (iter != end) {
432
0
            const size_t diff = std::distance(iter, end);
433
0
            const size_t room = _max_unprocessed - _unprocessed.size();
434
0
            auto mid = iter + std::min(diff, room);
435
0
            while (iter != mid) {
436
0
                _unprocessed.push_back(*(iter++));
437
0
            }
438
0
            if (_unprocessed.size() >= _max_unprocessed) {
439
0
                process();
440
0
            }
441
0
        }
442
0
    }
443
444
1.02k
    uint32_t serialized_size() {
445
1.02k
        return static_cast<uint32_t>(sizeof(uint32_t) + sizeof(Value) * 5 + sizeof(Index) * 2 +
446
1.02k
                                     sizeof(uint32_t) * 3 + _processed.size() * sizeof(Centroid) +
447
1.02k
                                     _unprocessed.size() * sizeof(Centroid) +
448
1.02k
                                     _cumulative.size() * sizeof(Weight));
449
1.02k
    }
450
451
493
    size_t serialize(uint8_t* writer) {
452
493
        uint8_t* dst = writer;
453
493
        uint32_t total_size = serialized_size();
454
493
        memcpy(writer, &total_size, sizeof(uint32_t));
455
493
        writer += sizeof(uint32_t);
456
493
        memcpy(writer, &_compression, sizeof(Value));
457
493
        writer += sizeof(Value);
458
493
        memcpy(writer, &_min, sizeof(Value));
459
493
        writer += sizeof(Value);
460
493
        memcpy(writer, &_max, sizeof(Value));
461
493
        writer += sizeof(Value);
462
493
        memcpy(writer, &_max_processed, sizeof(Index));
463
493
        writer += sizeof(Index);
464
493
        memcpy(writer, &_max_unprocessed, sizeof(Index));
465
493
        writer += sizeof(Index);
466
493
        memcpy(writer, &_processed_weight, sizeof(Value));
467
493
        writer += sizeof(Value);
468
493
        memcpy(writer, &_unprocessed_weight, sizeof(Value));
469
493
        writer += sizeof(Value);
470
471
493
        auto size = static_cast<uint32_t>(_processed.size());
472
493
        memcpy(writer, &size, sizeof(uint32_t));
473
493
        writer += sizeof(uint32_t);
474
377k
        for (int i = 0; i < size; i++) {
475
377k
            memcpy(writer, &_processed[i], sizeof(Centroid));
476
377k
            writer += sizeof(Centroid);
477
377k
        }
478
479
493
        size = static_cast<uint32_t>(_unprocessed.size());
480
493
        memcpy(writer, &size, sizeof(uint32_t));
481
493
        writer += sizeof(uint32_t);
482
        //TODO(weixiang): may be once memcpy is enough!
483
5.20k
        for (int i = 0; i < size; i++) {
484
4.70k
            memcpy(writer, &_unprocessed[i], sizeof(Centroid));
485
4.70k
            writer += sizeof(Centroid);
486
4.70k
        }
487
488
493
        size = static_cast<uint32_t>(_cumulative.size());
489
493
        memcpy(writer, &size, sizeof(uint32_t));
490
493
        writer += sizeof(uint32_t);
491
378k
        for (int i = 0; i < size; i++) {
492
377k
            memcpy(writer, &_cumulative[i], sizeof(Weight));
493
377k
            writer += sizeof(Weight);
494
377k
        }
495
493
        return writer - dst;
496
493
    }
497
498
476
    void unserialize(const uint8_t* type_reader) {
499
476
        uint32_t total_length = 0;
500
476
        memcpy(&total_length, type_reader, sizeof(uint32_t));
501
476
        type_reader += sizeof(uint32_t);
502
476
        memcpy(&_compression, type_reader, sizeof(Value));
503
476
        type_reader += sizeof(Value);
504
476
        memcpy(&_min, type_reader, sizeof(Value));
505
476
        type_reader += sizeof(Value);
506
476
        memcpy(&_max, type_reader, sizeof(Value));
507
476
        type_reader += sizeof(Value);
508
509
476
        memcpy(&_max_processed, type_reader, sizeof(Index));
510
476
        type_reader += sizeof(Index);
511
476
        memcpy(&_max_unprocessed, type_reader, sizeof(Index));
512
476
        type_reader += sizeof(Index);
513
476
        memcpy(&_processed_weight, type_reader, sizeof(Value));
514
476
        type_reader += sizeof(Value);
515
476
        memcpy(&_unprocessed_weight, type_reader, sizeof(Value));
516
476
        type_reader += sizeof(Value);
517
518
476
        uint32_t size;
519
476
        memcpy(&size, type_reader, sizeof(uint32_t));
520
476
        type_reader += sizeof(uint32_t);
521
476
        _processed.resize(size);
522
424k
        for (int i = 0; i < size; i++) {
523
423k
            memcpy(&_processed[i], type_reader, sizeof(Centroid));
524
423k
            type_reader += sizeof(Centroid);
525
423k
        }
526
476
        memcpy(&size, type_reader, sizeof(uint32_t));
527
476
        type_reader += sizeof(uint32_t);
528
476
        _unprocessed.resize(size);
529
5.08k
        for (int i = 0; i < size; i++) {
530
4.60k
            memcpy(&_unprocessed[i], type_reader, sizeof(Centroid));
531
4.60k
            type_reader += sizeof(Centroid);
532
4.60k
        }
533
476
        memcpy(&size, type_reader, sizeof(uint32_t));
534
476
        type_reader += sizeof(uint32_t);
535
476
        _cumulative.resize(size);
536
427k
        for (int i = 0; i < size; i++) {
537
426k
            memcpy(&_cumulative[i], type_reader, sizeof(Weight));
538
426k
            type_reader += sizeof(Weight);
539
426k
        }
540
476
    }
541
542
private:
543
    Value _compression;
544
545
    Value _min = std::numeric_limits<Value>::max();
546
547
    Value _max = std::numeric_limits<Value>::min();
548
549
    Index _max_processed;
550
551
    Index _max_unprocessed;
552
553
    Value _processed_weight = 0.0;
554
555
    Value _unprocessed_weight = 0.0;
556
557
    std::vector<Centroid> _processed;
558
559
    std::vector<Centroid> _unprocessed;
560
561
    std::vector<Weight> _cumulative;
562
563
    // return mean of i-th centroid
564
1.94k
    Value mean(int64_t i) const noexcept { return _processed[i].mean(); }
565
566
    // return weight of i-th centroid
567
146M
    Weight weight(int64_t i) const noexcept { return _processed[i].weight(); }
568
569
    // append all unprocessed centroids into current unprocessed vector
570
718
    void mergeUnprocessed(const std::vector<const TDigest*>& tdigests) {
571
718
        if (tdigests.size() == 0) {
572
0
            return;
573
0
        }
574
575
718
        size_t total = _unprocessed.size();
576
718
        for (const auto& td : tdigests) {
577
718
            total += td->_unprocessed.size();
578
718
        }
579
580
718
        _unprocessed.reserve(total);
581
718
        for (const auto& td : tdigests) {
582
718
            _unprocessed.insert(_unprocessed.end(), td->_unprocessed.cbegin(),
583
718
                                td->_unprocessed.cend());
584
718
            _unprocessed_weight += td->_unprocessed_weight;
585
718
        }
586
718
    }
587
588
    // merge all processed centroids together into a single sorted vector
589
717
    void mergeProcessed(const std::vector<const TDigest*>& tdigests) {
590
717
        if (tdigests.size() == 0) {
591
0
            return;
592
0
        }
593
594
717
        size_t total = 0;
595
717
        CentroidListQueue pq(CentroidListComparator {});
596
718
        for (const auto& td : tdigests) {
597
718
            const auto& sorted = td->_processed;
598
718
            auto size = sorted.size();
599
718
            if (size > 0) {
600
423
                pq.push(CentroidList(sorted));
601
423
                total += size;
602
423
                _processed_weight += td->_processed_weight;
603
423
            }
604
718
        }
605
717
        if (total == 0) {
606
295
            return;
607
295
        }
608
609
422
        if (_processed.size() > 0) {
610
422
            pq.push(CentroidList(_processed));
611
422
            total += _processed.size();
612
422
        }
613
614
422
        std::vector<Centroid> sorted;
615
422
        VLOG_CRITICAL << "total " << total;
616
422
        sorted.reserve(total);
617
618
2.04M
        while (!pq.empty()) {
619
2.04M
            auto best = pq.top();
620
2.04M
            pq.pop();
621
2.04M
            sorted.push_back(*(best.iter));
622
2.04M
            if (best.advance()) {
623
1.96M
                pq.push(best);
624
1.96M
            }
625
2.04M
        }
626
422
        _processed = std::move(sorted);
627
423
        if (_processed.size() > 0) {
628
423
            _min = std::min(_min, _processed[0].mean());
629
423
            _max = std::max(_max, (_processed.cend() - 1)->mean());
630
423
        }
631
422
    }
632
633
230k
    void processIfNecessary() {
634
230k
        if (isDirty()) {
635
423
            process();
636
423
        }
637
230k
    }
638
639
25.7k
    void updateCumulative() {
640
25.7k
        const auto n = _processed.size();
641
25.7k
        _cumulative.clear();
642
25.7k
        _cumulative.reserve(n + 1);
643
25.7k
        Weight previous = 0.0;
644
146M
        for (Index i = 0; i < n; i++) {
645
146M
            Weight current = weight(i);
646
146M
            auto halfCurrent = static_cast<Weight>(current / 2.0);
647
146M
            _cumulative.push_back(previous + halfCurrent);
648
146M
            previous = previous + current;
649
146M
        }
650
25.7k
        _cumulative.push_back(previous);
651
25.7k
    }
652
653
    // merges _unprocessed centroids and _processed centroids together and processes them
654
    // when complete, _unprocessed will be empty and _processed will have at most _max_processed centroids
655
25.0k
    void process() {
656
25.0k
        CentroidComparator cc;
657
        // select percentile_approx(lo_orderkey,0.5) from lineorder;
658
        // have test pdqsort and RadixSort, find here pdqsort performance is better when data is struct Centroid
659
        // But when sort plain type like int/float of std::vector<T>, find RadixSort is better
660
25.0k
        pdqsort(_unprocessed.begin(), _unprocessed.end(), cc);
661
25.0k
        auto count = _unprocessed.size();
662
25.0k
        _unprocessed.insert(_unprocessed.end(), _processed.cbegin(), _processed.cend());
663
25.0k
        std::inplace_merge(_unprocessed.begin(), _unprocessed.begin() + count, _unprocessed.end(),
664
25.0k
                           cc);
665
666
25.0k
        _processed_weight += _unprocessed_weight;
667
25.0k
        _unprocessed_weight = 0;
668
25.0k
        _processed.clear();
669
670
25.0k
        _processed.push_back(_unprocessed[0]);
671
25.0k
        Weight wSoFar = _unprocessed[0].weight();
672
25.0k
        Weight wLimit = _processed_weight * integratedQ(1.0);
673
674
25.0k
        auto end = _unprocessed.end();
675
146M
        for (auto iter = _unprocessed.cbegin() + 1; iter < end; iter++) {
676
146M
            const auto& centroid = *iter;
677
146M
            Weight projectedW = wSoFar + centroid.weight();
678
146M
            if (projectedW <= wLimit) {
679
1.07M
                wSoFar = projectedW;
680
1.07M
                (_processed.end() - 1)->add(centroid);
681
145M
            } else {
682
145M
                auto k1 = integratedLocation(wSoFar / _processed_weight);
683
145M
                wLimit = _processed_weight * integratedQ(static_cast<Value>(k1 + 1.0));
684
145M
                wSoFar += centroid.weight();
685
145M
                _processed.emplace_back(centroid);
686
145M
            }
687
146M
        }
688
25.0k
        _unprocessed.clear();
689
25.0k
        _min = std::min(_min, _processed[0].mean());
690
25.0k
        VLOG_CRITICAL << "new _min " << _min;
691
25.0k
        _max = std::max(_max, (_processed.cend() - 1)->mean());
692
25.0k
        VLOG_CRITICAL << "new _max " << _max;
693
25.0k
        updateCumulative();
694
25.0k
    }
695
696
0
    size_t checkWeights(const std::vector<Centroid>& sorted, Value total) {
697
0
        size_t badWeight = 0;
698
0
        auto k1 = 0.0;
699
0
        auto q = 0.0;
700
0
        for (auto iter = sorted.cbegin(); iter != sorted.cend(); iter++) {
701
0
            auto w = iter->weight();
702
0
            auto dq = w / total;
703
0
            auto k2 = integratedLocation(static_cast<Value>(q + dq));
704
0
            if (k2 - k1 > 1 && w != 1) {
705
0
                VLOG_CRITICAL << "Oversize centroid at " << std::distance(sorted.cbegin(), iter)
706
0
                              << " k1 " << k1 << " k2 " << k2 << " dk " << (k2 - k1) << " w " << w
707
0
                              << " q " << q;
708
0
                badWeight++;
709
0
            }
710
0
            if (k2 - k1 > 1.5 && w != 1) {
711
0
                VLOG_CRITICAL << "Egregiously Oversize centroid at "
712
0
                              << std::distance(sorted.cbegin(), iter) << " k1 " << k1 << " k2 "
713
0
                              << k2 << " dk " << (k2 - k1) << " w " << w << " q " << q;
714
0
                badWeight++;
715
0
            }
716
0
            q += dq;
717
0
            k1 = k2;
718
0
        }
719
0
720
0
        return badWeight;
721
0
    }
722
723
    /**
724
    * Converts a quantile into a centroid scale value.  The centroid scale is nomin_ally
725
    * the number k of the centroid that a quantile point q should belong to.  Due to
726
    * round-offs, however, we can't align things perfectly without splitting points
727
    * and sorted.  We don't want to do that, so we have to allow for offsets.
728
    * In the end, the criterion is that any quantile range that spans a centroid
729
    * scale range more than one should be split across more than one centroid if
730
    * possible.  This won't be possible if the quantile range refers to a single point
731
    * or an already existing centroid.
732
    * <p/>
733
    * This mapping is steep near q=0 or q=1 so each centroid there will correspond to
734
    * less q range.  Near q=0.5, the mapping is flatter so that sorted there will
735
    * represent a larger chunk of quantiles.
736
    *
737
    * @param q The quantile scale value to be mapped.
738
    * @return The centroid scale value corresponding to q.
739
    */
740
145M
    Value integratedLocation(Value q) const {
741
145M
        return static_cast<Value>(_compression * (std::asin(2.0 * q - 1.0) + M_PI / 2) / M_PI);
742
145M
    }
743
744
145M
    Value integratedQ(Value k) const {
745
145M
        return static_cast<Value>(
746
145M
                (std::sin(std::min(k, _compression) * M_PI / _compression - M_PI / 2) + 1) / 2);
747
145M
    }
748
749
    /**
750
     * Same as {@link #weightedAverageSorted(Value, Value, Value, Value)} but flips
751
     * the order of the variables if <code>x2</code> is greater than
752
     * <code>x1</code>.
753
    */
754
1.21k
    static Value weightedAverage(Value x1, Value w1, Value x2, Value w2) {
755
1.21k
        return (x1 <= x2) ? weightedAverageSorted(x1, w1, x2, w2)
756
1.21k
                          : weightedAverageSorted(x2, w2, x1, w1);
757
1.21k
    }
758
759
    /**
760
    * Compute the weighted average between <code>x1</code> with a weight of
761
    * <code>w1</code> and <code>x2</code> with a weight of <code>w2</code>.
762
    * This expects <code>x1</code> to be less than or equal to <code>x2</code>
763
    * and is guaranteed to return a number between <code>x1</code> and
764
    * <code>x2</code>.
765
    */
766
1.21k
    static Value weightedAverageSorted(Value x1, Value w1, Value x2, Value w2) {
767
1.21k
        DCHECK_LE(x1, x2);
768
1.21k
        const Value x = (x1 * w1 + x2 * w2) / (w1 + w2);
769
1.21k
        return std::max(x1, std::min(x, x2));
770
1.21k
    }
771
772
0
    static Value interpolate(Value x, Value x0, Value x1) { return (x - x0) / (x1 - x0); }
773
774
    /**
775
    * Computes an interpolated value of a quantile that is between two sorted.
776
    *
777
    * Index is the quantile desired multiplied by the total number of samples - 1.
778
    *
779
    * @param index              Denormalized quantile desired
780
    * @param previousIndex      The denormalized quantile corresponding to the center of the previous centroid.
781
    * @param nextIndex          The denormalized quantile corresponding to the center of the following centroid.
782
    * @param previousMean       The mean of the previous centroid.
783
    * @param nextMean           The mean of the following centroid.
784
    * @return  The interpolated mean.
785
    */
786
    static Value quantile(Value index, Value previousIndex, Value nextIndex, Value previousMean,
787
0
                          Value nextMean) {
788
0
        const auto delta = nextIndex - previousIndex;
789
0
        const auto previousWeight = (nextIndex - index) / delta;
790
0
        const auto nextWeight = (index - previousIndex) / delta;
791
0
        return previousMean * previousWeight + nextMean * nextWeight;
792
0
    }
793
};
794
#include "common/compile_check_end.h"
795
} // namespace doris