Line data Source code
1 : /*
2 : Copyright (c) 2015 Loïc Séguin-Charbonneau
3 :
4 : Permission is hereby granted, free of charge, to any person obtaining a copy
5 : of this software and associated documentation files (the "Software"), to deal
6 : in the Software without restriction, including without limitation the rights
7 : to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8 : copies of the Software, and to permit persons to whom the Software is
9 : furnished to do so, subject to the following conditions:
10 :
11 : The above copyright notice and this permission notice shall be included in
12 : all copies or substantial portions of the Software.
13 :
14 : THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15 : IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16 : FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17 : AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18 : LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19 : OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
20 : THE SOFTWARE.
21 : */
22 :
23 : /*
24 : * This is an implementation of the t-digest algorithm by Ted Dunning and Otmar
25 : * Ertl. The algorithm is detailed in https://github.com/tdunning/t-digest/blob/master/docs/t-digest-paper/histo.pdf
26 : * and a reference implementation in Java is available at https://github.com/tdunning/t-digest.
27 : *
28 : * The current implementation has also been inspired by the work of Cam
29 : * Davidson-Pilon: https://github.com/CamDavidsonPilon/tdigest.
30 : */
31 :
32 : #include <float.h>
33 : #include <math.h>
34 : #include <stdbool.h>
35 : #include <stdlib.h>
36 : #include "tdigest.h"
37 : #include "tree.h"
38 : #include <assert.h>
39 : #include <stdio.h>
40 : #define arc4random_uniform(x) (rand() % x)
41 :
42 : struct Centroid {
43 : RB_ENTRY(Centroid) entry;
44 : double mean;
45 : size_t count;
46 : };
47 :
48 3 : int centroidcmp(Centroid *c1, Centroid *c2)
49 : {
50 3 : return (c1->mean < c2->mean ? -1 : 1);
51 : }
52 :
53 : struct TDigest {
54 : RB_HEAD(CentroidTree, Centroid) C;
55 : size_t count;
56 : size_t ncentroids;
57 : double delta;
58 : unsigned int K;
59 : size_t ncompressions;
60 : };
61 :
62 89 : RB_GENERATE(CentroidTree, Centroid, entry, centroidcmp)
63 :
64 3 : TDigest* TDigest_create(double delta, unsigned int K)
65 : {
66 3 : TDigest *digest = calloc(1, sizeof(TDigest));
67 3 : RB_INIT(&(digest->C));
68 3 : digest->delta = delta;
69 3 : digest->K = K;
70 3 : digest->count = 0;
71 3 : digest->ncentroids = 0;
72 3 : digest->ncompressions = 0;
73 3 : return digest;
74 : }
75 :
76 3 : void TDigest_destroy(TDigest* digest)
77 : {
78 : Centroid *c, *nxt;
79 :
80 8 : for (c = RB_MIN(CentroidTree, &(digest->C)); c != NULL; c = nxt) {
81 5 : nxt = RB_NEXT(CentroidTree, &(digest->C), c);
82 5 : RB_REMOVE(CentroidTree, &(digest->C), c);
83 5 : free(c);
84 : }
85 3 : free(digest);
86 3 : }
87 :
88 5 : TDigest * TDigest_add(TDigest *pdigest, double x, size_t w)
89 : {
90 5 : TDigest *rd = NULL;
91 : Centroid *cj;
92 5 : cj = TDigest_find_closest_centroid(pdigest, x, w);
93 :
94 5 : if (cj != NULL) {
95 : // Add the data point to the selected centroid.
96 0 : RB_REMOVE(CentroidTree, &((pdigest)->C), cj);
97 0 : Centroid_add(cj, x, w);
98 0 : RB_INSERT(CentroidTree, &((pdigest)->C), cj);
99 : } else {
100 5 : Centroid *c = Centroid_create(x, w);
101 5 : RB_INSERT(CentroidTree, &((pdigest)->C), c);
102 5 : pdigest->ncentroids += 1;
103 : }
104 :
105 5 : (pdigest)->count += w;
106 :
107 5 : if ((pdigest)->ncentroids > (pdigest)->K / (pdigest)->delta) {
108 0 : rd = TDigest_compress(pdigest);
109 : }
110 5 : return rd;
111 : }
112 :
113 5 : Centroid *TDigest_find_closest_centroid(TDigest *digest, double x, size_t w)
114 : {
115 : // Find all the centroids whose mean is the closest to x. Return the number
116 : // of centroids that are closest to x.
117 5 : if (digest->ncentroids == 0) {
118 3 : return NULL;
119 : }
120 :
121 2 : double z, min_distance = DBL_MAX;
122 2 : double sum = 0.0;
123 2 : Centroid *c, *lower_closest, *upper_closest = NULL;
124 2 : lower_closest = RB_MIN(CentroidTree, &(digest->C));
125 :
126 : // Start at the beginning of the tree keep going as long as the distance to
127 : // x decreases.
128 4 : for (c = lower_closest; c != NULL; c = RB_NEXT(CentroidTree, &(digest->C), c)) {
129 : // Sum the counts of centroids with mean smaller than x. This is used
130 : // to compute the quantile.
131 3 : if (c->mean < x) {
132 2 : sum += c->count;
133 : }
134 3 : z = fabs(c->mean - x);
135 3 : if (z < min_distance) {
136 2 : min_distance = z;
137 2 : lower_closest = c;
138 1 : } else if (z > min_distance) {
139 1 : upper_closest = c;
140 1 : break;
141 : }
142 : }
143 :
144 : // Start at the lower_closest and choose one of the closer centroids at
145 : // random.
146 : double qc, threshold;
147 2 : double n = 0.0;
148 2 : Centroid *closest = NULL;
149 :
150 4 : for (c = lower_closest; c != upper_closest; c = RB_NEXT(CentroidTree, &(digest->C), c)) {
151 2 : qc = (c->count / 2.0 + sum) / digest->count;
152 2 : sum += c->count;
153 2 : threshold = 4 * digest->count * digest->delta * qc * (1 - qc);
154 2 : if (c->count + w <= threshold) {
155 0 : n++;
156 0 : if (rand() / (double)RAND_MAX < 1.0 / n) {
157 0 : closest = c;
158 : }
159 : }
160 : }
161 2 : return closest;
162 : }
163 :
164 0 : TDigest * TDigest_compress(TDigest *digestp)
165 : {
166 0 : TDigest *new_digest = TDigest_create(digestp->delta, digestp->K);
167 : Centroid *c;
168 : int i, j;
169 0 : while (!RB_EMPTY(&(digestp->C))) {
170 0 : j = arc4random_uniform(digestp->ncentroids);
171 0 : c = RB_MIN(CentroidTree, &(digestp->C));
172 0 : for (i = 0; i < j; i++) {
173 0 : c = RB_NEXT(CentroidTree, &(digestp->C), c);
174 : }
175 0 : RB_REMOVE(CentroidTree, &(digestp->C), c);
176 0 : digestp->count -= c->count;
177 0 : digestp->ncentroids -= 1;
178 0 : TDigest * next_digest = TDigest_add(new_digest, c->mean, c->count);
179 0 : if (next_digest != NULL) {
180 0 : fprintf(stdout, "TDigest_compress recursive\n");
181 0 : TDigest_destroy(new_digest);
182 0 : new_digest = next_digest;
183 0 : assert(0);
184 : }
185 0 : free(c);
186 : }
187 :
188 0 : new_digest->ncompressions = digestp->ncompressions;
189 0 : new_digest->ncompressions++;
190 0 : return new_digest;
191 : }
192 :
193 0 : size_t TDigest_get_ncompressions(TDigest *digest)
194 : {
195 0 : return digest->ncompressions;
196 : }
197 :
198 7 : double TDigest_percentile(TDigest *digest, double q)
199 : {
200 7 : double delta, t = 0;
201 7 : bool first = true;
202 : Centroid *c;
203 7 : q *= digest->count;
204 14 : RB_FOREACH(c, CentroidTree, &(digest->C)) {
205 14 : if (q < t + c->count) {
206 7 : if (first) {
207 3 : return c->mean;
208 4 : } else if (c == RB_MAX(CentroidTree, &(digest->C))) {
209 3 : return c->mean;
210 : } else {
211 1 : double dprev = c->mean - RB_PREV(CentroidTree, &(digest->C),c )->mean;
212 1 : double dnext = RB_NEXT(CentroidTree, &(digest->C), c)->mean - c->mean;
213 1 : delta = (dprev < dnext ? 2*dprev : 2*dnext);
214 : }
215 1 : return c->mean + ((q - t) / c->count - 0.5) * delta;
216 : }
217 7 : t += c->count;
218 7 : first = false;
219 : }
220 0 : return RB_MAX(CentroidTree, &(digest->C))->mean;
221 : }
222 :
223 2 : size_t TDigest_get_ncentroids(TDigest *digest)
224 : {
225 2 : return digest->ncentroids;
226 : }
227 :
228 2 : Centroid *TDigest_get_centroid(TDigest *digest, size_t i)
229 : {
230 2 : Centroid *c = RB_MIN(CentroidTree, &(digest->C));
231 : size_t j;
232 :
233 2 : for (j = 0; j < i; j++) {
234 0 : c = RB_NEXT(CentroidTree, &(digest->C), c);
235 : }
236 :
237 2 : return c;
238 : }
239 :
240 0 : size_t TDigest_get_count(TDigest *digest)
241 : {
242 0 : return digest->count;
243 : }
244 :
245 11 : Centroid* Centroid_create(double x, size_t w)
246 : {
247 11 : Centroid *centroid = malloc(sizeof(Centroid));
248 11 : centroid->count = w;
249 11 : centroid->mean = x;
250 11 : return centroid;
251 : }
252 :
253 10 : void Centroid_add(Centroid *c, double x, size_t w)
254 : {
255 10 : c->count += w;
256 10 : c->mean += w * (x - c->mean) / c->count;
257 10 : }
258 :
259 0 : double Centroid_quantile(Centroid *c, TDigest *digest)
260 : {
261 : Centroid *cj;
262 0 : double quantile = c->count / 2.0;
263 0 : for (cj = RB_PREV(CentroidTree, &(digest->C), c); cj != NULL; cj = RB_PREV(CentroidTree, &(digest->C), cj)) {
264 0 : quantile += cj->count;
265 : }
266 0 : return quantile / digest->count;
267 : }
268 :
269 8 : double Centroid_get_mean(Centroid *c)
270 : {
271 8 : return c->mean;
272 : }
273 :
274 6 : size_t Centroid_get_count(Centroid *c)
275 : {
276 6 : return c->count;
277 : }
278 :
|
