-- tdigest_percentile() and tdigest_percentile_of() are inverse functions -
-- one interpolates a value from a percentile, the other a percentile from a
-- value - so composing them has to be an identity:
--
--     tdigest_percentile_of(d, tdigest_percentile(d, p)) = p
--
-- This only holds for percentiles that are actually interpolated. Both
-- functions place the mean of a centroid in the middle of it, so below
-- (first centroid count / 2) items and above (last centroid count / 2) items
-- there is nothing to interpolate between and the result is clamped to the
-- first/last mean. Each digest below therefore comes with the range of
-- percentiles that are inside the interpolated part.
--
-- The identity does not hold at the moment, because the two functions do not
-- agree on where the mean of a centroid is. tdigest_compute_quantiles_of()
-- computes the offset of the mean with an int64 division
--
--     count -= (prev->count / 2);
--
-- (note "count" is a double, and the very next statement uses / 2.0), so for
-- centroids with an odd count it is off by half an item, i.e. the round trip
-- is off by 0.5/count. On the ten item digest below that is a whole 0.05:
--
--     p       tdigest_percentile   round trip     expected
--     0.1     1.5                  0.15           0.1
--     0.5     5.5                  0.55           0.5
--     0.9     9.5                  0.95           0.9
\set VERBOSITY terse
SET extra_float_digits = 0;
CREATE TABLE tdigest_roundtrip_digests (id int, descr text, d tdigest,
                                        p_low double precision, p_high double precision);
INSERT INTO tdigest_roundtrip_digests VALUES
    (1, 'ten centroids, one item each',
        'flags 1 count 10 compression 10000 centroids 10 (1, 1) (2, 1) (3, 1) (4, 1) (5, 1) (6, 1) (7, 1) (8, 1) (9, 1) (10, 1)',
        0.05, 0.95),
    (2, 'two centroids, odd + odd',
        'flags 1 count 10 compression 10000 centroids 2 (0, 7) (100, 3)',
        0.35, 0.85),
    (3, 'three centroids, mixed counts',
        'flags 1 count 12 compression 10000 centroids 3 (0, 3) (50, 4) (100, 5)',
        0.125, 0.79),
    (4, 'a large centroid in the middle',
        'flags 1 count 22 compression 10000 centroids 3 (0, 1) (10, 20) (20, 1)',
        0.05, 0.95),
    (5, 'negative and positive means',
        'flags 1 count 8 compression 10000 centroids 4 (-100, 2) (-1, 2) (1, 2) (100, 2)',
        0.125, 0.875);
-- And digests built by the aggregate from actual data. The first one is not
-- compacted at all (every centroid holds a single item), the second one is.
INSERT INTO tdigest_roundtrip_digests
SELECT 6, 'built from 1000 distinct values, no compaction', tdigest(i::double precision, 10000), 0.005, 0.995
  FROM generate_series(1, 1000) s(i);
INSERT INTO tdigest_roundtrip_digests
SELECT 7, 'built from 10000 rows, compression 100', tdigest((i % 1000)::double precision, 100), 0.05, 0.95
  FROM generate_series(1, 10000) s(i);
CREATE TABLE tdigest_roundtrip_percentiles (p double precision);
INSERT INTO tdigest_roundtrip_percentiles
SELECT i / 100.0 FROM generate_series(1, 99) s(i);
-- Report one row per digest that fails the round trip, so that a failure
-- stays readable.
WITH v AS (
    SELECT g.id, g.descr, p.p, tdigest_percentile(g.d, p.p) AS value
      FROM tdigest_roundtrip_digests g, tdigest_roundtrip_percentiles p
     WHERE p.p >= g.p_low AND p.p <= g.p_high
     GROUP BY g.id, g.descr, p.p
), rt AS (
    SELECT v.id, v.descr, v.p, v.value, tdigest_percentile_of(g.d, v.value) AS back
      FROM v, tdigest_roundtrip_digests g
     WHERE g.id = v.id
     GROUP BY v.id, v.descr, v.p, v.value
)
SELECT id, descr, count(*) AS failures, max(abs(back - p)) AS max_error
  FROM rt WHERE abs(back - p) > 1e-12
 GROUP BY id, descr ORDER BY id;
 id | descr | failures | max_error 
----+-------+----------+-----------
(0 rows)

-- The three concrete percentiles from the comment at the top.
WITH v AS (
    SELECT p, tdigest_percentile('flags 1 count 10 compression 10000 centroids 10 (1, 1) (2, 1) (3, 1) (4, 1) (5, 1) (6, 1) (7, 1) (8, 1) (9, 1) (10, 1)'::tdigest, p) AS value
      FROM unnest(ARRAY[0.1, 0.5, 0.9]::double precision[]) AS p
     GROUP BY p
)
SELECT p, value,
       tdigest_percentile_of('flags 1 count 10 compression 10000 centroids 10 (1, 1) (2, 1) (3, 1) (4, 1) (5, 1) (6, 1) (7, 1) (8, 1) (9, 1) (10, 1)'::tdigest, value) AS round_trip
  FROM v GROUP BY p, value ORDER BY p;
  p  | value | round_trip 
-----+-------+------------
 0.1 |   1.5 |        0.1
 0.5 |   5.5 |        0.5
 0.9 |   9.5 |        0.9
(3 rows)

-- The same identity has to hold for the aggregates building the digest from
-- data. Ten items, so the interpolated range is [0.05, 0.95].
WITH v AS (
    SELECT p, tdigest_percentile(i::double precision, 10000, p) AS value
      FROM generate_series(1, 10) s(i), unnest(ARRAY[0.1, 0.5, 0.9]::double precision[]) AS p
     GROUP BY p
)
SELECT v.p, v.value, tdigest_percentile_of(i::double precision, 10000, v.value) AS round_trip
  FROM v, generate_series(1, 10) s(i)
 GROUP BY v.p, v.value ORDER BY v.p;
  p  | value | round_trip 
-----+-------+------------
 0.1 |   1.5 |        0.1
 0.5 |   5.5 |        0.5
 0.9 |   9.5 |        0.9
(3 rows)

WITH v AS (
    SELECT p, tdigest_percentile(i::double precision, 2::bigint, 10000, p) AS value
      FROM generate_series(1, 5) s(i), unnest(ARRAY[0.2, 0.5, 0.8]::double precision[]) AS p
     GROUP BY p
)
SELECT v.p, v.value, tdigest_percentile_of(i::double precision, 2::bigint, 10000, v.value) AS round_trip
  FROM v, generate_series(1, 5) s(i)
 GROUP BY v.p, v.value ORDER BY v.p;
  p  | value | round_trip 
-----+-------+------------
 0.2 |   1.5 |        0.2
 0.5 |     3 |        0.5
 0.8 |   4.5 |        0.8
(3 rows)

DROP TABLE tdigest_roundtrip_digests;
DROP TABLE tdigest_roundtrip_percentiles;