/* * tdigest.c - implementation of t-digest for PostgreSQL, useful for estimation * of quantiles, percentiles, trimmed means, and various similar metrics. * * Copyright (C) Tomas Vondra, 2019 */ #include #include #include #include #include #include #include #include #include "postgres.h" #include "common/int.h" #include "libpq/pqformat.h" #include "miscadmin.h" #include "utils/array.h" #include "utils/builtins.h" #include "utils/memutils.h" #if PG_VERSION_NUM >= 120000 #include "utils/float.h" /* float8out_internal */ #endif #include "catalog/pg_type.h" PG_MODULE_MAGIC; /* * A centroid, used both for in-memory and on-disk storage. */ typedef struct centroid_t { double mean; int64 count; } centroid_t; /* * On-disk representation of the t-digest. */ typedef struct tdigest_t { int32 vl_len_; /* varlena header (do not touch directly!) */ int32 flags; /* on-disk format flags */ int64 count; /* number of items added to the t-digest */ int compression; /* compression used to build the digest */ int ncentroids; /* number of centroids in the array */ centroid_t centroids[FLEXIBLE_ARRAY_MEMBER]; } tdigest_t; /* * Centroids used to store (sum,count), but we want to store (mean,count) * because that allows us to prevent rounding errors e.g. when merging * centroids with the same mean, or adding the same value to the centroid. * * To handle existing tdigest data in backwards-compatible way, we have * a flag marking the new ones with mean, and we convert the old values. */ #define TDIGEST_STORES_MEAN 0x0001 /* All valid flags, OR-ed. */ #define TDIGEST_VALID_FLAGS (TDIGEST_STORES_MEAN) /* * An aggregate state, representing the t-digest and some additional info * (requested percentiles, ...). * * When adding new values to the t-digest, we add them as centroids into a * separate "uncompacted" part of the array. While centroids need more space * than plain points (16B vs. 8B), making the aggregate state quite a bit * larger, it does simplify the code quite a bit as it only needs to deal * with a single struct type instead of two (centroids + points). But maybe * we should separate those two things in the future. * * XXX We only ever use one of values/percentiles, never both at the same * time. In the future the values may use a different data type than double * (e.g. numeric), so we keep both fields. */ typedef struct tdigest_aggstate_t { /* basic t-digest fields (centroids at the end) */ int64 count; /* number of samples in the digest */ int ncompactions; /* number of merges/compactions */ int compression; /* compression parameter */ int maxcentroids; /* capacity of the centroids buffer */ int ncentroids; /* number of centroids */ int ncompacted; /* compacted part */ /* array of requested percentiles and values */ int npercentiles; /* number of percentiles */ int nvalues; /* number of values */ double trim_low; /* low threshold (for trimmed aggs) */ double trim_high; /* high threshold (for trimmed aggs) */ double *percentiles; /* array of percentiles (if any) */ double *values; /* array of values (if any) */ centroid_t *centroids; /* centroids for the digest */ } tdigest_aggstate_t; static int centroid_cmp(const void *a, const void *b); /* * Detoast a tdigest, making sure the result is properly aligned. * * The tdigest type is int4-aligned - that's the default for CREATE TYPE, and * changing it now would break on-disk compatibility, as typalign determines * where the value is placed in a tuple. But tdigest_t contains int64 and * double fields, which can require 8B alignment. Reading digest->count or * digest->centroids[i].mean from a value that is merely 4B aligned is * undefined behavior - it happens to work on x86-64, but it's a crash on * platforms not allowing unaligned access, and UBSan complains about it. * * Detoasting out-of-line, compressed and short-header values already makes * a freshly palloc-ed (and thus properly aligned) copy. But an inline value * with a 4B varlena header may still point into the tuple, so we have to make * an aligned copy ourselves if necessary. */ static tdigest_t * tdigest_detoast(Datum datum) { tdigest_t *digest = (tdigest_t *) PG_DETOAST_DATUM(datum); /* * Detoasting always produces a 4B header, so VARSIZE is fine here, and * reading it only needs the 4B alignment the type already guarantees. */ if (!PointerIsAligned(digest, int64)) { tdigest_t *aligned = (tdigest_t *) palloc(VARSIZE(digest)); memcpy(aligned, digest, VARSIZE(digest)); digest = aligned; } return digest; } #define PG_GETARG_TDIGEST(x) tdigest_detoast(PG_GETARG_DATUM(x)) /* * Size of buffer for incoming data, as a multiple of the compression value. * Quoting from the t-digest paper: * * The constant of proportionality should be determined by experiment, but * micro-benchmarks indicate that C2/C1 is in the range from 5 to 20 for * a single core of an Intel i7 processor. In these micro-benchmarks, * increasing the buffer size to (10 * delta) dramatically improves the * average speed but further buffer size increases have much less effect. * * XXX Maybe make the coefficient user-defined, with some reasonable limits * (say 2 - 20), so that users can pick the right trade-off between speed * and memory usage. */ #define BUFFER_SIZE(compression) (10 * (compression)) #define BUFFER_INITIAL_SIZE 16 #define AssertBounds(index, length) Assert((index) >= 0 && (index) < (length)) #define MIN_COMPRESSION 10 #define MAX_COMPRESSION 10000 /* prototypes */ PG_FUNCTION_INFO_V1(tdigest_add_double_array); PG_FUNCTION_INFO_V1(tdigest_add_double_array_count); PG_FUNCTION_INFO_V1(tdigest_add_double_array_values); PG_FUNCTION_INFO_V1(tdigest_add_double_array_values_count); PG_FUNCTION_INFO_V1(tdigest_add_double); PG_FUNCTION_INFO_V1(tdigest_add_double_count); PG_FUNCTION_INFO_V1(tdigest_add_double_values); PG_FUNCTION_INFO_V1(tdigest_add_double_values_count); PG_FUNCTION_INFO_V1(tdigest_add_digest_array); PG_FUNCTION_INFO_V1(tdigest_add_digest_array_values); PG_FUNCTION_INFO_V1(tdigest_add_digest); PG_FUNCTION_INFO_V1(tdigest_add_digest_values); PG_FUNCTION_INFO_V1(tdigest_array_percentiles); PG_FUNCTION_INFO_V1(tdigest_array_percentiles_of); PG_FUNCTION_INFO_V1(tdigest_percentiles); PG_FUNCTION_INFO_V1(tdigest_percentiles_of); PG_FUNCTION_INFO_V1(tdigest_digest); PG_FUNCTION_INFO_V1(tdigest_serial); PG_FUNCTION_INFO_V1(tdigest_deserial); PG_FUNCTION_INFO_V1(tdigest_combine); PG_FUNCTION_INFO_V1(tdigest_in); PG_FUNCTION_INFO_V1(tdigest_out); PG_FUNCTION_INFO_V1(tdigest_send); PG_FUNCTION_INFO_V1(tdigest_recv); PG_FUNCTION_INFO_V1(tdigest_is_valid); PG_FUNCTION_INFO_V1(tdigest_count); PG_FUNCTION_INFO_V1(tdigest_to_json); PG_FUNCTION_INFO_V1(tdigest_to_array); PG_FUNCTION_INFO_V1(tdigest_add_double_increment); PG_FUNCTION_INFO_V1(tdigest_add_double_array_increment); PG_FUNCTION_INFO_V1(tdigest_union_double_increment); PG_FUNCTION_INFO_V1(tdigest_add_double_trimmed); PG_FUNCTION_INFO_V1(tdigest_add_double_count_trimmed); PG_FUNCTION_INFO_V1(tdigest_add_digest_trimmed); PG_FUNCTION_INFO_V1(tdigest_trimmed_avg); PG_FUNCTION_INFO_V1(tdigest_trimmed_sum); PG_FUNCTION_INFO_V1(tdigest_digest_sum); PG_FUNCTION_INFO_V1(tdigest_digest_avg); Datum tdigest_add_double_array(PG_FUNCTION_ARGS); Datum tdigest_add_double_array_count(PG_FUNCTION_ARGS); Datum tdigest_add_double_array_values(PG_FUNCTION_ARGS); Datum tdigest_add_double_array_values_count(PG_FUNCTION_ARGS); Datum tdigest_add_double(PG_FUNCTION_ARGS); Datum tdigest_add_double_count(PG_FUNCTION_ARGS); Datum tdigest_add_double_values(PG_FUNCTION_ARGS); Datum tdigest_add_double_values_count(PG_FUNCTION_ARGS); Datum tdigest_add_digest_array(PG_FUNCTION_ARGS); Datum tdigest_add_digest_array_values(PG_FUNCTION_ARGS); Datum tdigest_add_digest(PG_FUNCTION_ARGS); Datum tdigest_add_digest_values(PG_FUNCTION_ARGS); Datum tdigest_array_percentiles(PG_FUNCTION_ARGS); Datum tdigest_array_percentiles_of(PG_FUNCTION_ARGS); Datum tdigest_percentiles(PG_FUNCTION_ARGS); Datum tdigest_percentiles_of(PG_FUNCTION_ARGS); Datum tdigest_digest(PG_FUNCTION_ARGS); Datum tdigest_serial(PG_FUNCTION_ARGS); Datum tdigest_deserial(PG_FUNCTION_ARGS); Datum tdigest_combine(PG_FUNCTION_ARGS); Datum tdigest_in(PG_FUNCTION_ARGS); Datum tdigest_out(PG_FUNCTION_ARGS); Datum tdigest_send(PG_FUNCTION_ARGS); Datum tdigest_recv(PG_FUNCTION_ARGS); Datum tdigest_is_valid(PG_FUNCTION_ARGS); Datum tdigest_count(PG_FUNCTION_ARGS); Datum tdigest_add_double_increment(PG_FUNCTION_ARGS); Datum tdigest_add_double_array_increment(PG_FUNCTION_ARGS); Datum tdigest_union_double_increment(PG_FUNCTION_ARGS); Datum tdigest_to_json(PG_FUNCTION_ARGS); Datum tdigest_to_array(PG_FUNCTION_ARGS); Datum tdigest_add_double_trimmed(PG_FUNCTION_ARGS); Datum tdigest_add_double_count_trimmed(PG_FUNCTION_ARGS); Datum tdigest_add_digest_trimmed(PG_FUNCTION_ARGS); Datum tdigest_trimmed_avg(PG_FUNCTION_ARGS); Datum tdigest_trimmed_sum(PG_FUNCTION_ARGS); Datum tdigest_digest_sum(PG_FUNCTION_ARGS); Datum tdigest_digest_avg(PG_FUNCTION_ARGS); static ArrayType *double_array_allocate(int nitems); static const double *array_to_double(ArrayType *v, const char *what, int *len); static int64 double_to_int64(double value, int64 maxvalue); static tdigest_aggstate_t *tdigest_copy(tdigest_aggstate_t *state); static void tdigest_aggstate_enlarge(tdigest_aggstate_t *state); static void tdigest_aggstate_shrink(tdigest_aggstate_t *state); #if PG_VERSION_NUM < 150000 /* * Thin wrappers that convert strings to exactly 64-bit integers, matching our * definition of int64. (For the naming, compare that POSIX has * strtoimax()/strtoumax() which return intmax_t/uintmax_t.) * * XXX Backward compatibility */ #ifdef HAVE_LONG_INT_64 /* int64 is "long int", so strtol() returns exactly the right width */ #define strtoi64(str, endptr, base) ((int64) strtol(str, endptr, base)) #else /* int64 is "long long int" (C99 guarantees it is at least 64 bits) */ #define strtoi64(str, endptr, base) ((int64) strtoll(str, endptr, base)) #endif #endif /* PG_VERSION_NUM < 150000 */ /* basic checks on the t-digest (proper sum of counts, ...) */ static void AssertCheckTDigest(tdigest_t *digest) { #ifdef USE_ASSERT_CHECKING int i; int64 cnt; Assert(digest->flags == 0 || digest->flags == TDIGEST_STORES_MEAN); Assert((digest->compression >= MIN_COMPRESSION) && (digest->compression <= MAX_COMPRESSION)); Assert(digest->count >= 0); Assert(digest->ncentroids >= 0); Assert(digest->ncentroids <= BUFFER_SIZE(digest->compression)); cnt = 0; for (i = 0; i < digest->ncentroids; i++) { Assert(digest->centroids[i].count > 0); Assert(isfinite(digest->centroids[i].mean)); cnt += digest->centroids[i].count; /* FIXME also check this does work with the scale function */ } Assert(VARSIZE_ANY(digest) == offsetof(tdigest_t, centroids) + digest->ncentroids * sizeof(centroid_t)); Assert(digest->count == cnt); #endif } static void AssertCheckTDigestAggState(tdigest_aggstate_t *state) { #ifdef USE_ASSERT_CHECKING int i; int64 cnt; Assert(state->npercentiles >= 0); Assert(((state->npercentiles == 0) && (state->percentiles == NULL)) || ((state->npercentiles > 0) && (state->percentiles != NULL))); for (i = 0; i < state->npercentiles; i++) Assert((state->percentiles[i] >= 0.0) && (state->percentiles[i] <= 1.0)); Assert((state->compression >= MIN_COMPRESSION) && (state->compression <= MAX_COMPRESSION)); Assert(state->count >= 0); Assert(state->ncentroids >= 0); Assert(state->ncentroids <= state->maxcentroids); Assert(state->maxcentroids <= BUFFER_SIZE(state->compression)); Assert(state->ncentroids <= BUFFER_SIZE(state->compression)); cnt = 0; for (i = 0; i < state->ncentroids; i++) { Assert(state->centroids[i].count > 0); Assert(isfinite(state->centroids[i].mean)); cnt += state->centroids[i].count; /* XXX maybe check this does work with the scale function */ } Assert(state->count == cnt); #endif } static void reverse_centroids(centroid_t *centroids, int ncentroids) { int start = 0, end = (ncentroids - 1); while (start < end) { centroid_t tmp = centroids[start]; centroids[start] = centroids[end]; centroids[end] = tmp; start++; end--; } } static void rebalance_centroids(centroid_t *centroids, int ncentroids, int64 weight_before, int64 weight_after) { double ratio = weight_before / (double) weight_after; int64 count_before = 0; int64 count_after = 0; int start = 0; int end = (ncentroids - 1); int i; centroid_t *scratch = palloc(sizeof(centroid_t) * ncentroids); Assert(weight_after > 0); i = 0; while (i < ncentroids) { while (i < ncentroids) { scratch[start] = centroids[i]; count_before += centroids[i].count; i++; start++; if (count_before > count_after * ratio) break; } while (i < ncentroids) { scratch[end] = centroids[i]; count_after += centroids[i].count; i++; end--; if (count_before < count_after * ratio) break; } } memcpy(centroids, scratch, sizeof(centroid_t) * ncentroids); pfree(scratch); } /* * Sort an array of centroids, with the given total count. * * We have to sort the whole array, because we don't just simply sort the * centroids - we do the rebalancing of items with the same mean too. */ static void tdigest_sort_centroids(centroid_t *centroids, int ncentroids, int64 count) { int i; int64 count_so_far; int64 next_group; int64 median_count; /* sort the whole array before rebalancing equal-mean groups */ pg_qsort(centroids, ncentroids, sizeof(centroid_t), centroid_cmp); /* * The centroids are sorted by (mean,count). That's fine for centroids up * to median, but above median this ordering is incorrect for centroids * with the same mean (or for groups crossing the median boundary). To fix * this we 'rebalance' those groups. Those entirely above median can be * simply sorted in the opposite order, while those crossing the median * need to be rebalanced depending on what part is below/above median. */ count_so_far = 0; next_group = 0; /* includes count_so_far */ median_count = (count / 2); /* * Split the centroids into groups with the same mean, process each group * depending on whether it falls before/after median. */ i = 0; while (i < ncentroids) { int j = i; int group_size = 0; CHECK_FOR_INTERRUPTS(); /* * Consume the first centroid of the group unconditionally. It is what * defines the group, so comparing it against itself decides nothing, * and for a NaN mean the comparison would be false, leaving "j" at "i" * and the outer loop without any way to advance. This way guarantees * forward progress. */ next_group += centroids[j].count; group_size++; j++; /* determine the end of the group */ while ((j < ncentroids) && (centroids[i].mean == centroids[j].mean)) { next_group += centroids[j].count; group_size++; j++; } /* * We can ignore groups of size 1 (number of centroids, not counts), as * those are trivially sorted. */ if (group_size > 1) { if (count_so_far >= median_count) { /* group fully above median - reverse the order */ reverse_centroids(¢roids[i], group_size); } else if (next_group > median_count) /* group split by median */ { rebalance_centroids(¢roids[i], group_size, median_count - count_so_far, next_group - median_count); } } /* * We should be making forward progress. If not, we're in an infinite * loop. With properly formed digests that should not happen. */ Assert(i < j); i = j; count_so_far = next_group; } } /* * Sort centroids in the aggregate state. */ static void tdigest_sort(tdigest_aggstate_t *state) { tdigest_sort_centroids(state->centroids, state->ncentroids, state->count); } /* * Perform compaction of the t-digest, i.e. merge the centroids as required * by the compression parameter. * * We always keep the data sorted in ascending order. This way we can reuse * the sort between compactions, and also when computing the quantiles. * * The regular compaction is not guaranteed to make any progress. The size * limits are calculated in double, and may end up too low to allow merging any * centroids. The limits are computed from exact integer remainders, which * makes this very unlikely, but if a compaction still leaves the buffer full, * we raise an error rather than continue with a digest that has no room left. * * XXX Switch the direction regularly, to eliminate possible bias and improve * accuracy, as mentioned in the paper. * * XXX This initially used the k1 scale function, but the implementation was * not limiting the number of centroids for some reason (it might have been * a bug in the implementation, of course). The current code is a modified * copy from ajwerner [1], and AFAIK it's the k2 function, it's much simpler * and generally works quite nicely. * * [1] https://github.com/ajwerner/tdigestc/blob/master/go/tdigest.c */ static void tdigest_compact(tdigest_aggstate_t *state) { int i; int cur; /* current centroid */ int64 count_so_far; int64 total_count; double denom; double normalizer; int start; int step; int n; AssertCheckTDigestAggState(state); /* if the digest is fully compacted, it's been already compacted */ if (state->ncompacted == state->ncentroids) return; /* * If there's just a single centroid, there's nothing to compact (or * sort). And we'd also end up with a division by zero below, because * log(1) = 0. It'd work out in the end, because 1/0 = infinity, and * so we'd merge nothing. But it's sloppy. * * XXX We're checking ncentroids, while the log() is on total_count. * But that's fine. Compaction/sort is pointless no matter how large * the single centroid is. And with 2+ centroids, the total_count has * to be 2+ too. */ if (state->ncentroids == 1) { state->ncompacted = state->ncentroids; return; } tdigest_sort(state); state->ncompactions++; if (state->ncompactions % 2 == 0) { start = 0; step = 1; } else { start = state->ncentroids - 1; step = -1; } total_count = state->count; denom = 2 * M_PI * total_count * log(total_count); normalizer = state->compression / denom; cur = start; count_so_far = 0; n = 1; for (i = start + step; (i >= 0) && (i < state->ncentroids); i += step) { int64 proposed_count; double q0; double q2; double z; bool should_add; CHECK_FOR_INTERRUPTS(); proposed_count = state->centroids[cur].count + state->centroids[i].count; z = proposed_count * normalizer; q0 = count_so_far / (double) total_count; q2 = (count_so_far + proposed_count) / (double) total_count; /* * Calculate the (1 - q) factors from the exact integer remainders, * instead of subtracting the quotients from 1. * * The two are equivalent with exact arithmetic, but not in double. If * a single centroid holds almost the whole weight of the digest, the * quotient rounds to exactly 1.0, and (1 - q) cancels to exactly 0. * The size limit then says no two centroids may be merged, even * though the exact limit is small but positive - and the compaction * ends up not making any progress at all. */ should_add = (z <= (q0 * ((double) (total_count - count_so_far) / (double) total_count))) && (z <= (q2 * ((double) (total_count - count_so_far - proposed_count) / (double) total_count))); if (should_add) { /* * If both centroids have the same mean, don't calculate it again. * The recalculation may cause rounding errors, so that the means * would drift apart over time. We want to keep them equal for as * long as possible. */ if (state->centroids[cur].mean != state->centroids[i].mean) { double mean; int64 count; /* count can't overflow int64 (total is within INT64_MAX) */ count = state->centroids[i].count; count += state->centroids[cur].count; /* calculate the mean in a way that should not overflow */ mean = state->centroids[i].mean * (state->centroids[i].count / (double) count); mean += state->centroids[cur].mean * (state->centroids[cur].count / (double) count); /* should not happen for finite inputs */ Assert(!isnan(mean)); /* * paranoia: clamp to not underflow/overflow the inputs * * We may be walking the centroids forward or backwards, which * determines whether (cur < i) or (cur > i). */ if (step > 0) { Assert(cur < i); mean = Min(Max(state->centroids[cur].mean, mean), state->centroids[i].mean); } else { Assert(cur > i); mean = Min(Max(state->centroids[i].mean, mean), state->centroids[cur].mean); } state->centroids[cur].mean = mean; } /* XXX Do this after possibly recalculating the mean. */ state->centroids[cur].count += state->centroids[i].count; } else { count_so_far += state->centroids[cur].count; cur += step; n++; state->centroids[cur] = state->centroids[i]; } if (cur != i) { state->centroids[i].count = 0; state->centroids[i].mean = 0; } } state->ncentroids = n; state->ncompacted = state->ncentroids; if (step < 0) memmove(state->centroids, &state->centroids[cur], n * sizeof(centroid_t)); /* * The compaction should have reduced the number of digest centroids. If * that did not succeed, and the buffer is still full, something must have * gone wrong. It's not safe to proceed. * * XXX Could be useful to include the digest through errdetail(), but with * large digests it might generate a significant log volume. */ if (state->ncentroids == BUFFER_SIZE(state->compression)) ereport(ERROR, (errcode(ERRCODE_INTERNAL_ERROR), errmsg("digest compaction failed to reduce number of centroids"), errhint("This should be impossible. Please report this to maintainers " "of the extension (ideally with a reproducer)."))); /* Maybe reclaim some of the centroid buffer. */ tdigest_aggstate_shrink(state); AssertCheckTDigestAggState(state); /* Must have freed some space in the buffer. */ Assert(state->ncentroids < BUFFER_SIZE(state->compression)); } /* * Interpolate between ordered finite endpoints without losing monotonicity. * * For endpoints of the same sign, separately rounded weighted products can * make (1 - q) * a + q * b decrease as q increases. Interpolating the * difference avoids that. * * Anchor negative endpoints at b to avoid subtracting nearly equal large * magnitudes when q is close to 1. That's necessary because of widely * separated negative means, like centroids at -2^54 and -1. Subtracting * almost equal large magnitudes exposes the rounding error in (b - a). * The right-anchored formula (b - (1 - q) * (b - a)) solves that. It's * monotone and does not have the cancellation issue. * * Positive intervals keep the existing left-anchored formula. * * For opposite signs, use the overflow-safe weighted sum. The difference may * overflow, but both weighted terms are nondecreasing, so use that. * * Return the endpoints explicitly when q is zero or one. A rounded difference * may not reconstruct either endpoint, and clamping cannot repair a result * that remains inside the bracket. */ static double tdigest_interpolate(double a, double b, double q) { double result; Assert(isfinite(a) && isfinite(b) && a <= b); Assert(q >= 0.0 && q <= 1.0); /* Rounded differences need not reproduce the endpoints exactly. */ if (q == 0.0) return a; if (q == 1.0) return b; /* * There are two approaches to calculating linear interpolation: * * (1 - q) * a + q * b and a + q * (b - a) * * The first formula is monotonic only when (a * b < 0). We know that * (a <= b), so this means (a < 0) and (b > 0). We use it also when * either value is 0. * * The second formula is monotonic, but does not guarantee the result * to be "b" when "q = 1", due to floating-point arithmetic errors. It * also has a risk of cancellation when subtracting almost equal large * magnitudes. To deal with that, we use a variation anchored to the * right endpoint. * * We switch between the options, to pick the better one. */ if (a <= 0.0 && b >= 0.0) result = (1 - q) * a + q * b; else if (b < 0.0) /* anchored to right endpoint */ result = b - (1 - q) * (b - a); else /* anchored to left endpoint */ result = a + q * (b - a); /* Final rounding must not move the result outside the endpoints. */ return Max(a, Min(b, result)); } /* * Estimate requested quantiles from the t-digest agg state. */ static void tdigest_compute_quantiles(tdigest_aggstate_t *state, double *result) { int i, j; AssertCheckTDigestAggState(state); /* * Trigger a compaction, which also sorts the data. * * XXX It might seem we can just do a sort here, and maybe get a bit more * accurate results. However, sorting single-item centroids in between much * larger centroids seems to interfere with the interpolation, making the * errors much larger. */ tdigest_compact(state); /* * Determine the two centroids the quantile lies between, and calculate the * estimate using linear interpolation. * * XXX All of this works fine for t-digests with non-extreme counts, up to * about 2^52. At that point the double precision ULP gets > 1.0, and some * of the calculations here start misbehaving a little. For example the * (count * 0.9999...) can get higher than count, etc. We try to prevent * obviously bogus results, but it's futile to try to fix this perfectly. * The cases are extremely rare, and we're calculating estimates anyway. * If we wanted to fix this properly, we'd need to use some sort of large * float data type (there seems to be "long double" and binary128). * * XXX The rounding/precision issues affect only accuracy of results, not * correctness of the code. For example, it must not result in OOB access * to bogus centroids etc. */ for (i = 0; i < state->npercentiles; i++) { int64 count; double goal = (state->count * state->percentiles[i]); bool is_before = false; centroid_t *c = NULL, *prev, *next; /* integer and fractional parts of half-centroids before/after */ double distance, total_distance, q; /* first centroid for percentile 0.0 */ if (state->percentiles[i] == 0.0) { c = &state->centroids[0]; result[i] = c->mean; continue; } /* last centroid for percentile 1.0 */ if (state->percentiles[i] == 1.0) { c = &state->centroids[state->ncentroids - 1]; result[i] = c->mean; continue; } /* * Walk the centroids and calculate running sum of counts. Stop before * adding a centroid that would exceed the goal - we don't know if the * goal falls before/after the mean yet. * * FIXME There can be multiple centroids with the same mean, in which * case we should use the total count for all of them. Not sure how * likely it's to have centroids with exactly the same mean. But it * might affect the interpolation later. */ count = 0; for (j = 0; j < state->ncentroids; j++) { CHECK_FOR_INTERRUPTS(); c = &state->centroids[j]; /* Adding the centroid would exceed the goal, so stop. */ if (count + c->count >= goal) break; count += c->count; } /* * Adding the whole centroid would exceed the goal, but we don't know * on which side of the mean the value lies yet. We might have also * hit the mean exactly. Let's figure that out. * * This will determine which centroids we'll look at for linear * interpolation (previous/following one) later. * * We know centroid "c" exceeds the goal, but did we hit the mean, * or are we to the left/right? We assume half the items are before * the mean, half after. */ is_before = goal < (count + c->count / 2.0); /* * Pick centroids for linear interpolation, depending on which side * of the "current" centroid we fell on. Either use the previous or * the following centroid. * * For extreme percentile values (or somehow weird digests) we can * end up before/after the last centroid, in which case we need to * be careful to not access OOB. */ if (is_before) { /* no previous centroid, use the current (first) one */ if (j == 0) { result[i] = c->mean; continue; } prev = &state->centroids[j - 1]; next = &state->centroids[j]; Assert(next == c); /* * Undo the centroid already added above (count is integer, * so we can't undo half of it without possibly losing half * of the count). We'll deal with that later. */ count -= prev->count; } else { /* no following centroid, use the current (last) one */ if (j == (state->ncentroids - 1)) { result[i] = c->mean; continue; } prev = &state->centroids[j]; next = &state->centroids[j + 1]; Assert(prev == c); } /* paranoia: make sure the prev/next centroids are valid */ Assert(prev >= &state->centroids[0]); Assert(next <= &state->centroids[state->ncentroids - 1]); Assert((prev + 1) == next); /* * Now we know the quantile lies somewhere between the centroids, * we need to calculate the correct value. (We know it's not at * either mean, that's what the above branches are for.) * * We will calculate the distance from the first mean, the total * distance between the means. And we'll do linear interpolation. */ /* distance to the first mean (of the previous centroid) */ distance = (double) (goal - count) - prev->count / 2.0; /* distance between the means of the two centroids */ total_distance = (prev->count / 2.0) + (next->count / 2.0); /* * We should be "to the right" of the first centroid, and should not * be so far ahead to exceed the next one. So in principle, this * should be true: * * Assert((distance >= 0) && (distance <= total_distance)); * * But, it's tricky due to precision and rounding. We're switching * from int64 to double, and double has much lower precision close * to INT64_MAX (ULP >> 1.0). With high goal and/or count values we * can end up with distance outside the [0, total_distance] range, * or just hit the centroids exactly. * * XXX Try uncommenting the assert, there's a test triggering it. */ /* * Clamp distance to [0, total_distance], to mitigate unexpected * rounding / precision errors. */ distance = Max(0.0, Min(total_distance, distance)); /* Keep the interpolation fraction within the centroid bracket. */ q = Max(0.0, Min(1.0, distance / total_distance)); result[i] = tdigest_interpolate(prev->mean, next->mean, q); } } /* * Calculate the rank at the midpoint of a centroid or equal-mean group. * * Form twice the midpoint exactly before converting to double. Rounding the * prefix and half-weight separately can reverse the order of nearby ranks. * Since count_before + count_at <= total_count <= INT64_MAX, the doubled * midpoint fits in uint64, even when it would overflow int64. */ static double tdigest_compute_rank(int64 count_before, int64 count_at, int64 total_count) { uint64 count_twice; double rank; Assert(count_at > 0 && count_at <= total_count); Assert(count_before >= 0 && count_before <= total_count - count_at); count_twice = 2 * (uint64) count_before + (uint64) count_at; rank = ((double) count_twice / 2.0) / (double) total_count; return Max(0.0, Min(1.0, rank)); } /* * Estimate inverse quantiles for values using a t-digest agg state. * * Essentially an inverse to tdigest_compute_quantiles. * * XXX I suspect this may not be doing the interpolation quite right when * there are multiple centroids with the same mean (compute_quantiles has * essentially the same issue in the opposite direction, but it's easier * to describe it for rank). The current code walks centroids until we find * a centroid with a higher mean - let's assume the value is between two * centroids (not equal to any centroid mean). It takes the previous and * following centroids, and interpolates those. * * But there can be multiple centroids with the same mean and different * counts. If the value falls before the median, we could see e.g. this: * * (1.0, 1000), (2.0, 1), (2.0, 999) * * Both means have the same weight, but the interprolation will use only * the (1.0, 1000) and (2.0, 1) centroids, i.e. 2.0 will have much less * weight. And then at 2.0 the CDF jumps up suddenly. * * This may not be very common, due to merging/compaction. We usually merge * the small centroids into the large ones quickly, so the differences * tend to be less extreme, which limits the effect on interpolation. Also, * the merging does floating point math with rounding errors, so means get * to diverge at some point. But it can happen e.g. when merging digests. * * There's also a bit of inconsistency (both here and in the reference * implementation), because we sum all the counts for the case when the * value matches a mean exactly, but not in the other cases. */ static void tdigest_compute_quantiles_of(tdigest_aggstate_t *state, double *result) { int i; AssertCheckTDigestAggState(state); /* * Trigger a compaction, which also sorts the data. * * XXX It might seem we can just do a sort here, and maybe get a bit more * accurate results. However, sorting single-item centroids in between much * larger centroids seems to interfere with the interpolation, making the * errors much larger. */ tdigest_compact(state); for (i = 0; i < state->nvalues; i++) { int j; int64 count; double value = state->values[i]; double d, q, q1, q2; /* next and previous centroids */ centroid_t *curr = NULL; centroid_t *prev = NULL; /* handle infinity/NaN values by mapping them to 0.0, 1.0 and NaN */ if (!isfinite(value)) { if (isnan(value)) result[i] = NAN; else if (value < 0) /* -infinity */ result[i] = 0.0; else /* infinity */ result[i] = 1.0; continue; } /* * Find the first centroid with (mean >= value), and remember the * last centroid before that - if the value is in between, we will * be calculating the percentile by linear approximation. */ count = 0; for (j = 0; j < state->ncentroids; j++) { CHECK_FOR_INTERRUPTS(); /* remember the previous centroid, grab the next one */ prev = curr; curr = &state->centroids[j]; if (curr->mean >= value) break; count += curr->count; } /* * If the value exactly matches the mean of the current centroid, * we're almost there. There may be multiple centroids with the same * mean, so we just need to advance past those. */ if (value == curr->mean) { int64 count_at_value = 0; /* * There may be multiple centroids with this mean (i.e. containing * this value), so find all of them and sum their weights. */ while ((j < state->ncentroids) && (state->centroids[j].mean == value)) { count_at_value += state->centroids[j].count; j++; } result[i] = tdigest_compute_rank(count, count_at_value, state->count); /* the next centroid has a higher mean, so we're done */ continue; } /* * If (value > curr->mean), it means we went through all centroids * without finding one with a larger mean. So the value is above * all centroids, and so it's 1.0 percentile. */ if (value > curr->mean) /* past the largest centroid */ { result[i] = 1; continue; } /* * It's also possible even the first centroid has a higher mean, in * which case the value is 0.0 percentile. */ if (prev == NULL) /* before the smallest centroid */ { result[i] = 0; continue; } /* we have two distinct centroids */ Assert((prev != NULL) && (curr != NULL) && (prev != curr)); Assert(prev->mean < curr->mean); /* * The value lies somewhere between two centroids. We want to figure out * where along the line from the prev node to this node the value is. * * FIXME What if there are multiple centroids with the same mean as the * prev/curr centroid? This probably needs to look up all of them and sum * their counts, just like we did in case of the exact mean equality, no? * Both for the current and previous centroids, so that the approximation * works well. */ d = (curr->mean - prev->mean); /* Use the same exact midpoint arithmetic as the exact-mean branch. */ q1 = tdigest_compute_rank(count - prev->count, prev->count, state->count); q2 = tdigest_compute_rank(count, curr->count, state->count); Assert(q1 <= q2); /* * Calculate the linear interpolation of q1/q2 percentiles. * * We need to be careful about infinity/NaN during calculation. The * means may be so close to +/- DBL_MAX that "d" becomes infinite. * That's equivalent to 0 slope, but we can do a bit better - if this * happens, we halve the values, which makes the difference finite * again (in exchange for loss of precision, but that's acceptable). */ if (isfinite(d)) q = (value - prev->mean) / d; else { /* trick - halve the means, so the difference can't overflow */ q = (value / 2.0 - prev->mean / 2.0) / (curr->mean / 2.0 - prev->mean / 2.0); } result[i] = tdigest_interpolate(q1, q2, q); } } /* * Make sure the aggregate state has space for more centroids. Double the * capacity up to the allowed maximum determined by the compression. */ static void tdigest_aggstate_enlarge(tdigest_aggstate_t *state) { Assert(state->ncentroids <= state->maxcentroids); Assert(state->maxcentroids <= BUFFER_SIZE(state->compression)); /* only enlarge if actually full */ if (state->ncentroids < state->maxcentroids) return; /* double the capacity, but cap it to BUFFER_SIZE */ state->maxcentroids = Min(2 * state->maxcentroids, BUFFER_SIZE(state->compression)); /* repalloc keeps the buffer in its original memory context */ state->centroids = repalloc(state->centroids, state->maxcentroids * sizeof(centroid_t)); /* make sure we have space for the value */ Assert(state->ncentroids < state->maxcentroids); Assert(state->maxcentroids <= BUFFER_SIZE(state->compression)); } /* * Try to reclaim some of the centroid buffer after compaction. * * After compaction, the centroid buffer is guaranteed to be mostly empty, * and there's no guarantee it'll ever be used again. Reclaim large, mostly * unused allocations, keeping headroom for new centroids. * * We only do this for large buffers, so that the chunks are allocted as * separate oversized chunks, which means we actually do free() on them, * instead of stashing them to a freelist. * * And we only do that for buffers that are at least 75% empty, and thus * outside the usual "doubling" strategy. * * This is not free, but it only happens after compaction, which is rather * rare and already fairly expensive, so the additional cost should be * rather acceptable. * * We have to enforce the power-of-2 sizing even during shrinking. */ static void tdigest_aggstate_shrink(tdigest_aggstate_t *state) { /* the buffer has to be at least 75% empty */ if (state->maxcentroids <= 4 * state->ncentroids) return; /* the buffer has to be at least 8KB */ if (state->maxcentroids * sizeof(centroid_t) <= 8192) return; /* * Find the first power-of-2 capacity above ncentroids. It has to be * strictly above, so that the buffer we just resized is not immediately * full again - the next centroid would have to grow it right back. */ state->maxcentroids = BUFFER_INITIAL_SIZE; while (state->maxcentroids <= state->ncentroids) state->maxcentroids *= 2; Assert(state->ncentroids < state->maxcentroids); Assert(state->maxcentroids <= BUFFER_SIZE(state->compression)); state->centroids = repalloc(state->centroids, state->maxcentroids * sizeof(centroid_t)); } /* add a value to the t-digest, trigger a compaction if full */ static void tdigest_add(tdigest_aggstate_t *state, double v) { int compression = state->compression; /* make sure we're not adding bogus NaN/infinity values as centroids */ if (!isfinite(v)) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("all values added to t-digest must be finite"))); /* * If the buffer is full, trigger compaction here so that we have * free space for the new value. */ if (state->ncentroids == BUFFER_SIZE(compression)) tdigest_compact(state); /* ensure there's free space in the aggregate state */ tdigest_aggstate_enlarge(state); /* for a single point, the value is both sum and mean */ state->centroids[state->ncentroids].count = 1; state->centroids[state->ncentroids].mean = v; state->ncentroids++; /* make sure the total does not overflow */ if (pg_add_s64_overflow(state->count, 1, &state->count)) ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("tdigest count overflow"))); } /* * Add a centroid (possibly with count not equal to 1) to the t-digest, * triggering a compaction when the buffer is full. */ static void tdigest_add_centroid(tdigest_aggstate_t *state, double mean, int64 count) { int compression = state->compression; /* we should not have an infinite/NaN mean in a digest */ Assert(isfinite(mean)); /* * If the buffer is full, trigger compaction here so that we have * free space for the new value. */ if (state->ncentroids == BUFFER_SIZE(compression)) tdigest_compact(state); /* ensure there's free space in the aggregate state */ tdigest_aggstate_enlarge(state); state->centroids[state->ncentroids].count = count; state->centroids[state->ncentroids].mean = mean; state->ncentroids++; /* make sure the total does not overflow */ if (pg_add_s64_overflow(state->count, count, &state->count)) ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("tdigest count overflow"))); } /* allocate t-digest with enough space for a requested number of centroids */ static tdigest_t * tdigest_allocate(int ncentroids) { Size len; tdigest_t *digest; char *ptr; len = offsetof(tdigest_t, centroids) + ncentroids * sizeof(centroid_t); ptr = palloc(len); SET_VARSIZE(ptr, len); digest = (tdigest_t *) ptr; digest->flags = 0; digest->ncentroids = 0; digest->count = 0; digest->compression = 0; /* new t-digests automatically store means */ digest->flags |= TDIGEST_STORES_MEAN; return digest; } /* * tdigest_update_format * Update t-digest format to represent centroids as (mean,count). * * Switches the centroids from (sum,count) to (mean,count), so that all * the places processing centroids can use just the new format. * * The caller must own the value before converting a legacy digest, and * the digest is modified in place (without creating a copy). Current * format values are left untouched. */ static void tdigest_update_format(tdigest_t *digest) { int i; /* if already new format, we're done */ if (digest->flags & TDIGEST_STORES_MEAN) return; for (i = 0; i < digest->ncentroids; i++) { CHECK_FOR_INTERRUPTS(); digest->centroids[i].mean = digest->centroids[i].mean / digest->centroids[i].count; } digest->flags |= TDIGEST_STORES_MEAN; } /* * tdigest_prepare * Prepare the digest for additional processing (adding values, ...). * * Detoast, align and normalize a digest, optionally sorting its centroids. * * Digests with the centroids in an arbitrary order are perfectly valid - the * incremental API can leave digests uncompacted (and possibly unsorted) when * compact is false, and the input functions accept such digests too. So the * places walking the centroids in mean order have to do the sort themselves. * * Reuse an owned detoasted/aligned copy for conversion and sorting. Otherwise * copy only before the first modification. The returned value is either the * original datum or one owned allocation, releasable with PG_FREE_IF_COPY. */ static tdigest_t * tdigest_prepare(Datum datum, bool sort) { tdigest_t *digest = tdigest_detoast(datum); int i; /* * If the digest uses the old format, switch to the new one (and make * sure we have a copy, as required by tdigest_update_format). */ if (!(digest->flags & TDIGEST_STORES_MEAN)) { /* if not a copy already, make one */ if ((Pointer) digest == DatumGetPointer(datum)) digest = (tdigest_t *) PG_DETOAST_DATUM_COPY(datum); tdigest_update_format(digest); } /* if not requested to sort centroids, we're done */ if (!sort) return digest; /* if the centroids are already sorted, we're done */ for (i = 1; i < digest->ncentroids; i++) { CHECK_FOR_INTERRUPTS(); /* * XXX Not quite right, it needs to consider the count too, if * the centroids have the same mean (and whether we're below or * above the mean of the whole digest. */ if (digest->centroids[i - 1].mean > digest->centroids[i].mean) break; } /* if the digest is already sorted, bail out */ if (i >= digest->ncentroids) return digest; /* if not a copy already, make one */ if ((Pointer) digest == DatumGetPointer(datum)) digest = (tdigest_t *) PG_DETOAST_DATUM_COPY(datum); tdigest_sort_centroids(digest->centroids, digest->ncentroids, digest->count); return digest; } /* * Allocate a tdigest aggregate state and its query parameters. The centroid * buffer grows separately, leaving the state pointer stable for the executor. * ncentroids is an initial capacity hint, not the number of populated slots. */ static tdigest_aggstate_t * tdigest_aggstate_allocate(int npercentiles, int nvalues, int compression, int ncentroids) { Size len; tdigest_aggstate_t *state; char *ptr; /* at least one of those values is 0 */ Assert(nvalues == 0 || npercentiles == 0); Assert(ncentroids >= 0 && ncentroids <= BUFFER_SIZE(compression)); /* * Allocate a single chunk for the struct, the optional percentile or * hypothetical-value array. */ len = MAXALIGN(sizeof(tdigest_aggstate_t)) + MAXALIGN(sizeof(double) * npercentiles) + MAXALIGN(sizeof(double) * nvalues); ptr = palloc0(len); state = (tdigest_aggstate_t *) ptr; ptr += MAXALIGN(sizeof(tdigest_aggstate_t)); state->nvalues = nvalues; state->npercentiles = npercentiles; state->compression = compression; if (npercentiles > 0) { state->percentiles = (double *) ptr; ptr += MAXALIGN(sizeof(double) * npercentiles); } if (nvalues > 0) { state->values = (double *) ptr; ptr += MAXALIGN(sizeof(double) * nvalues); } Assert(ptr == (char *) state + len); /* * ncentroids is an arbitrary value, but we want to stick to power-of-2 * sizes, to match the size classes used by AllocSet */ state->maxcentroids = BUFFER_INITIAL_SIZE; while (state->maxcentroids < ncentroids) state->maxcentroids *= 2; /* don't use buffers larger than BUFFER_SIZE */ state->maxcentroids = Min(BUFFER_SIZE(compression), state->maxcentroids); state->centroids = palloc(state->maxcentroids * sizeof(centroid_t)); return state; } static void tdigest_aggstate_free(tdigest_aggstate_t *state) { pfree(state->centroids); pfree(state); } static tdigest_t * tdigest_aggstate_to_digest(tdigest_aggstate_t *state, bool compact) { int i; tdigest_t *digest; if (compact) tdigest_compact(state); digest = tdigest_allocate(state->ncentroids); digest->count = state->count; digest->ncentroids = state->ncentroids; digest->compression = state->compression; for (i = 0; i < state->ncentroids; i++) { CHECK_FOR_INTERRUPTS(); digest->centroids[i].mean = state->centroids[i].mean; digest->centroids[i].count = state->centroids[i].count; } return digest; } /* check that the requested percentiles are valid */ static void check_percentiles(const double *percentiles, int npercentiles) { int i; for (i = 0; i < npercentiles; i++) { if (!((percentiles[i] >= 0.0) && (percentiles[i] <= 1.0))) elog(ERROR, "invalid percentile value %f, should be in [0.0, 1.0]", percentiles[i]); } } static void check_compression(int compression) { if (compression < MIN_COMPRESSION || compression > MAX_COMPRESSION) elog(ERROR, "invalid compression value %d", compression); } static void check_trim_values(double low, double high) { if (!((low >= 0.0) && (low <= 1.0))) elog(ERROR, "invalid low percentile value %f, should be in [0.0, 1.0]", low); if (!((high >= 0.0) && (high <= 1.0))) elog(ERROR, "invalid high percentile value %f, should be in [0.0, 1.0]", high); if (low > high) elog(ERROR, "invalid low/high percentile values %f/%f, should be low <= high", low, high); } /* * Add an input value to the aggregate state, creating it if needed. * Shared by tdigest and tdigest_percentile with a single percentile. * * Transition functions capture compression and query parameters when the * first non-NULL input creates the state. Callers should keep those * parameters constant within a group; later rows are not checked for changes. */ Datum tdigest_add_double(PG_FUNCTION_ARGS) { tdigest_aggstate_t *state; MemoryContext aggcontext; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_add_double called in non-aggregate context"); /* * We want to skip NULL values altogether - we return either the existing * t-digest (if it already exists) or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* if there's no digest allocated, create it now */ if (PG_ARGISNULL(0)) { int compression; double *percentiles = NULL; int npercentiles = 0; MemoryContext oldcontext; if (PG_ARGISNULL(2)) elog(ERROR, "compression must not be NULL"); compression = PG_GETARG_INT32(2); check_compression(compression); oldcontext = MemoryContextSwitchTo(aggcontext); if (PG_NARGS() >= 4) { /* the percentile is required to create the aggregate state */ if (PG_ARGISNULL(3)) elog(ERROR, "percentile must not be NULL"); percentiles = (double *) palloc(sizeof(double)); percentiles[0] = PG_GETARG_FLOAT8(3); npercentiles = 1; check_percentiles(percentiles, npercentiles); } state = tdigest_aggstate_allocate(npercentiles, 0, compression, 0); if (percentiles) { memcpy(state->percentiles, percentiles, sizeof(double) * npercentiles); pfree(percentiles); } MemoryContextSwitchTo(oldcontext); } else state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); tdigest_add(state, PG_GETARG_FLOAT8(1)); AssertCheckTDigestAggState(state); PG_RETURN_POINTER(state); } /* * Add a value with a given count to the t-digest, as a sequence of properly * sized centroids. * * This is an alternative to adding a single centroid, representing all the * points with the same value. It follows all the rules on centroid sizes, * etc. * * The centroids are handed over to the aggregate state as they are computed, * instead of building a t-digest first. The number of centroids the loop * produces is not bounded by the compression (on the tails the calculated * size drops below 1, and gets clamped), so there is no size of a centroid * array that would be guaranteed to be sufficient. */ static void tdigest_add_generated(tdigest_aggstate_t *state, double value, int64 count) { int64 count_so_far; int64 count_remaining; double denom; double normalizer; int compression = state->compression; /* make sure we're not adding bogus NaN/infinity values as centroids */ if (!isfinite(value)) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("all values added to t-digest must be finite"))); denom = 2 * M_PI * count * log(count); normalizer = compression / denom; count_so_far = 0; /* does not include current centroid */ count_remaining = count; /* * Create largest possible centroids, until we run out of items. In each * step we need to find the largest possible well-formed centroid, i.e. one * that matches the two conditions: * * z <= q0 * (1 - q0) where q0 = (count_so_far / count) * * z <= q2 * (1 - q2) where q2 = (count_so_far + X) / count; * * with z = (X * normalizer). X being the value we need to determine. Solving * q0 is trivial, while q2 leads to a quadratic equation with two roots. */ while (count_remaining > 0) { int64 proposed_count; double q0; double b, c, d; double r1, r2; CHECK_FOR_INTERRUPTS(); /* * Solving z <= q0 * (1 - q0) is trivial. * * Just like in tdigest_compact, we must not calculate (1 - q0) by * subtracting the two doubles - for q0 close to 1 that cancels out all * the significant digits. We already have count_remaining, which is * exactly (count - count_so_far), so use that remainder instead. */ q0 = count_so_far / (double) count; r1 = (q0 * (count_remaining / (double) count) / normalizer); /* * Solve z <= q2 * (1 - q2) as a quadratic equation. The inequality we * need to solve is * * 0 <= a * x^2 + b * x + c * * with (a = -1) and the following coefficients. * * XXX The counts may be very high values (int64), so we need to be * careful to prevent overflows by doing everything with double. * * XXX c is mathematically (count_so_far * (count - count_so_far)), so * calculate it as a plain product of the two exact integers. The * expanded form (count_so_far * count - count_so_far * count_so_far) * is a difference of two huge and nearly equal values, which loses * almost all the precision. */ b = ((double) count - 2 * (double) count_so_far - (double) count * (double) count * normalizer); c = ((double) count_so_far * (double) count_remaining); /* * As this is an "upside down" parabola, the values between the roots * are positive - we're looking for the larger of the two roots, which * for a = -1 is (b + sqrt(b*b + 4*c)) / 2. * * XXX Evaluating that expression directly is only safe for b >= 0. * For b < 0 the sqrt is very close to -b, so the addition cancels out * all the significant digits (and often yields exactly zero, forcing * us to emit a single-item centroid). Use the equivalent "conjugate" * form 2*c / (sqrt(b*b + 4*c) - b) in that case, which only ever adds * values of the same sign. Both branches are hit in practice - b is * positive whenever compression < 2*pi*ln(count). * * XXX c is never negative, so the discriminant is a sum of two * non-negative values and the sqrt is always well defined. */ d = sqrt(b * b + 4 * c); if (b >= 0) r2 = (b + d) / 2; else r2 = (2 * c) / (d - b); /* * paranoia: We should not be dealing with NaN values here. Crash in * debug build, double_to_int64 will mitigate it in regular builds. */ Assert(isfinite(r1) && isfinite(r2)); /* * We need to meet both conditions, so use the smaller solution. The * value may be large (or NaN), so clamp it - we must not add more * than what remains anyway. */ proposed_count = double_to_int64(floor(Min(r1, r2)), count_remaining); /* * It's possible to get very low values on the tails, but we must add * at least something, otherwise we'd get infinite loops. */ proposed_count = Max(proposed_count, 1); tdigest_add_centroid(state, value, proposed_count); count_so_far += proposed_count; count_remaining -= proposed_count; } } /* * Add an input value with a count to the aggregate state, creating it if * needed. Shared by tdigest and tdigest_percentile with a single percentile. */ Datum tdigest_add_double_count(PG_FUNCTION_ARGS) { int64 i; int64 count; tdigest_aggstate_t *state; MemoryContext aggcontext; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_add_double_count called in non-aggregate context"); /* * We want to skip NULL values altogether - we return either the existing * t-digest (if it already exists) or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* if there's no digest allocated, create it now */ if (PG_ARGISNULL(0)) { int compression; double *percentiles = NULL; int npercentiles = 0; MemoryContext oldcontext; if (PG_ARGISNULL(3)) elog(ERROR, "compression must not be NULL"); compression = PG_GETARG_INT32(3); check_compression(compression); oldcontext = MemoryContextSwitchTo(aggcontext); if (PG_NARGS() >= 5) { /* the percentile is required to create the aggregate state */ if (PG_ARGISNULL(4)) elog(ERROR, "percentile must not be NULL"); percentiles = (double *) palloc(sizeof(double)); percentiles[0] = PG_GETARG_FLOAT8(4); npercentiles = 1; check_percentiles(percentiles, npercentiles); } state = tdigest_aggstate_allocate(npercentiles, 0, compression, 0); if (percentiles) { memcpy(state->percentiles, percentiles, sizeof(double) * npercentiles); pfree(percentiles); } MemoryContextSwitchTo(oldcontext); } else state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); if (PG_ARGISNULL(2)) { count = 1; } else count = PG_GETARG_INT64(2); /* can't add values with non-positive counts */ if (count <= 0) elog(ERROR, "invalid count value %lld, must be a positive value", (long long) count); /* * When adding more values than would fit into an empty buffer (and * thus likely causing too many compactions), we instead add them as * properly sized centroids. * * This is much faster, because the centroids can be generated in one go, * so there are only very few compactions. */ if (count > BUFFER_SIZE(state->compression)) { tdigest_add_generated(state, PG_GETARG_FLOAT8(1), count); count = 0; } /* * If there are only a few values, just add them one by one, so that * we do proper compaction and sizing of centroids. Otherwise we might end * up with oversized centroids on the tails etc. */ for (i = 0; i < count; i++) tdigest_add(state, PG_GETARG_FLOAT8(1)); AssertCheckTDigestAggState(state); PG_RETURN_POINTER(state); } /* * Add an input value to the aggregate state for tdigest_percentile_of * with a single hypothetical value. */ Datum tdigest_add_double_values(PG_FUNCTION_ARGS) { tdigest_aggstate_t *state; MemoryContext aggcontext; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_add_double_values called in non-aggregate context"); /* * We want to skip NULL values altogether - we return either the existing * t-digest (if it already exists) or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* if there's no digest allocated, create it now */ if (PG_ARGISNULL(0)) { int compression; double *values = NULL; int nvalues = 0; MemoryContext oldcontext; if (PG_ARGISNULL(2)) elog(ERROR, "compression must not be NULL"); compression = PG_GETARG_INT32(2); check_compression(compression); oldcontext = MemoryContextSwitchTo(aggcontext); if (PG_NARGS() >= 4) { /* the value is required to create the aggregate state */ if (PG_ARGISNULL(3)) elog(ERROR, "value must not be NULL"); values = (double *) palloc(sizeof(double)); values[0] = PG_GETARG_FLOAT8(3); nvalues = 1; } state = tdigest_aggstate_allocate(0, nvalues, compression, 0); if (values) { memcpy(state->values, values, sizeof(double) * nvalues); pfree(values); } MemoryContextSwitchTo(oldcontext); } else state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); tdigest_add(state, PG_GETARG_FLOAT8(1)); AssertCheckTDigestAggState(state); PG_RETURN_POINTER(state); } /* * Add an input value with a count to the aggregate state for * tdigest_percentile_of with a single hypothetical value. */ Datum tdigest_add_double_values_count(PG_FUNCTION_ARGS) { int64 i; int64 count; tdigest_aggstate_t *state; MemoryContext aggcontext; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_add_double_values_count called in non-aggregate context"); /* * We want to skip NULL values altogether - we return either the existing * t-digest (if it already exists) or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* if there's no digest allocated, create it now */ if (PG_ARGISNULL(0)) { int compression; double *values = NULL; int nvalues = 0; MemoryContext oldcontext; if (PG_ARGISNULL(3)) elog(ERROR, "compression must not be NULL"); compression = PG_GETARG_INT32(3); check_compression(compression); oldcontext = MemoryContextSwitchTo(aggcontext); if (PG_NARGS() >= 5) { /* the value is required to create the aggregate state */ if (PG_ARGISNULL(4)) elog(ERROR, "value must not be NULL"); values = (double *) palloc(sizeof(double)); values[0] = PG_GETARG_FLOAT8(4); nvalues = 1; } state = tdigest_aggstate_allocate(0, nvalues, compression, 0); if (values) { memcpy(state->values, values, sizeof(double) * nvalues); pfree(values); } MemoryContextSwitchTo(oldcontext); } else state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); if (PG_ARGISNULL(2)) { count = 1; } else count = PG_GETARG_INT64(2); /* can't add values with non-positive counts */ if (count <= 0) elog(ERROR, "invalid count value %lld, must be a positive value", (long long) count); /* * When adding more values than would fit into an empty buffer (and * thus likely causing too many compactions), we instead add them as * properly sized centroids. * * This is much faster, because the centroids can be generated in one go, * so there are only very few compactions. */ if (count > BUFFER_SIZE(state->compression)) { tdigest_add_generated(state, PG_GETARG_FLOAT8(1), count); count = 0; } /* * If there are only a few values, just add them one by one, so that * we do proper compaction and sizing of centroids. Otherwise we might end * up with oversized centroids on the tails etc. */ for (i = 0; i < count; i++) tdigest_add(state, PG_GETARG_FLOAT8(1)); AssertCheckTDigestAggState(state); PG_RETURN_POINTER(state); } /* * Merge an input digest into the aggregate state, creating it if needed. * Shared by tdigest and tdigest_percentile with a single percentile. */ Datum tdigest_add_digest(PG_FUNCTION_ARGS) { int i; tdigest_aggstate_t *state; tdigest_t *digest; MemoryContext aggcontext; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_add_digest called in non-aggregate context"); /* * We want to skip NULL values altogether - we return either the existing * t-digest (if it already exists) or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* make sure we have a detoasted copy */ digest = tdigest_prepare(PG_GETARG_DATUM(1), false); /* if there's no aggregate state allocated, create it now */ if (PG_ARGISNULL(0)) { double *percentiles = NULL; int npercentiles = 0; MemoryContext oldcontext; oldcontext = MemoryContextSwitchTo(aggcontext); if (PG_NARGS() >= 3) { /* the percentile is required to create the aggregate state */ if (PG_ARGISNULL(2)) elog(ERROR, "percentile must not be NULL"); percentiles = (double *) palloc(sizeof(double)); percentiles[0] = PG_GETARG_FLOAT8(2); npercentiles = 1; check_percentiles(percentiles, npercentiles); } state = tdigest_aggstate_allocate(npercentiles, 0, digest->compression, digest->ncentroids); if (percentiles) { memcpy(state->percentiles, percentiles, sizeof(double) * npercentiles); pfree(percentiles); } MemoryContextSwitchTo(oldcontext); } else state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); /* * Keep the compression chosen when the state was created, even if the * input uses a different setting. Compaction can merge incoming * centroids, but cannot split them to recover lost detail. */ /* copy data from the tdigest into the aggstate */ for (i = 0; i < digest->ncentroids; i++) { CHECK_FOR_INTERRUPTS(); tdigest_add_centroid(state, digest->centroids[i].mean, digest->centroids[i].count); } AssertCheckTDigestAggState(state); PG_FREE_IF_COPY(digest, 1); PG_RETURN_POINTER(state); } /* * Merge an input digest into the aggregate state for tdigest_percentile_of * with a single hypothetical value. */ Datum tdigest_add_digest_values(PG_FUNCTION_ARGS) { int i; tdigest_aggstate_t *state; tdigest_t *digest; MemoryContext aggcontext; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_add_digest_values called in non-aggregate context"); /* * We want to skip NULL values altogether - we return either the existing * t-digest (if it already exists) or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* make sure we have a detoasted copy */ digest = tdigest_prepare(PG_GETARG_DATUM(1), false); /* if there's no aggregate state allocated, create it now */ if (PG_ARGISNULL(0)) { double *values = NULL; int nvalues = 0; MemoryContext oldcontext; oldcontext = MemoryContextSwitchTo(aggcontext); if (PG_NARGS() >= 3) { /* the value is required to create the aggregate state */ if (PG_ARGISNULL(2)) elog(ERROR, "value must not be NULL"); values = (double *) palloc(sizeof(double)); values[0] = PG_GETARG_FLOAT8(2); nvalues = 1; } state = tdigest_aggstate_allocate(0, nvalues, digest->compression, digest->ncentroids); if (values) { memcpy(state->values, values, sizeof(double) * nvalues); pfree(values); } MemoryContextSwitchTo(oldcontext); } else state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); /* * Keep the compression chosen when the state was created, even if the * input uses a different setting. Compaction can merge incoming * centroids, but cannot split them to recover lost detail. */ for (i = 0; i < digest->ncentroids; i++) { CHECK_FOR_INTERRUPTS(); tdigest_add_centroid(state, digest->centroids[i].mean, digest->centroids[i].count); } AssertCheckTDigestAggState(state); PG_FREE_IF_COPY(digest, 1); PG_RETURN_POINTER(state); } /* * Add an input value to the aggregate state for tdigest_percentile with * an array of requested percentiles. */ Datum tdigest_add_double_array(PG_FUNCTION_ARGS) { tdigest_aggstate_t *state; MemoryContext aggcontext; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_add_double_array called in non-aggregate context"); /* * We want to skip NULL values altogether - we return either the existing * t-digest or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* if there's no digest allocated, create it now */ if (PG_ARGISNULL(0)) { int compression; const double *percentiles; int npercentiles; MemoryContext oldcontext; ArrayType *array; if (PG_ARGISNULL(2)) elog(ERROR, "compression must not be NULL"); compression = PG_GETARG_INT32(2); check_compression(compression); /* Percentiles are required in order to create the aggregate state. */ if (PG_ARGISNULL(3)) elog(ERROR, "percentiles must not be NULL"); /* process the input array in caller's memory context */ array = PG_GETARG_ARRAYTYPE_P(3); percentiles = array_to_double(array, "a percentile value", &npercentiles); check_percentiles(percentiles, npercentiles); oldcontext = MemoryContextSwitchTo(aggcontext); state = tdigest_aggstate_allocate(npercentiles, 0, compression, 0); memcpy(state->percentiles, percentiles, sizeof(double) * npercentiles); MemoryContextSwitchTo(oldcontext); /* free the parsed data, possibly detoasted copy of the array */ PG_FREE_IF_COPY(array, 3); } else state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); tdigest_add(state, PG_GETARG_FLOAT8(1)); AssertCheckTDigestAggState(state); PG_RETURN_POINTER(state); } /* * Add an input value with a count to the aggregate state for * tdigest_percentile with an array of requested percentiles. */ Datum tdigest_add_double_array_count(PG_FUNCTION_ARGS) { int64 i; int64 count; tdigest_aggstate_t *state; MemoryContext aggcontext; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_add_double_array_count called in non-aggregate context"); /* * We want to skip NULL values altogether - we return either the existing * t-digest or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* if there's no digest allocated, create it now */ if (PG_ARGISNULL(0)) { int compression; const double *percentiles; int npercentiles; MemoryContext oldcontext; ArrayType *array; if (PG_ARGISNULL(3)) elog(ERROR, "compression must not be NULL"); compression = PG_GETARG_INT32(3); check_compression(compression); /* Percentiles are required in order to create the aggregate state. */ if (PG_ARGISNULL(4)) elog(ERROR, "percentiles must not be NULL"); /* process the input array in caller's memory context */ array = PG_GETARG_ARRAYTYPE_P(4); percentiles = array_to_double(array, "a percentile value", &npercentiles); check_percentiles(percentiles, npercentiles); oldcontext = MemoryContextSwitchTo(aggcontext); state = tdigest_aggstate_allocate(npercentiles, 0, compression, 0); memcpy(state->percentiles, percentiles, sizeof(double) * npercentiles); MemoryContextSwitchTo(oldcontext); /* free the parsed data, possibly detoasted copy of the array */ PG_FREE_IF_COPY(array, 4); } else state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); if (PG_ARGISNULL(2)) { count = 1; } else count = PG_GETARG_INT64(2); /* can't add values with non-positive counts */ if (count <= 0) elog(ERROR, "invalid count value %lld, must be a positive value", (long long) count); /* * When adding more values than would fit into an empty buffer (and * thus likely causing too many compactions), we instead add them as * properly sized centroids. * * This is much faster, because the centroids can be generated in one go, * so there are only very few compactions. */ if (count > BUFFER_SIZE(state->compression)) { tdigest_add_generated(state, PG_GETARG_FLOAT8(1), count); count = 0; } /* * Add the values one by one, not as one large centroid with the count. * We do it like this to allow proper compaction and sizing of centroids, * otherwise we might end up with oversized centroids on the tails etc. * * XXX If this turns out a bit too expensive, we may try determining the * size by looking for the smallest centroid covering this value. */ for (i = 0; i < count; i++) tdigest_add(state, PG_GETARG_FLOAT8(1)); AssertCheckTDigestAggState(state); PG_RETURN_POINTER(state); } /* * Add an input value to the aggregate state for tdigest_percentile_of * with an array of hypothetical values. */ Datum tdigest_add_double_array_values(PG_FUNCTION_ARGS) { tdigest_aggstate_t *state; MemoryContext aggcontext; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_add_double_array called in non-aggregate context"); /* * We want to skip NULL values altogether - we return either the existing * t-digest or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* if there's no digest allocated, create it now */ if (PG_ARGISNULL(0)) { int compression; const double *values; int nvalues; MemoryContext oldcontext; ArrayType *array; if (PG_ARGISNULL(2)) elog(ERROR, "compression must not be NULL"); compression = PG_GETARG_INT32(2); check_compression(compression); /* Values are required in order to create the aggregate state. */ if (PG_ARGISNULL(3)) elog(ERROR, "values must not be NULL"); /* process the input array in caller's memory context */ array = PG_GETARG_ARRAYTYPE_P(3); values = array_to_double(array, "a value", &nvalues); oldcontext = MemoryContextSwitchTo(aggcontext); state = tdigest_aggstate_allocate(0, nvalues, compression, 0); memcpy(state->values, values, sizeof(double) * nvalues); MemoryContextSwitchTo(oldcontext); /* free the parsed data, possibly detoasted copy of the array */ PG_FREE_IF_COPY(array, 3); } else state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); tdigest_add(state, PG_GETARG_FLOAT8(1)); AssertCheckTDigestAggState(state); PG_RETURN_POINTER(state); } /* * Add an input value with a count to the aggregate state for * tdigest_percentile_of with an array of hypothetical values. */ Datum tdigest_add_double_array_values_count(PG_FUNCTION_ARGS) { int64 i; int64 count; tdigest_aggstate_t *state; MemoryContext aggcontext; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_add_double_array_values_count called in non-aggregate context"); /* * We want to skip NULL values altogether - we return either the existing * t-digest or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* if there's no digest allocated, create it now */ if (PG_ARGISNULL(0)) { int compression; const double *values; int nvalues; MemoryContext oldcontext; ArrayType *array; if (PG_ARGISNULL(3)) elog(ERROR, "compression must not be NULL"); compression = PG_GETARG_INT32(3); check_compression(compression); /* Values are required in order to create the aggregate state. */ if (PG_ARGISNULL(4)) elog(ERROR, "values must not be NULL"); /* process the input array in caller's memory context */ array = PG_GETARG_ARRAYTYPE_P(4); values = array_to_double(array, "a value", &nvalues); oldcontext = MemoryContextSwitchTo(aggcontext); state = tdigest_aggstate_allocate(0, nvalues, compression, 0); memcpy(state->values, values, sizeof(double) * nvalues); MemoryContextSwitchTo(oldcontext); /* free the parsed data, possibly detoasted copy of the array */ PG_FREE_IF_COPY(array, 4); } else state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); if (PG_ARGISNULL(2)) { count = 1; } else count = PG_GETARG_INT64(2); /* can't add values with non-positive counts */ if (count <= 0) elog(ERROR, "invalid count value %lld, must be a positive value", (long long) count); /* * When adding more values than would fit into an empty buffer (and * thus likely causing too many compactions), we instead add them as * properly sized centroids. * * This is much faster, because the centroids can be generated in one go, * so there are only very few compactions. */ if (count > BUFFER_SIZE(state->compression)) { tdigest_add_generated(state, PG_GETARG_FLOAT8(1), count); count = 0; } /* * Add the values one by one, not as one large centroid with the count. * We do it like this to allow proper compaction and sizing of centroids, * otherwise we might end up with oversized centroids on the tails etc. * * XXX If this turns out a bit too expensive, we may try determining the * size by looking for the smallest centroid covering this value. */ for (i = 0; i < count; i++) tdigest_add(state, PG_GETARG_FLOAT8(1)); AssertCheckTDigestAggState(state); PG_RETURN_POINTER(state); } /* * Merge an input digest into the aggregate state for tdigest_percentile * with an array of requested percentiles. */ Datum tdigest_add_digest_array(PG_FUNCTION_ARGS) { int i; tdigest_aggstate_t *state; tdigest_t *digest; MemoryContext aggcontext; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_add_digest_array called in non-aggregate context"); /* * We want to skip NULL values altogether - we return either the existing * t-digest (if it already exists) or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* make sure we have a detoasted copy */ digest = tdigest_prepare(PG_GETARG_DATUM(1), false); /* if there's no aggregate state allocated, create it now */ if (PG_ARGISNULL(0)) { const double *percentiles; int npercentiles; MemoryContext oldcontext; ArrayType *array; /* Percentiles are required in order to create the aggregate state. */ if (PG_ARGISNULL(2)) elog(ERROR, "percentiles must not be NULL"); /* process the input array in caller's memory context */ array = PG_GETARG_ARRAYTYPE_P(2); percentiles = array_to_double(array, "a percentile value", &npercentiles); check_percentiles(percentiles, npercentiles); oldcontext = MemoryContextSwitchTo(aggcontext); state = tdigest_aggstate_allocate(npercentiles, 0, digest->compression, digest->ncentroids); memcpy(state->percentiles, percentiles, sizeof(double) * npercentiles); MemoryContextSwitchTo(oldcontext); /* free the parsed data, possibly detoasted copy of the array */ PG_FREE_IF_COPY(array, 2); } else state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); /* * Keep the compression chosen when the state was created, even if the * input uses a different setting. Compaction can merge incoming * centroids, but cannot split them to recover lost detail. */ for (i = 0; i < digest->ncentroids; i++) { CHECK_FOR_INTERRUPTS(); tdigest_add_centroid(state, digest->centroids[i].mean, digest->centroids[i].count); } AssertCheckTDigestAggState(state); PG_FREE_IF_COPY(digest, 1); PG_RETURN_POINTER(state); } /* * Merge an input digest into the aggregate state for tdigest_percentile_of * with an array of hypothetical values. */ Datum tdigest_add_digest_array_values(PG_FUNCTION_ARGS) { int i; tdigest_aggstate_t *state; tdigest_t *digest; MemoryContext aggcontext; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_add_digest_array_values called in non-aggregate context"); /* * We want to skip NULL values altogether - we return either the existing * t-digest (if it already exists) or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* make sure we have a detoasted copy */ digest = tdigest_prepare(PG_GETARG_DATUM(1), false); /* if there's no aggregate state allocated, create it now */ if (PG_ARGISNULL(0)) { const double *values; int nvalues; MemoryContext oldcontext; ArrayType *array; /* Values are required in order to create the aggregate state. */ if (PG_ARGISNULL(2)) elog(ERROR, "values must not be NULL"); /* process the input array in caller's memory context */ array = PG_GETARG_ARRAYTYPE_P(2); values = array_to_double(array, "a value", &nvalues); oldcontext = MemoryContextSwitchTo(aggcontext); state = tdigest_aggstate_allocate(0, nvalues, digest->compression, digest->ncentroids); memcpy(state->values, values, sizeof(double) * nvalues); MemoryContextSwitchTo(oldcontext); /* free the parsed data, possibly detoasted copy of the array */ PG_FREE_IF_COPY(array, 2); } else state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); /* * Keep the compression chosen when the state was created, even if the * input uses a different setting. Compaction can merge incoming * centroids, but cannot split them to recover lost detail. */ for (i = 0; i < digest->ncentroids; i++) { CHECK_FOR_INTERRUPTS(); tdigest_add_centroid(state, digest->centroids[i].mean, digest->centroids[i].count); } AssertCheckTDigestAggState(state); PG_FREE_IF_COPY(digest, 1); PG_RETURN_POINTER(state); } /* * Compute a percentile from the aggregate state. Final function for * tdigest_percentile with a single requested percentile. */ Datum tdigest_percentiles(PG_FUNCTION_ARGS) { tdigest_aggstate_t *state; MemoryContext aggcontext; double ret; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_percentiles called in non-aggregate context"); /* if there's no digest, return NULL */ if (PG_ARGISNULL(0)) PG_RETURN_NULL(); state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); /* windows and shared aggregates may need the original state again. */ if (AggStateIsShared(fcinfo)) { tdigest_aggstate_t *copy = tdigest_copy(state); tdigest_compute_quantiles(copy, &ret); tdigest_aggstate_free(copy); } else tdigest_compute_quantiles(state, &ret); PG_RETURN_FLOAT8(ret); } /* * Compute an inverse percentile from the aggregate state. Final function * for tdigest_percentile_of with a single hypothetical value. */ Datum tdigest_percentiles_of(PG_FUNCTION_ARGS) { tdigest_aggstate_t *state; MemoryContext aggcontext; double ret; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_percentiles_of called in non-aggregate context"); /* if there's no digest, return NULL */ if (PG_ARGISNULL(0)) PG_RETURN_NULL(); state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); /* windows and shared aggregates may need the original state again. */ if (AggStateIsShared(fcinfo)) { tdigest_aggstate_t *copy = tdigest_copy(state); tdigest_compute_quantiles_of(copy, &ret); tdigest_aggstate_free(copy); } else tdigest_compute_quantiles_of(state, &ret); PG_RETURN_FLOAT8(ret); } /* * Build a t-digest varlena value from the aggregate state. */ Datum tdigest_digest(PG_FUNCTION_ARGS) { tdigest_t *digest; tdigest_aggstate_t *state; MemoryContext aggcontext; bool shared; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_digest called in non-aggregate context"); /* if there's no digest, return NULL */ if (PG_ARGISNULL(0)) PG_RETURN_NULL(); state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); /* windows and shared aggregates may need the original state again. */ shared = AggStateIsShared(fcinfo); if (shared) state = tdigest_copy(state); digest = tdigest_aggstate_to_digest(state, true); if (shared) tdigest_aggstate_free(state); PG_RETURN_POINTER(digest); } /* * Compute percentiles from the aggregate state. Final function for * tdigest_percentile with an array of requested percentiles. */ Datum tdigest_array_percentiles(PG_FUNCTION_ARGS) { double *result; ArrayType *array; MemoryContext aggcontext; tdigest_aggstate_t *state; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_array_percentiles called in non-aggregate context"); if (PG_ARGISNULL(0)) PG_RETURN_NULL(); state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); /* write results directly into a correctly sized ArrayType */ array = double_array_allocate(state->npercentiles); result = (double *) ARR_DATA_PTR(array); /* windows and shared aggregates may need the original state again. */ if (AggStateIsShared(fcinfo)) { tdigest_aggstate_t *copy = tdigest_copy(state); tdigest_compute_quantiles(copy, result); tdigest_aggstate_free(copy); } else tdigest_compute_quantiles(state, result); PG_RETURN_ARRAYTYPE_P(array); } /* * Compute inverse percentiles from the aggregate state. Final function * for tdigest_percentile_of with an array of hypothetical values. */ Datum tdigest_array_percentiles_of(PG_FUNCTION_ARGS) { double *result; ArrayType *array; MemoryContext aggcontext; tdigest_aggstate_t *state; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_array_percentiles_of called in non-aggregate context"); if (PG_ARGISNULL(0)) PG_RETURN_NULL(); state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); /* write results directly into a correctly sized ArrayType */ array = double_array_allocate(state->nvalues); result = (double *) ARR_DATA_PTR(array); /* windows and shared aggregates may need the original state again. */ if (AggStateIsShared(fcinfo)) { tdigest_aggstate_t *copy = tdigest_copy(state); tdigest_compute_quantiles_of(copy, result); tdigest_aggstate_free(copy); } else tdigest_compute_quantiles_of(state, result); PG_RETURN_ARRAYTYPE_P(array); } Datum tdigest_serial(PG_FUNCTION_ARGS) { bytea *v; tdigest_aggstate_t *state; Size len; char *ptr; state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); len = offsetof(tdigest_aggstate_t, percentiles) + state->npercentiles * sizeof(double) + state->nvalues * sizeof(double) + state->ncentroids * sizeof(centroid_t); v = palloc(len + VARHDRSZ); SET_VARSIZE(v, len + VARHDRSZ); ptr = VARDATA(v); memcpy(ptr, state, offsetof(tdigest_aggstate_t, percentiles)); ptr += offsetof(tdigest_aggstate_t, percentiles); if (state->npercentiles > 0) { memcpy(ptr, state->percentiles, sizeof(double) * state->npercentiles); ptr += sizeof(double) * state->npercentiles; } if (state->nvalues > 0) { memcpy(ptr, state->values, sizeof(double) * state->nvalues); ptr += sizeof(double) * state->nvalues; } memcpy(ptr, state->centroids, sizeof(centroid_t) * state->ncentroids); ptr += sizeof(centroid_t) * state->ncentroids; Assert(VARDATA(v) + len == ptr); PG_RETURN_POINTER(v); } /* * XXX Unlike the other "input" functions (tdigest_in/tdigest_recv), this * does not validate the digest at all. We assume this function is used only * on data we created in the same query (possibly in a parallel worker), * and not on untrusted values controlled by the user (which is why the other * input functions need the validation). */ Datum tdigest_deserial(PG_FUNCTION_ARGS) { bytea *v = (bytea *) PG_GETARG_POINTER(0); char *ptr = VARDATA_ANY(v); tdigest_aggstate_t tmp; tdigest_aggstate_t *state; /* copy aggstate header into a local variable */ memcpy(&tmp, ptr, offsetof(tdigest_aggstate_t, percentiles)); ptr += offsetof(tdigest_aggstate_t, percentiles); state = tdigest_aggstate_allocate(tmp.npercentiles, tmp.nvalues, tmp.compression, tmp.ncentroids); if (tmp.npercentiles > 0) { memcpy(state->percentiles, ptr, tmp.npercentiles * sizeof(double)); ptr += tmp.npercentiles * sizeof(double); } if (tmp.nvalues > 0) { memcpy(state->values, ptr, tmp.nvalues * sizeof(double)); ptr += tmp.nvalues * sizeof(double); } /* * copy the aggstate header (except for the maxcentroids, which would * corrupt the aggstate we just allocated) */ state->count = tmp.count; state->ncompactions = tmp.ncompactions; state->compression = tmp.compression; /* skip maxcentroids */ state->ncentroids = tmp.ncentroids; state->ncompacted = tmp.ncompacted; state->npercentiles = tmp.npercentiles; state->nvalues = tmp.nvalues; state->trim_low = tmp.trim_low; state->trim_high = tmp.trim_high; /* we don't need to move the pointer */ /* copy the centroids back */ memcpy(state->centroids, ptr, sizeof(centroid_t) * state->ncentroids); ptr += sizeof(centroid_t) * state->ncentroids; Assert(ptr == VARDATA_ANY(v) + VARSIZE_ANY_EXHDR(v)); PG_RETURN_POINTER(state); } static tdigest_aggstate_t * tdigest_copy(tdigest_aggstate_t *state) { tdigest_aggstate_t *copy; copy = tdigest_aggstate_allocate(state->npercentiles, state->nvalues, state->compression, state->ncentroids); /* * copy the aggstate header (except for the maxcentroids, which would * corrupt the aggstate we just allocated) */ copy->count = state->count; copy->ncompactions = state->ncompactions; copy->compression = state->compression; /* skip maxcentroids */ copy->ncentroids = state->ncentroids; copy->ncompacted = state->ncompacted; copy->npercentiles = state->npercentiles; copy->nvalues = state->nvalues; copy->trim_low = state->trim_low; copy->trim_high = state->trim_high; if (state->nvalues > 0) memcpy(copy->values, state->values, sizeof(double) * state->nvalues); if (state->npercentiles > 0) memcpy(copy->percentiles, state->percentiles, sizeof(double) * state->npercentiles); memcpy(copy->centroids, state->centroids, state->ncentroids * sizeof(centroid_t)); return copy; } Datum tdigest_combine(PG_FUNCTION_ARGS) { tdigest_aggstate_t *src; tdigest_aggstate_t *dst; MemoryContext aggcontext; MemoryContext oldcontext; int i; if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_combine called in non-aggregate context"); /* if no "merged" state yet, try creating it */ if (PG_ARGISNULL(0)) { /* nope, the second argument is NULL too, so return NULL */ if (PG_ARGISNULL(1)) PG_RETURN_NULL(); /* the second argument is not NULL, so copy it */ src = (tdigest_aggstate_t *) PG_GETARG_POINTER(1); /* copy the digest into the right long-lived memory context */ oldcontext = MemoryContextSwitchTo(aggcontext); src = tdigest_copy(src); MemoryContextSwitchTo(oldcontext); PG_RETURN_POINTER(src); } /* * If the second argument is NULL, just return the first one (we know * it's not NULL at this point). */ if (PG_ARGISNULL(1)) PG_RETURN_DATUM(PG_GETARG_DATUM(0)); /* Now we know neither argument is NULL, so merge them. */ src = (tdigest_aggstate_t *) PG_GETARG_POINTER(1); dst = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); AssertCheckTDigestAggState(dst); AssertCheckTDigestAggState(src); /* * Keep the compression chosen when the state was created, even if the * input uses a different setting. Compaction can merge incoming * centroids, but cannot split them to recover lost detail. */ /* copy data from the tdigest into the aggstate */ for (i = 0; i < src->ncentroids; i++) { CHECK_FOR_INTERRUPTS(); tdigest_add_centroid(dst, src->centroids[i].mean, src->centroids[i].count); } AssertCheckTDigestAggState(dst); PG_RETURN_POINTER(dst); } /* API for incremental updates */ /* * expand the t-digest into an in-memory aggregate state */ static tdigest_aggstate_t * tdigest_digest_to_aggstate(tdigest_t *digest) { int i; tdigest_aggstate_t *state; /* * The digest should have gone through tdigest_prepare(), which means * we should have a local copy in the new format. */ Assert(digest->flags & TDIGEST_STORES_MEAN); state = tdigest_aggstate_allocate(0, 0, digest->compression, digest->ncentroids); /* copy data from the tdigest into the aggstate */ for (i = 0; i < digest->ncentroids; i++) { CHECK_FOR_INTERRUPTS(); tdigest_add_centroid(state, digest->centroids[i].mean, digest->centroids[i].count); } AssertCheckTDigestAggState(state); return state; } /* * Add a single value to the t-digest. This is not very efficient, as it has * to deserialize the t-digest into the in-memory aggstate representation * and serialize it back for each call, but it's convenient and acceptable * for some use cases. * * When efficiency is important, it may be possible to use the batch variant * by first aggregating the updates into a t-digest, and then merging that * into an existing t-digest in one step using tdigest_union_double_increment. * * The compact flag controls final compaction only. Adding to a full buffer * still triggers compaction through tdigest_add or tdigest_add_centroid. * * This is similar to hll_add, while the "union" is more like hll_union. */ Datum tdigest_add_double_increment(PG_FUNCTION_ARGS) { tdigest_aggstate_t *state; tdigest_t *result; bool compact; /* the flag determines whether the result gets compacted */ if (PG_ARGISNULL(3)) elog(ERROR, "compact flag must not be NULL"); compact = PG_GETARG_BOOL(3); /* * We want to skip NULL values altogether - we return either the existing * t-digest (if it already exists) or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* if there's no digest allocated, create it now */ if (PG_ARGISNULL(0)) { int compression; /* * We don't require compression, but only when there is an existing * t-digest value. Make sure the value was supplied. */ if (PG_ARGISNULL(2)) elog(ERROR, "compression value not supplied, but t-digest is NULL"); compression = PG_GETARG_INT32(2); check_compression(compression); state = tdigest_aggstate_allocate(0, 0, compression, 0); } else { tdigest_t *digest; digest = tdigest_prepare(PG_GETARG_DATUM(0), false); state = tdigest_digest_to_aggstate(digest); PG_FREE_IF_COPY(digest, 0); } tdigest_add(state, PG_GETARG_FLOAT8(1)); AssertCheckTDigestAggState(state); result = tdigest_aggstate_to_digest(state, compact); tdigest_aggstate_free(state); PG_RETURN_POINTER(result); } /* * Add an array of values to the t-digest. This amortizes the overhead of * deserializing and serializing the t-digest, compared to the per-value * version. * * When efficiency is important, it may be possible to use the batch variant * by first aggregating the updates into a t-digest, and then merging that * into an existing t-digest in one step using tdigest_union_double_increment. * * The compact flag controls final compaction only. Adding to a full buffer * still triggers compaction through tdigest_add or tdigest_add_centroid. * * This is similar to hll_add, while the "union" is more like hll_union. */ Datum tdigest_add_double_array_increment(PG_FUNCTION_ARGS) { tdigest_aggstate_t *state; tdigest_t *result; bool compact; ArrayType *array; const double *values; int nvalues; int i; /* the flag determines whether the result gets compacted */ if (PG_ARGISNULL(3)) elog(ERROR, "compact flag must not be NULL"); compact = PG_GETARG_BOOL(3); /* * We want to skip NULL values altogether - we return either the existing * t-digest (if it already exists) or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* if there's no digest allocated, create it now */ if (PG_ARGISNULL(0)) { int compression; /* * We don't require compression, but only when there is an existing * t-digest value. Make sure the value was supplied. */ if (PG_ARGISNULL(2)) elog(ERROR, "compression value not supplied, but t-digest is NULL"); compression = PG_GETARG_INT32(2); check_compression(compression); state = tdigest_aggstate_allocate(0, 0, compression, 0); } else { tdigest_t *digest; digest = tdigest_prepare(PG_GETARG_DATUM(0), false); state = tdigest_digest_to_aggstate(digest); PG_FREE_IF_COPY(digest, 0); } array = PG_GETARG_ARRAYTYPE_P(1); values = array_to_double(array, "an element", &nvalues); for (i = 0; i < nvalues; i++) tdigest_add(state, values[i]); PG_FREE_IF_COPY(array, 1); AssertCheckTDigestAggState(state); result = tdigest_aggstate_to_digest(state, compact); tdigest_aggstate_free(state); PG_RETURN_POINTER(result); } /* * Merge a t-digest into another t-digest. This is somewhat inefficient, as * it has to deserialize the t-digests into the in-memory aggstate values, * and serialize it back for each call, but it's better than doing it for * each individual value (like tdigest_add_double_increment). * * The compact flag controls final compaction only. Adding to a full buffer * still triggers compaction through tdigest_add or tdigest_add_centroid. * * This is similar to hll_union. */ Datum tdigest_union_double_increment(PG_FUNCTION_ARGS) { int i; tdigest_aggstate_t *state; tdigest_t *digest; tdigest_t *result; bool compact; /* the flag determines whether the result gets compacted */ if (PG_ARGISNULL(2)) elog(ERROR, "compact flag must not be NULL"); compact = PG_GETARG_BOOL(2); if (PG_ARGISNULL(0) && PG_ARGISNULL(1)) PG_RETURN_NULL(); else if (PG_ARGISNULL(0)) PG_RETURN_POINTER(PG_GETARG_POINTER(1)); else if (PG_ARGISNULL(1)) PG_RETURN_POINTER(PG_GETARG_POINTER(0)); /* now we know both arguments are non-null */ /* parse the first digest (we'll merge the other one into this) */ digest = tdigest_prepare(PG_GETARG_DATUM(0), false); state = tdigest_digest_to_aggstate(digest); PG_FREE_IF_COPY(digest, 0); AssertCheckTDigestAggState(state); /* parse the second digest */ digest = tdigest_prepare(PG_GETARG_DATUM(1), false); AssertCheckTDigest(digest); /* copy data from the tdigest into the aggstate */ for (i = 0; i < digest->ncentroids; i++) { CHECK_FOR_INTERRUPTS(); tdigest_add_centroid(state, digest->centroids[i].mean, digest->centroids[i].count); } PG_FREE_IF_COPY(digest, 1); AssertCheckTDigestAggState(state); result = tdigest_aggstate_to_digest(state, compact); tdigest_aggstate_free(state); PG_RETURN_POINTER(result); } /* * Comparator, ordering the centroids by mean value. * * When the mean is the same, we try ordering the centroids by count. * * In principle, centroids with the same mean represent the same value, * but we still need to care about the count to allow rebalancing the * centroids later. */ static int centroid_cmp(const void *a, const void *b) { double ma, mb; centroid_t *ca = (centroid_t *) a; centroid_t *cb = (centroid_t *) b; ma = ca->mean; mb = cb->mean; if (ma < mb) return -1; else if (ma > mb) return 1; if (ca->count < cb->count) return -1; else if (ca->count > cb->count) return 1; return 0; } /* * Parsing of the textual t-digest representation. * * We can't use sscanf, because it does not report overflows in any way - the * value simply saturates to the maximum for the data type. That's a problem * for the count, where the saturated value is a perfectly valid count, so we * can't detect it after the fact. Use strtoll/strtod, which do set errno. * * All of these advance the pointer past the parsed part on success, and never * return on failure. */ /* * Match a literal string, after skipping (optional) leading space. */ static void parse_str(char **ptr, const char *value, bool space) { char *str = *ptr; size_t len = strlen(value); /* if requested, skip the one initial space character */ if (space) { if (isspace((unsigned char) *str)) str++; else ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("failed to parse t-digest value, missing space"))); } /* at this point there must be no whitespace */ if (isspace((unsigned char) *str)) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("failed to parse t-digest value, unexpected space"))); /* the prefix should match our string */ if (strncmp(str, value, len) != 0) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("failed to parse t-digest value, expected \"%s\"", value))); *ptr = str + len; } /* * Parse an int64 value, and make sure it's in range. */ static int64 parse_int64(char **ptr, const char *field) { char *endptr; int64 value; errno = 0; value = strtoi64(*ptr, &endptr, 10); if (endptr == *ptr) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("failed to parse %s of a t-digest", field))); if (errno == ERANGE) ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("%s of a t-digest is out of range for bigint", field))); *ptr = endptr; return value; } /* * Parse an int32 value, and make sure it's in range. * * Parse it as int64 first, so that we can range check it before narrowing it * down, instead of relying on an implementation-defined narrowing conversion. */ static int32 parse_int32(char **ptr, const char *field) { int64 value = parse_int64(ptr, field); if (value < PG_INT32_MIN || value > PG_INT32_MAX) ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("%s of a t-digest is out of range for integer", field))); return (int32) value; } /* * Parse a double value, and make sure it's in range. */ static double parse_double(char **ptr, const char *field) { char *endptr; double value; errno = 0; value = strtod(*ptr, &endptr); if (endptr == *ptr) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("failed to parse %s of a t-digest", field))); /* ERANGE may also signal a nonzero subnormal, which we can store. */ if ((errno == ERANGE) && ((value == 0.0) || !isfinite(value))) ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("%s of a t-digest is out of range for double precision", field))); *ptr = endptr; return value; } Datum tdigest_in(PG_FUNCTION_ARGS) { int i; char *str = PG_GETARG_CSTRING(0); tdigest_t *digest = NULL; size_t slen; /* t-digest header fields */ int32 flags; int64 count, total_count; int compression; int ncentroids; char *ptr; slen = strlen(str); ptr = str; parse_str(&ptr, "flags", false); flags = parse_int32(&ptr, "flags"); parse_str(&ptr, "count", true); count = parse_int64(&ptr, "count"); parse_str(&ptr, "compression", true); compression = parse_int32(&ptr, "compression"); parse_str(&ptr, "centroids", true); ncentroids = parse_int32(&ptr, "number of centroids"); if ((flags & ~TDIGEST_VALID_FLAGS) != 0) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("invalid flags for t-digest"))); if ((compression < MIN_COMPRESSION) || (compression > MAX_COMPRESSION)) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("compression for t-digest must be in [%d, %d]", MIN_COMPRESSION, MAX_COMPRESSION))); if (count <= 0) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("count value for the t-digest must be positive"))); if (ncentroids <= 0) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("number of centroids for the t-digest must be positive"))); if (ncentroids > BUFFER_SIZE(compression)) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("number of centroids for the t-digest exceeds buffer size"))); digest = tdigest_allocate(ncentroids); digest->flags = flags; digest->count = count; digest->ncentroids = ncentroids; digest->compression = compression; total_count = 0; ncentroids = 0; for (i = 0; i < digest->ncentroids; i++) { double mean; CHECK_FOR_INTERRUPTS(); parse_str(&ptr, "(", true); mean = parse_double(&ptr, "mean of a centroid"); parse_str(&ptr, ",", false); count = parse_int64(&ptr, "count of a centroid"); parse_str(&ptr, ")", false); if (!isfinite(mean)) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("mean value for all centroids in a t-digest must be valid"))); if (count <= 0) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("count value for all centroids in a t-digest must be positive"))); else if (count > digest->count) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("count value of a centroid exceeds total count"))); digest->centroids[i].count = count; digest->centroids[i].mean = mean; /* * track the total count so that we can check later * * Make sure the count does not overflow at any point. It could * overflow and then wrap around to the expected total, but it would * still cause an issue. */ if (pg_add_s64_overflow(total_count, count, &total_count)) ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("tdigest count overflow"))); /* * This can't overflow - each centroid has a positive count, and if * the total_count does not overflow, this can't either. */ ncentroids++; /* * The parsing already moved the pointer past the closing parenthesis. * If this is the end of the string, stop parsing, even if we failed to * parse the right number of centroids. */ if (*ptr == '\0') break; /* must not scan past the end of the input string */ Assert(ptr <= str + slen); } /* * Malformed inputs may have the wrong number of centroids, in which case * we either don't consume the whole input (ncentroids too low), or we * don't get all the expected centroids (ncentroids too high). */ if (ptr < str + slen) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("input t-digest value too long"))); if (ncentroids != digest->ncentroids) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("input t-digest value too short"))); /* * If we consumed just the right number of centroids, we must have read * the whole input value exactly. */ Assert(ptr == str + slen); /* check that the total matches */ if (total_count != digest->count) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("total count does not match the data (%lld != %lld)", (long long) total_count, (long long) digest->count))); /* * Make sure we return digest with the new format (it might be the * old format, in which case "mean" fields actually store "sum"). */ tdigest_update_format(digest); AssertCheckTDigest(digest); PG_RETURN_POINTER(digest); } Datum tdigest_out(PG_FUNCTION_ARGS) { int i; tdigest_t *digest = PG_GETARG_TDIGEST(0); StringInfoData str; AssertCheckTDigest(digest); initStringInfo(&str); appendStringInfo(&str, "flags %d count " INT64_FORMAT " compression %d centroids %d", digest->flags, digest->count, digest->compression, digest->ncentroids); /* * If this is an old tdigest with sum values, we'll send those, and * it's up to the reader to fix it. It'll be indicated by not having * the TDIGEST_STORES_MEAN flag. */ for (i = 0; i < digest->ncentroids; i++) { char *tmp = float8out_internal(digest->centroids[i].mean); CHECK_FOR_INTERRUPTS(); appendStringInfo(&str, " (%s, " INT64_FORMAT ")", tmp, digest->centroids[i].count); pfree(tmp); } PG_FREE_IF_COPY(digest, 0); PG_RETURN_CSTRING(str.data); } Datum tdigest_recv(PG_FUNCTION_ARGS) { StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); tdigest_t *digest; int i; int64 count; int64 total_count; int32 flags; int32 compression; int32 ncentroids; flags = pq_getmsgint(buf, sizeof(int32)); /* make sure the t-digest format is supported */ if ((flags != 0) && (flags != TDIGEST_STORES_MEAN)) elog(ERROR, "unsupported t-digest on-disk format"); count = pq_getmsgint64(buf); compression = pq_getmsgint(buf, sizeof(int32)); ncentroids = pq_getmsgint(buf, sizeof(int32)); if ((compression < MIN_COMPRESSION) || (compression > MAX_COMPRESSION)) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("compression for t-digest must be in [%d, %d]", MIN_COMPRESSION, MAX_COMPRESSION))); if (count <= 0) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("count value for the t-digest must be positive"))); if (ncentroids <= 0) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("number of centroids for the t-digest must be positive"))); if (ncentroids > BUFFER_SIZE(compression)) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("number of centroids for the t-digest exceeds buffer size"))); digest = tdigest_allocate(ncentroids); digest->flags = flags; digest->count = count; digest->compression = compression; digest->ncentroids = ncentroids; total_count = 0; for (i = 0; i < digest->ncentroids; i++) { CHECK_FOR_INTERRUPTS(); digest->centroids[i].mean = pq_getmsgfloat8(buf); digest->centroids[i].count = pq_getmsgint64(buf); if (!isfinite(digest->centroids[i].mean)) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("mean value for all centroids in a t-digest must be valid"))); if (digest->centroids[i].count <= 0) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("count value for all centroids in a t-digest must be positive"))); else if (digest->centroids[i].count > digest->count) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("count value of a centroid exceeds total count"))); /* * track the total count so that we can check later * * Make sure the count does not overflow at any point. It could * overflow and then wrap around to the expected total, but it would * still cause an issue. */ if (pg_add_s64_overflow(total_count, digest->centroids[i].count, &total_count)) ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("tdigest count overflow"))); } /* check that the total matches */ if (total_count != digest->count) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("total count does not match the data (%lld != %lld)", (long long) total_count, (long long) digest->count))); /* * Make sure we return digest with the new format (it might be the * old format, in which case "mean" fields actually store "sum"). */ tdigest_update_format(digest); AssertCheckTDigest(digest); PG_RETURN_POINTER(digest); } Datum tdigest_send(PG_FUNCTION_ARGS) { tdigest_t *digest = PG_GETARG_TDIGEST(0); StringInfoData buf; int i; pq_begintypsend(&buf); pq_sendint(&buf, digest->flags, 4); pq_sendint64(&buf, digest->count); pq_sendint(&buf, digest->compression, 4); pq_sendint(&buf, digest->ncentroids, 4); for (i = 0; i < digest->ncentroids; i++) { CHECK_FOR_INTERRUPTS(); pq_sendfloat8(&buf, digest->centroids[i].mean); pq_sendint64(&buf, digest->centroids[i].count); } PG_FREE_IF_COPY(digest, 0); PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); } /* * tdigest_is_valid * Check the t-digest passes the same sanity checks as the input funcs. * * The input functions (tdigest_in and tdigest_recv) validate the values, so * it's not possible to construct an invalid digest through them. But digests * stored by older versions of the extension (which did not have all those * checks) are not re-validated when read back, so they may be broken in * various ways - bogus flags, compression out of range, centroid counts not * adding up to the total count, and so on. This function performs the same * checks on an existing value, allowing users to find such digests. * * Returns true if the digest passes all the checks, false otherwise. Never * raises an error for an invalid digest (that's the whole point). * * The one check the input functions don't need to do explicitly is on the * length of the value - they build the digest from the parsed header, so it * always matches. For an existing value we have to check the value really is * long enough for the centroids the header promises, otherwise we'd read past * the end of it. * * XXX Keep this in sync with the checks in tdigest_in and tdigest_recv. */ static bool tdigest_is_valid_internal(tdigest_t *digest) { int i; int64 total_count; Size vlen; Size expected; vlen = VARSIZE_ANY(digest); /* * We need at least the header, otherwise we can't even look at the fields * describing the rest of the value. */ if (vlen < offsetof(tdigest_t, centroids)) return false; /* make sure the t-digest format is supported */ if ((digest->flags & ~TDIGEST_VALID_FLAGS) != 0) return false; if ((digest->compression < MIN_COMPRESSION) || (digest->compression > MAX_COMPRESSION)) return false; if (digest->count <= 0) return false; if (digest->ncentroids <= 0) return false; if (digest->ncentroids > BUFFER_SIZE(digest->compression)) return false; /* * The header determines how long the value has to be, so make sure it * really is that long before reading any of the centroids. * * The number of centroids is limited by the buffer size (checked above), * so this can't overflow. */ expected = offsetof(tdigest_t, centroids) + digest->ncentroids * sizeof(centroid_t); if (vlen != expected) return false; total_count = 0; for (i = 0; i < digest->ncentroids; i++) { CHECK_FOR_INTERRUPTS(); if (!isfinite(digest->centroids[i].mean)) return false; if (digest->centroids[i].count <= 0) return false; else if (digest->centroids[i].count > digest->count) return false; /* * track the total count so that we can check later * * Make sure the count does not overflow at any point. It could * overflow and then wrap around to the expected total, but it would * still cause an issue. */ if (pg_add_s64_overflow(total_count, digest->centroids[i].count, &total_count)) return false; } /* check that the total matches */ if (total_count != digest->count) return false; return true; } /* * Wrapper for tdigest_is_valid_internal(), so that we can free the digest. */ Datum tdigest_is_valid(PG_FUNCTION_ARGS) { tdigest_t *digest = PG_GETARG_TDIGEST(0); bool valid = tdigest_is_valid_internal(digest); PG_FREE_IF_COPY(digest, 0); PG_RETURN_BOOL(valid); } Datum tdigest_count(PG_FUNCTION_ARGS) { int64 result; tdigest_t *digest = PG_GETARG_TDIGEST(0); result = digest->count; PG_FREE_IF_COPY(digest, 0); PG_RETURN_INT64(result); } /* * tdigest_to_json * Transform the tdigest into a JSON value. * * We make sure to always print mean, even for tdigests in the older format * storing sum for centroids. Otherwise the "means" key would be confusing. * But we don't call tdigest_update_format, and instead we simply update the * flags and convert the sum/mean values. * * The centroids are stored in two separate arrays - one for means, one for * counts. That makes it easier to process, because it's clear the i-th * in each array is for i-th centroid. We might store it in a single array, * but then we'd have to walk it in pairs. And it'd mix float and int * values in the same array. * * The arrays are named in plural ("means" and "counts"), so that the array * of per-centroid counts does not collide with the total count. Duplicate * keys are not strictly forbidden by JSON, but most parsers keep just one * of the values (jsonb keeps the last one), losing the other. */ Datum tdigest_to_json(PG_FUNCTION_ARGS) { int i; StringInfoData str; text *result; tdigest_t *digest = PG_GETARG_TDIGEST(0); int32 flags = digest->flags; initStringInfo(&str); appendStringInfoChar(&str, '{'); flags |= TDIGEST_STORES_MEAN; appendStringInfo(&str, "\"flags\": %d, ", flags); appendStringInfo(&str, "\"count\": " INT64_FORMAT ", ", digest->count); appendStringInfo(&str, "\"compression\": %d, ", digest->compression); appendStringInfo(&str, "\"centroids\": %d, ", digest->ncentroids); appendStringInfoString(&str, "\"means\": ["); for (i = 0; i < digest->ncentroids; i++) { double mean = digest->centroids[i].mean; char *tmp; CHECK_FOR_INTERRUPTS(); if (i > 0) appendStringInfoString(&str, ", "); /* * When the TDIGEST_STORES_MEAN flag is not set, the value is * actually a sum, so convert it to mean now. We have to check the * digest->flags, not the local variable. */ if (! (digest->flags & TDIGEST_STORES_MEAN)) mean = mean / digest->centroids[i].count; tmp = float8out_internal(mean); /* don't print insignificant zeroes to the right of decimal point */ appendStringInfo(&str, "%s", tmp); pfree(tmp); } appendStringInfoString(&str, "], "); appendStringInfoString(&str, "\"counts\": ["); for (i = 0; i < digest->ncentroids; i++) { CHECK_FOR_INTERRUPTS(); if (i > 0) appendStringInfoString(&str, ", "); appendStringInfo(&str, INT64_FORMAT, digest->centroids[i].count); } appendStringInfoString(&str, "]"); appendStringInfoChar(&str, '}'); result = cstring_to_text(str.data); /* free the local buffer, and possibly the detoasted digest copy */ pfree(str.data); PG_FREE_IF_COPY(digest, 0); PG_RETURN_TEXT_P(result); } /* * tdigest_to_array * Transform the tdigest into an array of double values. * * The whole digest is stored in a single "double precision" array, which * may be a bit confusing and perhaps fragile if more fields need to be * added in the future. The initial elements are flags, count (number of * items added to the digest), compression (determines the limit on number * of centroids) and current number of centroids, followed by a stream of values * encoding the centroids in pairs of (mean, count). * * We make sure to always print mean, even for tdigests in the older format * storing sum for centroids. Otherwise the "mean" key would be confusing. * But we don't call tdigest_update_format, and instead we simply update the * flags and convert the sum/mean values. */ Datum tdigest_to_array(PG_FUNCTION_ARGS) { int i, idx; tdigest_t *digest = PG_GETARG_TDIGEST(0); int32 flags = digest->flags; ArrayType *array; double *values; int nvalues; flags |= TDIGEST_STORES_MEAN; /* number of values to store in the array */ nvalues = 4 + (digest->ncentroids * 2); /* write results directly into a correctly sized ArrayType */ array = double_array_allocate(nvalues); values = (double *) ARR_DATA_PTR(array); idx = 0; values[idx++] = flags; values[idx++] = digest->count; values[idx++] = digest->compression; values[idx++] = digest->ncentroids; for (i = 0; i < digest->ncentroids; i++) { double mean = digest->centroids[i].mean; CHECK_FOR_INTERRUPTS(); /* * When the TDIGEST_STORES_MEAN flag is not set, the value is * actually a sum, so convert it to mean now. We have to check the * digest->flags, not the local variable. */ if (! (digest->flags & TDIGEST_STORES_MEAN)) mean = mean / digest->centroids[i].count; /* don't print insignificant zeroes to the right of decimal point */ values[idx++] = mean; values[idx++] = digest->centroids[i].count; } Assert(idx == nvalues); PG_FREE_IF_COPY(digest, 0); PG_RETURN_ARRAYTYPE_P(array); } Datum tdigest_add_double_trimmed(PG_FUNCTION_ARGS) { tdigest_aggstate_t *state; MemoryContext aggcontext; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_add_double_trimmed called in non-aggregate context"); /* * We want to skip NULL values altogether - we return either the existing * t-digest (if it already exists) or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* if there's no digest allocated, create it now */ if (PG_ARGISNULL(0)) { MemoryContext oldcontext; int compression; double low, high; /* the trim thresholds are required to create the aggregate state */ if (PG_ARGISNULL(3) || PG_ARGISNULL(4)) elog(ERROR, "trim thresholds must not be NULL"); low = PG_GETARG_FLOAT8(3); high = PG_GETARG_FLOAT8(4); if (PG_ARGISNULL(2)) elog(ERROR, "compression must not be NULL"); compression = PG_GETARG_INT32(2); check_compression(compression); check_trim_values(low, high); oldcontext = MemoryContextSwitchTo(aggcontext); state = tdigest_aggstate_allocate(0, 0, compression, 0); state->trim_low = low; state->trim_high = high; MemoryContextSwitchTo(oldcontext); } else state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); tdigest_add(state, PG_GETARG_FLOAT8(1)); AssertCheckTDigestAggState(state); PG_RETURN_POINTER(state); } Datum tdigest_add_double_count_trimmed(PG_FUNCTION_ARGS) { int64 i; int64 count; tdigest_aggstate_t *state; MemoryContext aggcontext; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_add_double_count_trimmed called in non-aggregate context"); /* * We want to skip NULL values altogether - we return either the existing * t-digest (if it already exists) or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* if there's no digest allocated, create it now */ if (PG_ARGISNULL(0)) { MemoryContext oldcontext; int compression; double low, high; /* the trim thresholds are required to create the aggregate state */ if (PG_ARGISNULL(4) || PG_ARGISNULL(5)) elog(ERROR, "trim thresholds must not be NULL"); low = PG_GETARG_FLOAT8(4); high = PG_GETARG_FLOAT8(5); if (PG_ARGISNULL(3)) elog(ERROR, "compression must not be NULL"); compression = PG_GETARG_INT32(3); check_compression(compression); check_trim_values(low, high); oldcontext = MemoryContextSwitchTo(aggcontext); state = tdigest_aggstate_allocate(0, 0, compression, 0); state->trim_low = low; state->trim_high = high; MemoryContextSwitchTo(oldcontext); } else state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); if (PG_ARGISNULL(2)) count = 1; else count = PG_GETARG_INT64(2); /* can't add values with non-positive counts */ if (count <= 0) elog(ERROR, "invalid count value %lld, must be a positive value", (long long) count); /* * When adding more values than would fit into an empty buffer (and * thus likely causing too many compactions), we instead add them as * properly sized centroids. * * This is much faster, because the centroids can be generated in one go, * so there are only very few compactions. */ if (count > BUFFER_SIZE(state->compression)) { tdigest_add_generated(state, PG_GETARG_FLOAT8(1), count); count = 0; } /* * If there are only a few values, just add them one by one, so that * we do proper compaction and sizing of centroids. Otherwise we might end * up with oversized centroids on the tails etc. */ for (i = 0; i < count; i++) tdigest_add(state, PG_GETARG_FLOAT8(1)); AssertCheckTDigestAggState(state); PG_RETURN_POINTER(state); } /* * Add a value to the tdigest (create one if needed). Transition function * for tdigest aggregate with a single value. */ Datum tdigest_add_digest_trimmed(PG_FUNCTION_ARGS) { int i; tdigest_aggstate_t *state; tdigest_t *digest; MemoryContext aggcontext; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_add_digest_trimmed called in non-aggregate context"); /* * We want to skip NULL values altogether - we return either the existing * t-digest (if it already exists) or NULL. */ if (PG_ARGISNULL(1)) { if (PG_ARGISNULL(0)) PG_RETURN_NULL(); /* if there already is a state accumulated, don't forget it */ PG_RETURN_DATUM(PG_GETARG_DATUM(0)); } /* make sure we have a detoasted copy */ digest = tdigest_prepare(PG_GETARG_DATUM(1), false); /* if there's no aggregate state allocated, create it now */ if (PG_ARGISNULL(0)) { MemoryContext oldcontext; double low, high; /* the trim thresholds are required to create the aggregate state */ if (PG_ARGISNULL(2) || PG_ARGISNULL(3)) elog(ERROR, "trim thresholds must not be NULL"); low = PG_GETARG_FLOAT8(2); high = PG_GETARG_FLOAT8(3); check_trim_values(low, high); oldcontext = MemoryContextSwitchTo(aggcontext); state = tdigest_aggstate_allocate(0, 0, digest->compression, digest->ncentroids); state->trim_low = low; state->trim_high = high; MemoryContextSwitchTo(oldcontext); } else state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); /* * Keep the compression chosen when the state was created, even if the * input uses a different setting. Compaction can merge incoming * centroids, but cannot split them to recover lost detail. */ for (i = 0; i < digest->ncentroids; i++) { CHECK_FOR_INTERRUPTS(); tdigest_add_centroid(state, digest->centroids[i].mean, digest->centroids[i].count); } AssertCheckTDigestAggState(state); PG_FREE_IF_COPY(digest, 1); PG_RETURN_POINTER(state); } /* * Convert a (count * fraction) product back to a count. * * The product is calculated in double, so for very high counts it may end up * outside the int64 range - even for frac = 1.0. For example * * 1.000000 * 9223372036854775296 = 9223372036854775808.000000 * * which is above both the count and INT64_MAX (9223372036854775807). That * makes the conversion undefined. It can't lead to underflow/overflow, as * we're not using this to access memory, but it might lead to bogus results * (e.g. NULL instead of the correct result). * * Clamp it to the [0, maxvalue] range, which is the only range that makes * sense anyway. * * Note: We know the frac value is in [0.0, 1.0], but we don't rely on that * here. We'll clamp it to [0, maxvalue]. */ static int64 double_to_int64(double value, int64 maxvalue) { /* paranoia: we should not get NaN values here */ if (isnan(value)) return 0; /* clamp it to the [0, count] range */ if (value < 0) return 0; /* * The comparison is done in double on purpose. If we did it as int64, * it might already overflow and wrap. Converting count to double may * round it up to 2^63, but that should be fine - the comparison is * still correct, and we return count for anything that large. */ if (value >= (double) maxvalue) return maxvalue; /* ok, should be safe to cast */ return (int64) value; } /* * How many items of a centroid to use for the aggregate? * * Calculates the number of items of a centroid that fall into the * [count_low, count_high) range of items, with count_done items preceding * the centroid. * * The centroids in the middle of the range are included as a whole, but the * first and last one may be cut in half by the boundary, in which case only * part of the centroid is included. */ static int64 tdigest_trimmed_count(centroid_t *centroid, int64 count_done, int64 count_low, int64 count_high) { int64 count_add; /* Assume the whole centroid falls into the range. */ count_add = centroid->count; /* * If we haven't reached the low threshold yet, skip appropriate * part of the centroid. * * (count_low - count_done) is how far we're from the low threshold, so a * positive value is how many items we still need to "skip" (capped to * size of the centroid). A negative value means we're past the low * threshold, and the centroid may be in the range (unless it's past the * high threshold too). */ count_add -= Min(Max(0, count_low - count_done), count_add); /* * If we have reached the upper threshold, ignore the overflowing * part of the centroid. * * The items we still have start at count_low (or at the beginning of * the centroid, whichever comes later), not at count_done - the part * below count_low was already removed by the preceding step. Don't * count that part a second time, i.e. don't start at the beginning * of the centroid. */ count_add = Min(Max(0, count_high - Max(count_done, count_low)), count_add); return count_add; } /* * Calculate trimmed aggregates from centroids. * * Returns the trimmed mean, the trimmed sum, and the number of items in the * trimmed range. * * The obvious way to calculate the mean is to add (mean * count) for all the * centroids in the range, and then divide the sum by the number of items. * That however may overflow to infinity, even when the mean is perfectly * representable - the average of 1e307 values is 1e307, but the sum of many * such values is not. And once one partial sum saturates to +Infinity and * another one to -Infinity, the accumulator turns into NaN. * * The mean of a digest is always within the range of the centroid means, so * it should never overflow. We just have to calculate it in a way that does * not overflow either, i.e. as a weighted average of the centroid means, * similarly to how tdigest_compact() does when merging centroids * * mean += centroids[i].mean * (count_add / total_count) * * The weights are non-negative and add up to 1.0 (albeit maybe not perfectly, * due to limited precision of float8), so the running mean stays within the * range of the centroid means. * * This needs the total number of items in the range up front, so we walk the * centroids twice. The first pass determines the first and last centroid of * the range, and adds up the (int64) counts - that can't overflow, as the * counts add up to the total count of the digest. * * The sum, on the other hand, may legitimately exceed the float8 range, so we * keep accumulating it the simple way (which is also exact), and leave it to * the caller to complain about the overflow. * * FIXME I believe this has two issues: * * First, the count_low/count_high round in opposite direction, but may round * to the same value when the product hits an integer exactly (when low==high). * If this happens, we end up returning NULL a bit later (count_total==0). * Maybe we should require (low < high)? Would that solve the issue, or are * there other ways to still hit this? Or maybe we should do something like * count_high = Max(count_high, count_low + 1), to force non-empty ranges? * * Second, the code assumes the whole centroid is located at the mean, and does * not use the interpolation (assuming half the items are below/above the mean). * Maybe it should do something like the quantile/quantile-of code? */ static void tdigest_trimmed_agg(centroid_t *centroids, int ncentroids, int64 count, double low, double high, double *meanp, double *sump, int64 *countp) { int i; int first = 0, last = -1; double mean = 0, sum = 0; int64 count_done = 0, count_first = 0, count_total = 0, count_low, count_high; /* translate the percentiles to counts */ count_low = double_to_int64(floor(count * low), count); count_high = double_to_int64(ceil(count * high), count); /* verify sane range */ Assert((count_low <= count_high) && (0 <= count_low) && (count_high <= count)); *meanp = 0; *sump = 0; *countp = 0; /* * Find the first and last centroid of the range, and the number of items * the range contains. */ for (i = 0; i < ncentroids; i++) { int64 count_add = tdigest_trimmed_count(¢roids[i], count_done, count_low, count_high); CHECK_FOR_INTERRUPTS(); if (count_add > 0) { /* remember where the range starts, including the item offset */ if (last == -1) { first = i; count_first = count_done; } last = i; count_total += count_add; } /* consider the whole centroid processed */ count_done += centroids[i].count; /* break once we cross the high threshold */ if (count_done >= count_high) break; } /* no items in the range, the callers return NULL in that case */ if (count_total == 0) return; /* now walk just the centroids in the range, and calculate mean / sum */ count_done = count_first; for (i = first; i <= last; i++) { int64 count_add = tdigest_trimmed_count(¢roids[i], count_done, count_low, count_high); CHECK_FOR_INTERRUPTS(); /* consider the whole centroid processed */ count_done += centroids[i].count; /* increment the mean / sum */ mean += centroids[i].mean * (count_add / (double) count_total); sum += centroids[i].mean * count_add; Assert(!isnan(mean)); } // Assert((centroids[first].mean <= mean) && (mean <= centroids[last].mean)); /* * paranoia: handle possible overflow of the mean by clamping it to the * valid range using the mean of the first/last centroid * * Maybe overflow is not the right term, but what can happen easily is * the value "drifting" outside the valid range with very high counts, * even if the centroids have the exact same mean. For example with * centroids like (1, 22522773787004704) and (1, 22627252715671912) the * mean will drift to 1.0000000000000002. * * To confirm, uncomment the assert above, and rerun tests. There's a * query that triggers it. */ mean = Max(Min(centroids[last].mean, mean), centroids[first].mean); *meanp = mean; *sump = sum; *countp = count_total; } /* * Return the trimmed sum, and complain if it overflowed - just like the * regular float8 arithmetic does, instead of silently returning infinity. * * Infinity is a perfectly valid result if some of the input values were * infinite, in which case the mean is infinite too. * * XXX This is what float8_mul() does, but that's only available on PG12+. */ static double tdigest_trimmed_sum_value(double sum, double mean, int64 count) { /* * The accumulator may overflow even when the sum itself is perfectly * representable (e.g. with values of the opposite sign). The mean can't * overflow, so recalculate the sum from that. */ if (!isfinite(sum) && isfinite(mean)) sum = mean * (double) count; if (isinf(sum) && !isinf(mean)) ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("value out of range: overflow"))); return sum; } /* * Compute percentile from a tdigest. Final function for tdigest aggregate * with a single percentile. */ Datum tdigest_trimmed_avg(PG_FUNCTION_ARGS) { tdigest_aggstate_t *state; MemoryContext aggcontext; double mean; double sum; int64 count; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_trimmed_avg called in non-aggregate context"); /* if there's no digest, return NULL */ if (PG_ARGISNULL(0)) PG_RETURN_NULL(); state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); /* make sure the centroids are sorted */ tdigest_sort(state); tdigest_trimmed_agg(state->centroids, state->ncentroids, state->count, state->trim_low, state->trim_high, &mean, &sum, &count); if (count > 0) PG_RETURN_FLOAT8(mean); PG_RETURN_NULL(); } /* * Compute percentile from a tdigest. Final function for tdigest aggregate * with a single percentile. */ Datum tdigest_trimmed_sum(PG_FUNCTION_ARGS) { tdigest_aggstate_t *state; MemoryContext aggcontext; double mean; double sum; int64 count; /* cannot be called directly because of internal-type argument */ if (!AggCheckCallContext(fcinfo, &aggcontext)) elog(ERROR, "tdigest_trimmed_sum called in non-aggregate context"); /* if there's no digest, return NULL */ if (PG_ARGISNULL(0)) PG_RETURN_NULL(); state = (tdigest_aggstate_t *) PG_GETARG_POINTER(0); /* make sure the centroids are sorted */ tdigest_sort(state); tdigest_trimmed_agg(state->centroids, state->ncentroids, state->count, state->trim_low, state->trim_high, &mean, &sum, &count); if (count > 0) PG_RETURN_FLOAT8(tdigest_trimmed_sum_value(sum, mean, count)); PG_RETURN_NULL(); } /* * Trimmed sum of a single digest (non-aggregate function). */ Datum tdigest_digest_sum(PG_FUNCTION_ARGS) { tdigest_t *digest; double low = PG_GETARG_FLOAT8(1); double high = PG_GETARG_FLOAT8(2); double mean; double sum; int64 count; check_trim_values(low, high); /* tdigest_trimmed_agg expects the centroids sorted by mean */ digest = tdigest_prepare(PG_GETARG_DATUM(0), true); AssertCheckTDigest(digest); tdigest_trimmed_agg(digest->centroids, digest->ncentroids, digest->count, low, high, &mean, &sum, &count); PG_FREE_IF_COPY(digest, 0); if (count > 0) PG_RETURN_FLOAT8(tdigest_trimmed_sum_value(sum, mean, count)); PG_RETURN_NULL(); } /* * Trimmed average of a single digest (non-aggregate function) */ Datum tdigest_digest_avg(PG_FUNCTION_ARGS) { tdigest_t *digest; double low = PG_GETARG_FLOAT8(1); double high = PG_GETARG_FLOAT8(2); double mean; double sum; int64 count; check_trim_values(low, high); /* tdigest_trimmed_agg expects the centroids sorted by mean */ digest = tdigest_prepare(PG_GETARG_DATUM(0), true); AssertCheckTDigest(digest); tdigest_trimmed_agg(digest->centroids, digest->ncentroids, digest->count, low, high, &mean, &sum, &count); PG_FREE_IF_COPY(digest, 0); if (count > 0) PG_RETURN_FLOAT8(mean); PG_RETURN_NULL(); } /* * Return a read-only view of an input FLOAT8 SQL array as C doubles. * * This expects a single-dimensional float8 array, fails otherwise. * * "what" names a single element of the array, with the article, and is used * when reporting a NULL element. The callers pass arrays of different things * (percentiles, hypothetical values, values to add to a digest), and a message * naming the wrong one points at the wrong argument. * * The caller must keep the detoasted array alive while using the view. */ static const double * array_to_double(ArrayType *v, const char *what, int *len) { int nitems, *dims, ndims; Oid element_type; ndims = ARR_NDIM(v); dims = ARR_DIMS(v); nitems = ArrayGetNItems(ndims, dims); /* * Reject empty arrays explicitly. An empty array has ndims = 0, so * without this it would be caught by the single-dimension check below * and reported as a dimensionality problem, which is misleading - the * array is well-formed, it just has nothing in it. */ if (nitems == 0) elog(ERROR, "the array must not be empty"); /* this is a special-purpose function for single-dimensional arrays */ if (ndims != 1) elog(ERROR, "expected a single-dimensional array (dims = %d)", ndims); element_type = ARR_ELEMTYPE(v); /* XXX not sure if really needed (can it actually happen?) */ if (element_type != FLOAT8OID) elog(ERROR, "array_to_double expects FLOAT8 array"); if (array_contains_nulls(v)) elog(ERROR, "NULL not allowed as %s", what); (*len) = nitems; /* Non-NULL float8 elements have the same layout as a C double array. */ return (const double *) ARR_DATA_PTR(v); } /* * Allocate a one-dimensional, non-NULL float8 array for direct result writes. * Array storage uses native doubles even on platforms with pass-by-reference * float8 Datums, so no per-element Datum allocations are needed. */ static ArrayType * double_array_allocate(int nitems) { ArrayType *array; Size size; /* should not happen */ if (nitems <= 0) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("invalid array size (%d)", nitems))); if (nitems > MaxArraySize) ereport(ERROR, (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED), errmsg("array size exceeds the maximum allowed (%d)", (int) MaxArraySize))); size = ARR_OVERHEAD_NONULLS(1) + nitems * sizeof(double); array = (ArrayType *) palloc(size); SET_VARSIZE(array, size); ARR_NDIM(array) = 1; array->dataoffset = 0; ARR_ELEMTYPE(array) = FLOAT8OID; ARR_DIMS(array)[0] = nitems; ARR_LBOUND(array)[0] = 1; return array; }