Line data Source code
1 : /*
2 : * SPDX-License-Identifier: MPL-2.0
3 : *
4 : * This Source Code Form is subject to the terms of the Mozilla Public
5 : * License, v. 2.0. If a copy of the MPL was not distributed with this
6 : * file, You can obtain one at http://mozilla.org/MPL/2.0/.
7 : *
8 : * Copyright 2024 MonetDB Foundation;
9 : * Copyright August 2008 - 2023 MonetDB B.V.;
10 : * Copyright 1997 - July 2008 CWI.
11 : */
12 :
13 : #include "monetdb_config.h"
14 : #include "gdk.h"
15 : #include "gdk_private.h"
16 : #include "gdk_calc_private.h"
17 :
18 : /* grouped aggregates
19 : *
20 : * The following functions take two to four input BATs and produce a
21 : * single output BAT.
22 : *
23 : * The input BATs are
24 : * - b, a dense-headed BAT with the values to work on in the tail;
25 : * - g, a dense-headed BAT, aligned with b, with group ids (OID) in
26 : * the tail;
27 : * - e, optional but recommended, a dense-headed BAT with the list of
28 : * group ids in the head(!) (the tail is completely ignored);
29 : * - s, optional, a dense-headed bat with a list of candidate ids in
30 : * the tail.
31 : *
32 : * The tail values of s refer to the head of b and g. Only entries at
33 : * the specified ids are taken into account for the grouped
34 : * aggregates. All other values are ignored. s is compatible with
35 : * the result of BATselect().
36 : *
37 : * If e is not specified, we need to do an extra scan over g to find
38 : * out the range of the group ids that are used. e is defined in such
39 : * a way that it can be either the extents or the histo result from
40 : * BATgroups().
41 : *
42 : * All functions calculate grouped aggregates. There are as many
43 : * groups as there are entries in e. If e is not specified, the
44 : * number of groups is equal to the difference between the maximum and
45 : * minimum values in g.
46 : *
47 : * If a group is empty, the result for that group is nil.
48 : *
49 : * If there is overflow during the calculation of an aggregate, the
50 : * whole operation fails.
51 : *
52 : * If skip_nils is non-zero, a nil value in b is ignored, otherwise a
53 : * nil in b results in a nil result for the group.
54 : */
55 :
56 : /* helper function
57 : *
58 : * This function finds the minimum and maximum group id (and the
59 : * number of groups) and initializes the variables for candidates
60 : * selection.
61 : *
62 : * In case of error, returns an error message.
63 : */
64 : const char *
65 27371 : BATgroupaggrinit(BAT *b, BAT *g, BAT *e, BAT *s,
66 : /* outputs: */
67 : oid *minp, oid *maxp, BUN *ngrpp,
68 : struct canditer *ci)
69 : {
70 27371 : oid min, max;
71 27371 : BUN i, ngrp;
72 27371 : const oid *restrict gids;
73 :
74 27371 : if (b == NULL)
75 : return "b must exist";
76 27371 : canditer_init(ci, b, s);
77 27369 : if (g) {
78 18686 : if (ci->ncand != BATcount(g) ||
79 11858 : (ci->ncand != 0 && ci->seq != g->hseqbase))
80 : return "b with s and g must be aligned";
81 18687 : assert(BATttype(g) == TYPE_oid);
82 : }
83 : if (g == NULL) {
84 : min = 0;
85 : max = 0;
86 : ngrp = 1;
87 18687 : } else if (e == NULL) {
88 3 : if (g->tmaxpos != BUN_NONE) {
89 0 : min = 0;
90 0 : max = BUNtoid(g, g->tmaxpos);
91 : } else {
92 3 : min = oid_nil; /* note that oid_nil > 0! (unsigned) */
93 3 : max = 0;
94 3 : if (BATtdense(g)) {
95 2 : min = g->tseqbase;
96 2 : max = g->tseqbase + BATcount(g) - 1;
97 1 : } else if (g->tsorted) {
98 1 : gids = (const oid *) Tloc(g, 0);
99 : /* find first non-nil */
100 1 : for (i = 0, ngrp = BATcount(g); i < ngrp; i++) {
101 1 : if (!is_oid_nil(gids[i])) {
102 : min = gids[i];
103 : break;
104 : }
105 : }
106 1 : if (!is_oid_nil(min)) {
107 : /* found a non-nil, max must be last
108 : * value (and there is one!) */
109 1 : max = gids[BATcount(g) - 1];
110 : }
111 : } else {
112 0 : const ValRecord *prop;
113 0 : prop = BATgetprop(g, GDK_MAX_BOUND);
114 0 : if (prop != NULL) {
115 0 : assert(prop->vtype == TYPE_oid);
116 0 : min = 0; /* just assume it starts at 0 */
117 0 : max = prop->val.oval - 1; /* bound is exclusive */
118 : } else {
119 : /* we'll do a complete scan */
120 0 : gids = (const oid *) Tloc(g, 0);
121 0 : for (i = 0, ngrp = BATcount(g); i < ngrp; i++) {
122 0 : if (!is_oid_nil(gids[i])) {
123 0 : if (gids[i] < min)
124 : min = gids[i];
125 0 : if (gids[i] > max)
126 : max = gids[i];
127 : }
128 : }
129 : /* note: max < min is possible
130 : * if all groups are nil (or
131 : * BATcount(g)==0) */
132 : }
133 : }
134 : }
135 3 : ngrp = max < min ? 0 : max - min + 1;
136 : } else {
137 18684 : ngrp = BATcount(e);
138 18684 : min = e->hseqbase;
139 18684 : max = e->hseqbase + ngrp - 1;
140 : }
141 27370 : *minp = min;
142 27370 : *maxp = max;
143 27370 : *ngrpp = ngrp;
144 :
145 27370 : return NULL;
146 : }
147 :
148 : /* ---------------------------------------------------------------------- */
149 : /* sum */
150 :
151 : static inline bool
152 293 : samesign(double x, double y)
153 : {
154 293 : return (x >= 0) == (y >= 0);
155 : }
156 :
157 : /* Add two values, producing the sum and the remainder due to limited
158 : * range of floating point arithmetic. This function depends on the
159 : * fact that the sum returns INFINITY in *hi of the correct sign
160 : * (i.e. isinf() returns TRUE) in case of overflow. */
161 : static inline void
162 20477 : twosum(volatile double *hi, volatile double *lo, double x, double y)
163 : {
164 20477 : volatile double yr;
165 :
166 20477 : assert(fabs(x) >= fabs(y));
167 :
168 20477 : *hi = x + y;
169 20477 : yr = *hi - x;
170 20477 : *lo = y - yr;
171 20477 : }
172 :
173 : static inline void
174 9446 : exchange(double *x, double *y)
175 : {
176 9446 : double t = *x;
177 9446 : *x = *y;
178 9446 : *y = t;
179 9446 : }
180 :
181 : /* this function was adapted from https://bugs.python.org/file10357/msum4.py */
182 : BUN
183 1108 : dofsum(const void *restrict values, oid seqb,
184 : struct canditer *restrict ci,
185 : void *restrict results, BUN ngrp, int tp1, int tp2,
186 : const oid *restrict gids,
187 : oid min, oid max, bool skip_nils, bool nil_if_empty)
188 : {
189 1108 : struct pergroup {
190 : int npartials;
191 : int maxpartials;
192 : bool valseen;
193 : #ifdef INFINITES_ALLOWED
194 : float infs;
195 : #else
196 : int infs;
197 : #endif
198 : double *partials;
199 : } *pergroup;
200 1108 : BUN listi;
201 1108 : int parti;
202 1108 : int i;
203 1108 : BUN grp;
204 1108 : double x, y;
205 1108 : volatile double lo, hi;
206 1108 : double twopow = pow((double) FLT_RADIX, (double) (DBL_MAX_EXP - 1));
207 1108 : BUN nils = 0;
208 1108 : volatile flt f;
209 :
210 1108 : lng timeoffset = 0;
211 1108 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
212 1108 : if (qry_ctx != NULL) {
213 1108 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
214 : }
215 :
216 : /* we only deal with the two floating point types */
217 1108 : assert(tp1 == TYPE_flt || tp1 == TYPE_dbl);
218 1108 : assert(tp2 == TYPE_flt || tp2 == TYPE_dbl);
219 : /* if input is dbl, then so it output */
220 1108 : assert(tp2 == TYPE_flt || tp1 == TYPE_dbl);
221 : /* if no gids, then we have a single group */
222 1108 : assert(ngrp == 1 || gids != NULL);
223 1108 : if (gids == NULL || ngrp == 1) {
224 1092 : min = max = 0;
225 1092 : ngrp = 1;
226 1092 : gids = NULL;
227 : }
228 1108 : pergroup = GDKmalloc(ngrp * sizeof(*pergroup));
229 1108 : if (pergroup == NULL)
230 : return BUN_NONE;
231 10737 : for (grp = 0; grp < ngrp; grp++) {
232 19258 : pergroup[grp] = (struct pergroup) {
233 : .maxpartials = 2,
234 9629 : .partials = GDKmalloc(2 * sizeof(double)),
235 : };
236 9629 : if (pergroup[grp].partials == NULL) {
237 0 : while (grp > 0)
238 0 : GDKfree(pergroup[--grp].partials);
239 0 : GDKfree(pergroup);
240 0 : return BUN_NONE;
241 : }
242 : }
243 23688 : TIMEOUT_LOOP(ci->ncand, timeoffset) {
244 21472 : listi = canditer_next(ci) - seqb;
245 21472 : grp = gids ? gids[listi] : 0;
246 21472 : if (grp < min || grp > max)
247 0 : continue;
248 21472 : if (pergroup[grp].partials == NULL)
249 0 : continue;
250 21472 : if (tp1 == TYPE_flt && !is_flt_nil(((const flt *) values)[listi]))
251 16 : x = ((const flt *) values)[listi];
252 21456 : else if (tp1 == TYPE_dbl && !is_dbl_nil(((const dbl *) values)[listi]))
253 : x = ((const dbl *) values)[listi];
254 : else {
255 : /* it's a nil */
256 78 : if (!skip_nils) {
257 0 : if (tp2 == TYPE_flt)
258 0 : ((flt *) results)[grp] = flt_nil;
259 : else
260 0 : ((dbl *) results)[grp] = dbl_nil;
261 0 : GDKfree(pergroup[grp].partials);
262 0 : pergroup[grp].partials = NULL;
263 0 : if (++nils == ngrp)
264 0 : TIMEOUT_LOOP_BREAK;
265 : }
266 78 : continue;
267 : }
268 21394 : pergroup[grp].valseen = true;
269 : #ifdef INFINITES_ALLOWED
270 : if (isinf(x)) {
271 : pergroup[grp].infs += x;
272 : continue;
273 : }
274 : #endif
275 21394 : i = 0;
276 41164 : for (parti = 0; parti < pergroup[grp].npartials; parti++) {
277 19770 : y = pergroup[grp].partials[parti];
278 19770 : if (fabs(x) < fabs(y))
279 9438 : exchange(&x, &y);
280 19770 : twosum(&hi, &lo, x, y);
281 19770 : if (isinf(hi)) {
282 58 : int sign = hi > 0 ? 1 : -1;
283 58 : hi = x - twopow * sign;
284 58 : x = hi - twopow * sign;
285 58 : pergroup[grp].infs += sign;
286 58 : if (fabs(x) < fabs(y))
287 8 : exchange(&x, &y);
288 58 : twosum(&hi, &lo, x, y);
289 : }
290 19770 : if (lo != 0)
291 8891 : pergroup[grp].partials[i++] = lo;
292 19770 : x = hi;
293 : }
294 21394 : if (x != 0) {
295 21323 : if (i == pergroup[grp].maxpartials) {
296 247 : double *temp;
297 247 : pergroup[grp].maxpartials += pergroup[grp].maxpartials;
298 247 : temp = GDKrealloc(pergroup[grp].partials, pergroup[grp].maxpartials * sizeof(double));
299 247 : if (temp == NULL)
300 0 : goto bailout;
301 247 : pergroup[grp].partials = temp;
302 : }
303 21323 : pergroup[grp].partials[i++] = x;
304 : }
305 21394 : pergroup[grp].npartials = i;
306 : }
307 1108 : TIMEOUT_CHECK(timeoffset, GOTO_LABEL_TIMEOUT_HANDLER(bailout));
308 10709 : for (grp = 0; grp < ngrp; grp++) {
309 9629 : if (pergroup[grp].partials == NULL)
310 0 : continue;
311 9629 : if (!pergroup[grp].valseen) {
312 31 : if (tp2 == TYPE_flt)
313 5 : ((flt *) results)[grp] = nil_if_empty ? flt_nil : 0;
314 : else
315 26 : ((dbl *) results)[grp] = nil_if_empty ? dbl_nil : 0;
316 31 : nils += nil_if_empty;
317 31 : GDKfree(pergroup[grp].partials);
318 31 : pergroup[grp].partials = NULL;
319 31 : continue;
320 : }
321 : #ifdef INFINITES_ALLOWED
322 : if (isinf(pergroup[grp].infs) || isnan(pergroup[grp].infs)) {
323 : goto overflow;
324 : }
325 : #endif
326 :
327 9598 : if ((pergroup[grp].infs == 1 || pergroup[grp].infs == -1) &&
328 42 : pergroup[grp].npartials > 0 &&
329 38 : !samesign(pergroup[grp].infs, pergroup[grp].partials[pergroup[grp].npartials - 1])) {
330 34 : twosum(&hi, &lo, pergroup[grp].infs * twopow, pergroup[grp].partials[pergroup[grp].npartials - 1] / 2);
331 34 : if (isinf(2 * hi)) {
332 22 : y = 2 * lo;
333 22 : x = hi + y;
334 22 : x -= hi;
335 22 : if (x == y &&
336 18 : pergroup[grp].npartials > 1 &&
337 12 : samesign(lo, pergroup[grp].partials[pergroup[grp].npartials - 2])) {
338 6 : GDKfree(pergroup[grp].partials);
339 6 : pergroup[grp].partials = NULL;
340 6 : x = 2 * (hi + y);
341 6 : if (tp2 == TYPE_flt) {
342 0 : f = (flt) x;
343 0 : if (isinf(f) ||
344 0 : isnan(f) ||
345 0 : is_flt_nil(f)) {
346 0 : goto overflow;
347 : } else {
348 0 : ((flt *) results)[grp] = f;
349 : }
350 6 : } else if (is_dbl_nil(x)) {
351 0 : goto overflow;
352 : } else
353 6 : ((dbl *) results)[grp] = x;
354 6 : continue;
355 : }
356 : } else {
357 12 : if (lo) {
358 2 : if (pergroup[grp].npartials == pergroup[grp].maxpartials) {
359 0 : double *temp;
360 : /* we need space for one more */
361 0 : pergroup[grp].maxpartials++;
362 0 : temp = GDKrealloc(pergroup[grp].partials, pergroup[grp].maxpartials * sizeof(double));
363 0 : if (temp == NULL)
364 0 : goto bailout;
365 0 : pergroup[grp].partials = temp;
366 : }
367 2 : pergroup[grp].partials[pergroup[grp].npartials - 1] = 2 * lo;
368 2 : pergroup[grp].partials[pergroup[grp].npartials++] = 2 * hi;
369 : } else {
370 10 : pergroup[grp].partials[pergroup[grp].npartials - 1] = 2 * hi;
371 : }
372 12 : pergroup[grp].infs = 0;
373 : }
374 : }
375 :
376 9592 : if (pergroup[grp].infs != 0)
377 28 : goto overflow;
378 :
379 9564 : if (pergroup[grp].npartials == 0) {
380 14 : GDKfree(pergroup[grp].partials);
381 14 : pergroup[grp].partials = NULL;
382 14 : if (tp2 == TYPE_flt)
383 0 : ((flt *) results)[grp] = 0;
384 : else
385 14 : ((dbl *) results)[grp] = 0;
386 14 : continue;
387 : }
388 :
389 : /* accumulate into hi */
390 9550 : hi = pergroup[grp].partials[--pergroup[grp].npartials];
391 9550 : while (pergroup[grp].npartials > 0) {
392 615 : twosum(&hi, &lo, hi, pergroup[grp].partials[--pergroup[grp].npartials]);
393 615 : if (lo) {
394 615 : pergroup[grp].partials[pergroup[grp].npartials++] = lo;
395 615 : break;
396 : }
397 : }
398 :
399 9550 : if (pergroup[grp].npartials >= 2 &&
400 243 : samesign(pergroup[grp].partials[pergroup[grp].npartials - 1], pergroup[grp].partials[pergroup[grp].npartials - 2]) &&
401 25 : hi + 2 * pergroup[grp].partials[pergroup[grp].npartials - 1] - hi == 2 * pergroup[grp].partials[pergroup[grp].npartials - 1]) {
402 20 : hi += 2 * pergroup[grp].partials[pergroup[grp].npartials - 1];
403 20 : pergroup[grp].partials[pergroup[grp].npartials - 1] = -pergroup[grp].partials[pergroup[grp].npartials - 1];
404 : }
405 :
406 9550 : GDKfree(pergroup[grp].partials);
407 9550 : pergroup[grp].partials = NULL;
408 9550 : if (tp2 == TYPE_flt) {
409 11 : f = (flt) hi;
410 11 : if (isinf(f) || isnan(f) || is_flt_nil(f)) {
411 0 : goto overflow;
412 : } else {
413 11 : ((flt *) results)[grp] = f;
414 : }
415 9539 : } else if (is_dbl_nil(hi)) {
416 0 : goto overflow;
417 : } else {
418 9539 : ((dbl *) results)[grp] = hi;
419 : }
420 : }
421 1080 : GDKfree(pergroup);
422 1080 : return nils;
423 :
424 28 : overflow:
425 28 : GDKerror("22003!overflow in sum aggregate.\n");
426 : bailout:
427 56 : for (grp = 0; grp < ngrp; grp++)
428 28 : GDKfree(pergroup[grp].partials);
429 28 : GDKfree(pergroup);
430 28 : return BUN_NONE;
431 : }
432 :
433 : #define AGGR_SUM(TYPE1, TYPE2) \
434 : do { \
435 : TYPE1 x; \
436 : const TYPE1 *restrict vals = (const TYPE1 *) values; \
437 : if (ngrp == 1 && ci->tpe == cand_dense) { \
438 : /* single group, no candidate list */ \
439 : TYPE2 sum; \
440 : *algo = "sum: no candidates, no groups"; \
441 : sum = 0; \
442 : if (nonil) { \
443 : *seen = ci->ncand > 0; \
444 : TIMEOUT_LOOP_IDX(i, ci->ncand, timeoffset) { \
445 : x = vals[ci->seq + i - seqb]; \
446 : ADD_WITH_CHECK(x, sum, \
447 : TYPE2, sum, \
448 : GDK_##TYPE2##_max, \
449 : goto overflow); \
450 : } \
451 : TIMEOUT_CHECK(timeoffset, \
452 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
453 : } else { \
454 : bool seenval = false; \
455 : TIMEOUT_LOOP_IDX(i, ci->ncand, timeoffset) { \
456 : x = vals[ci->seq + i - seqb]; \
457 : if (is_##TYPE1##_nil(x)) { \
458 : if (!skip_nils) { \
459 : sum = TYPE2##_nil; \
460 : nils = 1; \
461 : TIMEOUT_LOOP_BREAK; \
462 : } \
463 : } else { \
464 : ADD_WITH_CHECK(x, sum, \
465 : TYPE2, sum, \
466 : GDK_##TYPE2##_max, \
467 : goto overflow); \
468 : seenval = true; \
469 : } \
470 : } \
471 : TIMEOUT_CHECK(timeoffset, \
472 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
473 : *seen = seenval; \
474 : } \
475 : if (*seen) \
476 : *sums = sum; \
477 : } else if (ngrp == 1) { \
478 : /* single group, with candidate list */ \
479 : TYPE2 sum; \
480 : bool seenval = false; \
481 : *algo = "sum: with candidates, no groups"; \
482 : sum = 0; \
483 : TIMEOUT_LOOP_IDX(i, ci->ncand, timeoffset) { \
484 : x = vals[canditer_next(ci) - seqb]; \
485 : if (is_##TYPE1##_nil(x)) { \
486 : if (!skip_nils) { \
487 : sum = TYPE2##_nil; \
488 : nils = 1; \
489 : TIMEOUT_LOOP_BREAK; \
490 : } \
491 : } else { \
492 : ADD_WITH_CHECK(x, sum, \
493 : TYPE2, sum, \
494 : GDK_##TYPE2##_max, \
495 : goto overflow); \
496 : seenval = true; \
497 : } \
498 : } \
499 : TIMEOUT_CHECK(timeoffset, \
500 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
501 : if (seenval) \
502 : *sums = sum; \
503 : } else if (ci->tpe == cand_dense) { \
504 : /* multiple groups, no candidate list */ \
505 : *algo = "sum: no candidates, with groups"; \
506 : TIMEOUT_LOOP_IDX(i, ci->ncand, timeoffset) { \
507 : if (gids == NULL || \
508 : (gids[i] >= min && gids[i] <= max)) { \
509 : gid = gids ? gids[i] - min : (oid) i; \
510 : x = vals[ci->seq + i - seqb]; \
511 : if (is_##TYPE1##_nil(x)) { \
512 : if (!skip_nils) { \
513 : sums[gid] = TYPE2##_nil; \
514 : nils++; \
515 : } \
516 : } else { \
517 : if (nil_if_empty && \
518 : !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
519 : seen[gid >> 5] |= 1U << (gid & 0x1F); \
520 : sums[gid] = 0; \
521 : } \
522 : if (!is_##TYPE2##_nil(sums[gid])) { \
523 : ADD_WITH_CHECK( \
524 : x, \
525 : sums[gid], \
526 : TYPE2, \
527 : sums[gid], \
528 : GDK_##TYPE2##_max, \
529 : goto overflow); \
530 : } \
531 : } \
532 : } \
533 : } \
534 : TIMEOUT_CHECK(timeoffset, \
535 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
536 : } else { \
537 : /* multiple groups, with candidate list */ \
538 : *algo = "sum: with candidates, with groups"; \
539 : TIMEOUT_LOOP(ci->ncand, timeoffset) { \
540 : i = canditer_next(ci) - seqb; \
541 : if (gids == NULL || \
542 : (gids[i] >= min && gids[i] <= max)) { \
543 : gid = gids ? gids[i] - min : (oid) i; \
544 : x = vals[i]; \
545 : if (is_##TYPE1##_nil(x)) { \
546 : if (!skip_nils) { \
547 : sums[gid] = TYPE2##_nil; \
548 : nils++; \
549 : } \
550 : } else { \
551 : if (nil_if_empty && \
552 : !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
553 : seen[gid >> 5] |= 1U << (gid & 0x1F); \
554 : sums[gid] = 0; \
555 : } \
556 : if (!is_##TYPE2##_nil(sums[gid])) { \
557 : ADD_WITH_CHECK( \
558 : x, \
559 : sums[gid], \
560 : TYPE2, \
561 : sums[gid], \
562 : GDK_##TYPE2##_max, \
563 : goto overflow); \
564 : } \
565 : } \
566 : } \
567 : } \
568 : TIMEOUT_CHECK(timeoffset, \
569 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
570 : } \
571 : } while (0)
572 :
573 : #define AGGR_SUM_NOOVL(TYPE1, TYPE2) \
574 : do { \
575 : TYPE1 x; \
576 : const TYPE1 *restrict vals = (const TYPE1 *) values; \
577 : if (ngrp == 1 && ci->tpe == cand_dense) { \
578 : /* single group, no candidate list */ \
579 : TYPE2 sum; \
580 : sum = 0; \
581 : if (nonil) { \
582 : *algo = "sum: no candidates, no groups, no nils, no overflow"; \
583 : *seen = ci->ncand > 0; \
584 : TIMEOUT_LOOP_IDX(i, ci->ncand, timeoffset) { \
585 : sum += vals[ci->seq + i - seqb]; \
586 : } \
587 : TIMEOUT_CHECK(timeoffset, \
588 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
589 : } else { \
590 : bool seenval = false; \
591 : *algo = "sum: no candidates, no groups, no overflow"; \
592 : TIMEOUT_LOOP_IDX(i, ci->ncand, timeoffset) { \
593 : x = vals[ci->seq + i - seqb]; \
594 : if (is_##TYPE1##_nil(x)) { \
595 : if (!skip_nils) { \
596 : sum = TYPE2##_nil; \
597 : nils = 1; \
598 : TIMEOUT_LOOP_BREAK; \
599 : } \
600 : } else { \
601 : sum += x; \
602 : seenval = true; \
603 : } \
604 : } \
605 : TIMEOUT_CHECK(timeoffset, \
606 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
607 : *seen = seenval; \
608 : } \
609 : if (*seen) \
610 : *sums = sum; \
611 : } else if (ngrp == 1) { \
612 : /* single group, with candidate list */ \
613 : TYPE2 sum; \
614 : bool seenval = false; \
615 : *algo = "sum: with candidates, no groups, no overflow"; \
616 : sum = 0; \
617 : TIMEOUT_LOOP_IDX(i, ci->ncand, timeoffset) { \
618 : x = vals[canditer_next(ci) - seqb]; \
619 : if (is_##TYPE1##_nil(x)) { \
620 : if (!skip_nils) { \
621 : sum = TYPE2##_nil; \
622 : nils = 1; \
623 : TIMEOUT_LOOP_BREAK; \
624 : } \
625 : } else { \
626 : sum += x; \
627 : seenval = true; \
628 : } \
629 : } \
630 : TIMEOUT_CHECK(timeoffset, \
631 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
632 : if (seenval) \
633 : *sums = sum; \
634 : } else if (ci->tpe == cand_dense) { \
635 : /* multiple groups, no candidate list */ \
636 : if (nonil) { \
637 : *algo = "sum: no candidates, with groups, no nils, no overflow"; \
638 : TIMEOUT_LOOP_IDX(i, ci->ncand, timeoffset) { \
639 : if (gids == NULL || \
640 : (gids[i] >= min && gids[i] <= max)) { \
641 : gid = gids ? gids[i] - min : (oid) i; \
642 : x = vals[ci->seq + i - seqb]; \
643 : if (nil_if_empty && \
644 : !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
645 : seen[gid >> 5] |= 1U << (gid & 0x1F); \
646 : sums[gid] = 0; \
647 : } \
648 : sums[gid] += x; \
649 : } \
650 : } \
651 : TIMEOUT_CHECK(timeoffset, \
652 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
653 : } else { \
654 : *algo = "sum: no candidates, with groups, no overflow"; \
655 : TIMEOUT_LOOP_IDX(i, ci->ncand, timeoffset) { \
656 : if (gids == NULL || \
657 : (gids[i] >= min && gids[i] <= max)) { \
658 : gid = gids ? gids[i] - min : (oid) i; \
659 : x = vals[ci->seq + i - seqb]; \
660 : if (is_##TYPE1##_nil(x)) { \
661 : if (!skip_nils) { \
662 : sums[gid] = TYPE2##_nil; \
663 : nils++; \
664 : } \
665 : } else { \
666 : if (nil_if_empty && \
667 : !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
668 : seen[gid >> 5] |= 1U << (gid & 0x1F); \
669 : sums[gid] = 0; \
670 : } \
671 : if (!is_##TYPE2##_nil(sums[gid])) { \
672 : sums[gid] += x; \
673 : } \
674 : } \
675 : } \
676 : } \
677 : TIMEOUT_CHECK(timeoffset, \
678 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
679 : } \
680 : } else { \
681 : /* multiple groups, with candidate list */ \
682 : *algo = "sum: with candidates, with groups, no overflow"; \
683 : TIMEOUT_LOOP(ci->ncand, timeoffset) { \
684 : i = canditer_next(ci) - seqb; \
685 : if (gids == NULL || \
686 : (gids[i] >= min && gids[i] <= max)) { \
687 : gid = gids ? gids[i] - min : (oid) i; \
688 : x = vals[i]; \
689 : if (is_##TYPE1##_nil(x)) { \
690 : if (!skip_nils) { \
691 : sums[gid] = TYPE2##_nil; \
692 : nils++; \
693 : } \
694 : } else { \
695 : if (nil_if_empty && \
696 : !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
697 : seen[gid >> 5] |= 1U << (gid & 0x1F); \
698 : sums[gid] = 0; \
699 : } \
700 : if (!is_##TYPE2##_nil(sums[gid])) { \
701 : sums[gid] += x; \
702 : } \
703 : } \
704 : } \
705 : } \
706 : TIMEOUT_CHECK(timeoffset, \
707 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
708 : } \
709 : } while (0)
710 :
711 : static BUN
712 10666 : dosum(const void *restrict values, bool nonil, oid seqb,
713 : struct canditer *restrict ci,
714 : void *restrict results, BUN ngrp, int tp1, int tp2,
715 : const oid *restrict gids,
716 : oid min, oid max, bool skip_nils,
717 : bool nil_if_empty, const char *func, const char **algo)
718 : {
719 10666 : BUN nils = 0;
720 10666 : BUN i;
721 10666 : oid gid;
722 10666 : unsigned int *restrict seen = NULL; /* bitmask for groups that we've seen */
723 :
724 10666 : lng timeoffset = 0;
725 10666 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
726 10665 : if (qry_ctx != NULL) {
727 6892 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
728 : }
729 :
730 10665 : switch (tp2) {
731 12 : case TYPE_flt:
732 12 : if (tp1 != TYPE_flt)
733 0 : goto unsupported;
734 : /* fall through */
735 : case TYPE_dbl:
736 1091 : if (tp1 != TYPE_flt && tp1 != TYPE_dbl)
737 0 : goto unsupported;
738 1103 : *algo = "sum: floating point";
739 1103 : return dofsum(values, seqb, ci, results, ngrp, tp1, tp2,
740 : gids, min, max, skip_nils,
741 : nil_if_empty);
742 : }
743 :
744 : /* allocate bitmap for seen group ids */
745 9562 : seen = GDKzalloc(((ngrp + 31) / 32) * sizeof(int));
746 9563 : if (seen == NULL) {
747 : return BUN_NONE;
748 : }
749 :
750 9563 : switch (tp2) {
751 0 : case TYPE_bte: {
752 0 : bte *restrict sums = (bte *) results;
753 0 : switch (tp1) {
754 0 : case TYPE_bte:
755 0 : AGGR_SUM(bte, bte);
756 : break;
757 0 : default:
758 0 : goto unsupported;
759 : }
760 : break;
761 : }
762 0 : case TYPE_sht: {
763 0 : sht *restrict sums = (sht *) results;
764 0 : switch (tp1) {
765 0 : case TYPE_bte:
766 0 : if (ci->ncand < ((BUN) 1 << ((sizeof(sht) - sizeof(bte)) << 3)))
767 0 : AGGR_SUM_NOOVL(bte, sht);
768 : else
769 0 : AGGR_SUM(bte, sht);
770 : break;
771 0 : case TYPE_sht:
772 0 : AGGR_SUM(sht, sht);
773 : break;
774 0 : default:
775 0 : goto unsupported;
776 : }
777 : break;
778 : }
779 12 : case TYPE_int: {
780 12 : int *restrict sums = (int *) results;
781 12 : switch (tp1) {
782 0 : case TYPE_bte:
783 0 : if (ci->ncand < ((BUN) 1 << ((sizeof(int) - sizeof(bte)) << 3)))
784 0 : AGGR_SUM_NOOVL(bte, int);
785 : else
786 0 : AGGR_SUM(bte, int);
787 : break;
788 0 : case TYPE_sht:
789 0 : if (ci->ncand < ((BUN) 1 << ((sizeof(int) - sizeof(sht)) << 3)))
790 0 : AGGR_SUM_NOOVL(sht, int);
791 : else
792 0 : AGGR_SUM(sht, int);
793 : break;
794 12 : case TYPE_int:
795 54 : AGGR_SUM(int, int);
796 : break;
797 0 : default:
798 0 : goto unsupported;
799 : }
800 : break;
801 : }
802 7466 : case TYPE_lng: {
803 7466 : lng *restrict sums = (lng *) results;
804 7466 : switch (tp1) {
805 : #if SIZEOF_BUN == 4
806 : case TYPE_bte:
807 : AGGR_SUM_NOOVL(bte, lng);
808 : break;
809 : case TYPE_sht:
810 : AGGR_SUM_NOOVL(sht, lng);
811 : break;
812 : case TYPE_int:
813 : AGGR_SUM_NOOVL(int, lng);
814 : break;
815 : #else
816 0 : case TYPE_bte:
817 0 : if (ci->ncand < ((BUN) 1 << ((sizeof(lng) - sizeof(bte)) << 3)))
818 0 : AGGR_SUM_NOOVL(bte, lng);
819 : else
820 0 : AGGR_SUM(bte, lng);
821 : break;
822 0 : case TYPE_sht:
823 0 : if (ci->ncand < ((BUN) 1 << ((sizeof(lng) - sizeof(sht)) << 3)))
824 0 : AGGR_SUM_NOOVL(sht, lng);
825 : else
826 0 : AGGR_SUM(sht, lng);
827 : break;
828 0 : case TYPE_int:
829 0 : if (ci->ncand < ((BUN) 1 << ((sizeof(lng) - sizeof(int)) << 3)))
830 0 : AGGR_SUM_NOOVL(int, lng);
831 : else
832 0 : AGGR_SUM(int, lng);
833 : break;
834 : #endif
835 7466 : case TYPE_lng:
836 8470408 : AGGR_SUM(lng, lng);
837 : break;
838 0 : default:
839 0 : goto unsupported;
840 : }
841 : break;
842 : }
843 : #ifdef HAVE_HGE
844 2085 : case TYPE_hge: {
845 2085 : hge *sums = (hge *) results;
846 2085 : switch (tp1) {
847 115 : case TYPE_bte:
848 3834346 : AGGR_SUM_NOOVL(bte, hge);
849 : break;
850 10 : case TYPE_sht:
851 51 : AGGR_SUM_NOOVL(sht, hge);
852 : break;
853 1251 : case TYPE_int:
854 65674235 : AGGR_SUM_NOOVL(int, hge);
855 : break;
856 20 : case TYPE_lng:
857 164262 : AGGR_SUM_NOOVL(lng, hge);
858 : break;
859 689 : case TYPE_hge:
860 33754868 : AGGR_SUM(hge, hge);
861 : break;
862 0 : default:
863 0 : goto unsupported;
864 : }
865 : break;
866 : }
867 : #endif
868 0 : default:
869 0 : goto unsupported;
870 : }
871 :
872 9559 : if (nils == 0 && nil_if_empty) {
873 : /* figure out whether there were any empty groups
874 : * (that result in a nil value) */
875 9563 : if (ngrp & 0x1F) {
876 : /* fill last slot */
877 9546 : seen[ngrp >> 5] |= ~0U << (ngrp & 0x1F);
878 : }
879 226860 : for (i = 0, ngrp = (ngrp + 31) / 32; i < ngrp; i++) {
880 217665 : if (seen[i] != ~0U) {
881 : nils = 1;
882 : break;
883 : }
884 : }
885 : }
886 9559 : GDKfree(seen);
887 :
888 9559 : return nils;
889 :
890 0 : unsupported:
891 0 : GDKfree(seen);
892 0 : GDKerror("%s: type combination (sum(%s)->%s) not supported.\n",
893 : func, ATOMname(tp1), ATOMname(tp2));
894 0 : return BUN_NONE;
895 :
896 0 : overflow:
897 0 : GDKfree(seen);
898 0 : GDKerror("22003!overflow in sum aggregate.\n");
899 0 : return BUN_NONE;
900 :
901 0 : bailout:
902 0 : GDKfree(seen);
903 0 : return BUN_NONE;
904 : }
905 :
906 : /* calculate group sums with optional candidates list */
907 : BAT *
908 5443 : BATgroupsum(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
909 : {
910 5443 : const oid *restrict gids;
911 5443 : oid min, max;
912 5443 : BUN ngrp;
913 5443 : BUN nils;
914 5443 : BAT *bn;
915 5443 : struct canditer ci;
916 5443 : const char *err;
917 5443 : const char *algo = NULL;
918 5443 : lng t0 = 0;
919 :
920 5443 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
921 :
922 5443 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
923 0 : GDKerror("%s\n", err);
924 0 : return NULL;
925 : }
926 5441 : if (g == NULL) {
927 0 : GDKerror("b and g must be aligned\n");
928 0 : return NULL;
929 : }
930 :
931 5441 : if (ci.ncand == 0 || ngrp == 0) {
932 : /* trivial: no sums, so return bat aligned with g with
933 : * nil in the tail */
934 1594 : return BATconstant(ngrp == 0 ? 0 : min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
935 : }
936 :
937 3847 : if ((e == NULL ||
938 3847 : (BATcount(e) == ci.ncand && e->hseqbase == ci.hseq)) &&
939 1353 : (BATtdense(g) || (g->tkey && g->tnonil))) {
940 : /* trivial: singleton groups, so all results are equal
941 : * to the inputs (but possibly a different type) */
942 1353 : return BATconvert(b, s, tp, 0, 0, 0);
943 : }
944 :
945 2494 : bn = BATconstant(min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
946 2495 : if (bn == NULL) {
947 : return NULL;
948 : }
949 :
950 2495 : if (BATtdense(g))
951 : gids = NULL;
952 : else
953 2268 : gids = (const oid *) Tloc(g, 0);
954 :
955 2495 : BATiter bi = bat_iterator(b);
956 4990 : nils = dosum(bi.base, bi.nonil, b->hseqbase, &ci,
957 2495 : Tloc(bn, 0), ngrp, bi.type, tp, gids, min, max,
958 : skip_nils, true, __func__, &algo);
959 2495 : bat_iterator_end(&bi);
960 :
961 2495 : if (nils < BUN_NONE) {
962 2495 : BATsetcount(bn, ngrp);
963 2495 : bn->tkey = BATcount(bn) <= 1;
964 2495 : bn->tsorted = BATcount(bn) <= 1;
965 2495 : bn->trevsorted = BATcount(bn) <= 1;
966 2495 : bn->tnil = nils != 0;
967 2495 : bn->tnonil = nils == 0;
968 : } else {
969 0 : BBPunfix(bn->batCacheid);
970 0 : bn = NULL;
971 : }
972 :
973 2495 : if (algo)
974 2495 : MT_thread_setalgorithm(algo);
975 2495 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
976 : "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
977 : "; start " OIDFMT ", count " BUNFMT " (%s -- " LLFMT " usec)\n",
978 : ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
979 : ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
980 : ci.seq, ci.ncand, algo ? algo : "", GDKusec() - t0);
981 : return bn;
982 : }
983 :
984 : static BUN
985 30 : mskCountOnes(BAT *b, struct canditer *ci)
986 : {
987 30 : BUN cnt = 0, ncand = ci->ncand;
988 :
989 30 : if (ci->s == NULL && mask_cand(b))
990 30 : return BATcount(b);
991 0 : if (ci->tpe == cand_dense && ncand > 0 && !mask_cand(b)) {
992 0 : BATiter bi = bat_iterator(b);
993 0 : const uint32_t *restrict src = (const uint32_t *) bi.base + (ci->seq - b->hseqbase) / 32;
994 0 : int bits = (ci->seq - b->hseqbase) % 32;
995 0 : if (bits + ncand <= 32) {
996 0 : if (ncand == 32)
997 0 : cnt = candmask_pop(src[0]);
998 : else
999 0 : cnt = candmask_pop(src[0] & (((1U << ncand) - 1) << bits));
1000 0 : bat_iterator_end(&bi);
1001 0 : return cnt;
1002 : }
1003 0 : if (bits != 0) {
1004 0 : cnt = candmask_pop(src[0] & (~0U << bits));
1005 0 : src++;
1006 0 : ncand -= 32 - bits;
1007 : }
1008 0 : while (ncand >= 32) {
1009 0 : cnt += candmask_pop(*src);
1010 0 : src++;
1011 0 : ncand -= 32;
1012 : }
1013 0 : if (ncand > 0)
1014 0 : cnt += candmask_pop(*src & ((1U << ncand) - 1));
1015 0 : bat_iterator_end(&bi);
1016 0 : return cnt;
1017 : }
1018 0 : for (BUN i = 0; i < ncand; i++) {
1019 0 : BUN x = canditer_next(ci) - b->hseqbase;
1020 0 : cnt += mskGetVal(b, x);
1021 : }
1022 : return cnt;
1023 : }
1024 :
1025 : gdk_return
1026 8459 : BATsum(void *res, int tp, BAT *b, BAT *s, bool skip_nils, bool nil_if_empty)
1027 : {
1028 8459 : oid min, max;
1029 8459 : BUN ngrp;
1030 8459 : struct canditer ci;
1031 8459 : const char *err;
1032 8459 : const char *algo = NULL;
1033 8459 : lng t0 = 0;
1034 :
1035 8459 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
1036 :
1037 8459 : if ((err = BATgroupaggrinit(b, NULL, NULL, s, &min, &max, &ngrp, &ci)) != NULL) {
1038 0 : GDKerror("%s\n", err);
1039 0 : return GDK_FAIL;
1040 : }
1041 8459 : if (ATOMstorage(b->ttype) == TYPE_msk || mask_cand(b)) {
1042 30 : BUN n = mskCountOnes(b, &ci);
1043 30 : switch (tp) {
1044 0 : case TYPE_bte:
1045 0 : if (n > GDK_bte_max) {
1046 0 : GDKerror("22003!overflow in sum aggregate.\n");
1047 0 : return GDK_FAIL;
1048 : }
1049 0 : * (bte *) res = (bte) n;
1050 0 : break;
1051 0 : case TYPE_sht:
1052 0 : if (n > GDK_sht_max) {
1053 0 : GDKerror("22003!overflow in sum aggregate.\n");
1054 0 : return GDK_FAIL;
1055 : }
1056 0 : * (sht *) res = (sht) n;
1057 0 : break;
1058 0 : case TYPE_int:
1059 : #if SIZEOF_BUN > 4
1060 0 : if (n > GDK_int_max) {
1061 0 : GDKerror("22003!overflow in sum aggregate.\n");
1062 0 : return GDK_FAIL;
1063 : }
1064 : #endif
1065 0 : * (int *) res = (int) n;
1066 0 : break;
1067 30 : case TYPE_lng:
1068 30 : * (lng *) res = (lng) n;
1069 30 : break;
1070 : #ifdef HAVE_HGE
1071 0 : case TYPE_hge:
1072 0 : * (hge *) res = (hge) n;
1073 0 : break;
1074 : #endif
1075 0 : case TYPE_flt:
1076 0 : * (flt *) res = (flt) n;
1077 0 : break;
1078 0 : case TYPE_dbl:
1079 0 : * (dbl *) res = (dbl) n;
1080 0 : break;
1081 0 : default:
1082 0 : GDKerror("type combination (sum(%s)->%s) not supported.\n",
1083 : ATOMname(b->ttype), ATOMname(tp));
1084 0 : return GDK_FAIL;
1085 : }
1086 30 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",s=" ALGOOPTBATFMT "; "
1087 : "start " OIDFMT ", count " BUNFMT " (pop count -- " LLFMT " usec)\n",
1088 : ALGOBATPAR(b), ALGOOPTBATPAR(s),
1089 : ci.seq, ci.ncand, GDKusec() - t0);
1090 30 : return GDK_SUCCEED;
1091 : }
1092 8429 : switch (tp) {
1093 0 : case TYPE_bte:
1094 0 : * (bte *) res = nil_if_empty ? bte_nil : 0;
1095 0 : break;
1096 0 : case TYPE_sht:
1097 0 : * (sht *) res = nil_if_empty ? sht_nil : 0;
1098 0 : break;
1099 16 : case TYPE_int:
1100 16 : * (int *) res = nil_if_empty ? int_nil : 0;
1101 16 : break;
1102 6370 : case TYPE_lng:
1103 6370 : * (lng *) res = nil_if_empty ? lng_nil : 0;
1104 6370 : break;
1105 : #ifdef HAVE_HGE
1106 934 : case TYPE_hge:
1107 934 : * (hge *) res = nil_if_empty ? hge_nil : 0;
1108 934 : break;
1109 : #endif
1110 1109 : case TYPE_flt:
1111 : case TYPE_dbl:
1112 1109 : switch (b->ttype) {
1113 1 : case TYPE_bte:
1114 : case TYPE_sht:
1115 : case TYPE_int:
1116 : case TYPE_lng:
1117 : #ifdef HAVE_HGE
1118 : case TYPE_hge:
1119 : #endif
1120 : {
1121 : /* special case for summing integer types into
1122 : * a floating point: We calculate the average
1123 : * (which is done exactly), and multiply the
1124 : * result by the count to get the sum. Note
1125 : * that if we just summed into a floating
1126 : * point number, we could loose too much
1127 : * accuracy, and if we summed into lng first,
1128 : * we could get unnecessary overflow. */
1129 1 : dbl avg;
1130 1 : BUN cnt;
1131 :
1132 1 : if (BATcalcavg(b, s, &avg, &cnt, 0) != GDK_SUCCEED)
1133 : return GDK_FAIL;
1134 1 : if (cnt == 0) {
1135 0 : avg = nil_if_empty ? dbl_nil : 0;
1136 : }
1137 1 : if (cnt < ci.ncand && !skip_nils) {
1138 : /* there were nils */
1139 0 : avg = dbl_nil;
1140 : }
1141 1 : if (tp == TYPE_flt) {
1142 0 : if (is_dbl_nil(avg))
1143 0 : *(flt *) res = flt_nil;
1144 0 : else if (cnt > 0 &&
1145 0 : GDK_flt_max / cnt < fabs(avg)) {
1146 0 : GDKerror("22003!overflow in sum aggregate.\n");
1147 0 : return GDK_FAIL;
1148 : } else {
1149 0 : *(flt *) res = (flt) avg * cnt;
1150 : }
1151 : } else {
1152 1 : if (is_dbl_nil(avg)) {
1153 0 : *(dbl *) res = dbl_nil;
1154 1 : } else if (cnt > 0 &&
1155 1 : GDK_dbl_max / cnt < fabs(avg)) {
1156 0 : GDKerror("22003!overflow in sum aggregate.\n");
1157 0 : return GDK_FAIL;
1158 : } else {
1159 1 : *(dbl *) res = avg * cnt;
1160 : }
1161 : }
1162 : return GDK_SUCCEED;
1163 : }
1164 : default:
1165 1108 : break;
1166 : }
1167 1108 : if (b->ttype == TYPE_dbl)
1168 1095 : * (dbl *) res = nil_if_empty ? dbl_nil : 0;
1169 : else
1170 13 : * (flt *) res = nil_if_empty ? flt_nil : 0;
1171 : break;
1172 0 : default:
1173 0 : GDKerror("type combination (sum(%s)->%s) not supported.\n",
1174 : ATOMname(b->ttype), ATOMname(tp));
1175 0 : return GDK_FAIL;
1176 : }
1177 8428 : if (ci.ncand == 0)
1178 : return GDK_SUCCEED;
1179 8171 : BATiter bi = bat_iterator(b);
1180 16342 : BUN nils = dosum(bi.base, bi.nonil, b->hseqbase, &ci,
1181 8171 : res, true, bi.type, tp, &min, min, max,
1182 : skip_nils, nil_if_empty, __func__, &algo);
1183 8171 : bat_iterator_end(&bi);
1184 8169 : if (algo)
1185 8169 : MT_thread_setalgorithm(algo);
1186 8171 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",s=" ALGOOPTBATFMT "; "
1187 : "start " OIDFMT ", count " BUNFMT " (%s -- " LLFMT " usec)\n",
1188 : ALGOBATPAR(b), ALGOOPTBATPAR(s),
1189 : ci.seq, ci.ncand, algo ? algo : "", GDKusec() - t0);
1190 8171 : return nils < BUN_NONE ? GDK_SUCCEED : GDK_FAIL;
1191 : }
1192 :
1193 : /* ---------------------------------------------------------------------- */
1194 : /* product */
1195 :
1196 : #define AGGR_PROD(TYPE1, TYPE2, TYPE3) \
1197 : do { \
1198 : const TYPE1 *restrict vals = (const TYPE1 *) values; \
1199 : gid = 0; /* doesn't change if gidincr == false */ \
1200 : TIMEOUT_LOOP(ci->ncand, timeoffset) { \
1201 : i = canditer_next(ci) - seqb; \
1202 : if (gids == NULL || !gidincr || \
1203 : (gids[i] >= min && gids[i] <= max)) { \
1204 : if (gidincr) { \
1205 : if (gids) \
1206 : gid = gids[i] - min; \
1207 : else \
1208 : gid = (oid) i; \
1209 : } \
1210 : if (is_##TYPE1##_nil(vals[i])) { \
1211 : if (!skip_nils) { \
1212 : prods[gid] = TYPE2##_nil; \
1213 : nils++; \
1214 : } \
1215 : } else { \
1216 : if (nil_if_empty && \
1217 : !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
1218 : seen[gid >> 5] |= 1U << (gid & 0x1F); \
1219 : prods[gid] = 1; \
1220 : } \
1221 : if (!is_##TYPE2##_nil(prods[gid])) { \
1222 : MUL4_WITH_CHECK( \
1223 : vals[i], \
1224 : prods[gid], \
1225 : TYPE2, prods[gid], \
1226 : GDK_##TYPE2##_max, \
1227 : TYPE3, \
1228 : goto overflow); \
1229 : } \
1230 : } \
1231 : } \
1232 : } \
1233 : TIMEOUT_CHECK(timeoffset, \
1234 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
1235 : } while (0)
1236 :
1237 : #ifdef HAVE_HGE
1238 : #define AGGR_PROD_HGE(TYPE) \
1239 : do { \
1240 : const TYPE *vals = (const TYPE *) values; \
1241 : gid = 0; /* doesn't change if gidincr == false */ \
1242 : TIMEOUT_LOOP(ci->ncand, timeoffset) { \
1243 : i = canditer_next(ci) - seqb; \
1244 : if (gids == NULL || !gidincr || \
1245 : (gids[i] >= min && gids[i] <= max)) { \
1246 : if (gidincr) { \
1247 : if (gids) \
1248 : gid = gids[i] - min; \
1249 : else \
1250 : gid = (oid) i; \
1251 : } \
1252 : if (nil_if_empty && \
1253 : !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
1254 : seen[gid >> 5] |= 1U << (gid & 0x1F); \
1255 : prods[gid] = 1; \
1256 : } \
1257 : if (is_##TYPE##_nil(vals[i])) { \
1258 : if (!skip_nils) { \
1259 : prods[gid] = hge_nil; \
1260 : nils++; \
1261 : } \
1262 : } else if (!is_hge_nil(prods[gid])) { \
1263 : HGEMUL_CHECK(vals[i], \
1264 : prods[gid], \
1265 : prods[gid], \
1266 : GDK_hge_max, \
1267 : goto overflow); \
1268 : } \
1269 : } \
1270 : } \
1271 : TIMEOUT_CHECK(timeoffset, \
1272 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
1273 : } while (0)
1274 : #else
1275 : #define AGGR_PROD_LNG(TYPE) \
1276 : do { \
1277 : const TYPE *restrict vals = (const TYPE *) values; \
1278 : gid = 0; /* doesn't change if gidincr == false */ \
1279 : TIMEOUT_LOOP(ci->ncand, timeoffset) { \
1280 : i = canditer_next(ci) - seqb; \
1281 : if (gids == NULL || !gidincr || \
1282 : (gids[i] >= min && gids[i] <= max)) { \
1283 : if (gidincr) { \
1284 : if (gids) \
1285 : gid = gids[i] - min; \
1286 : else \
1287 : gid = (oid) i; \
1288 : } \
1289 : if (is_##TYPE##_nil(vals[i])) { \
1290 : if (!skip_nils) { \
1291 : prods[gid] = lng_nil; \
1292 : nils++; \
1293 : } \
1294 : } else { \
1295 : if (nil_if_empty && \
1296 : !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
1297 : seen[gid >> 5] |= 1U << (gid & 0x1F); \
1298 : prods[gid] = 1; \
1299 : } \
1300 : if (!is_lng_nil(prods[gid])) { \
1301 : LNGMUL_CHECK( \
1302 : vals[i], \
1303 : prods[gid], \
1304 : prods[gid], \
1305 : GDK_lng_max, \
1306 : goto overflow); \
1307 : } \
1308 : } \
1309 : } \
1310 : } \
1311 : TIMEOUT_CHECK(timeoffset, \
1312 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
1313 : } while (0)
1314 : #endif
1315 :
1316 : #define AGGR_PROD_FLOAT(TYPE1, TYPE2) \
1317 : do { \
1318 : const TYPE1 *restrict vals = (const TYPE1 *) values; \
1319 : gid = 0; /* doesn't change if gidincr == false */ \
1320 : TIMEOUT_LOOP(ci->ncand, timeoffset) { \
1321 : i = canditer_next(ci) - seqb; \
1322 : if (gids == NULL || !gidincr || \
1323 : (gids[i] >= min && gids[i] <= max)) { \
1324 : if (gidincr) { \
1325 : if (gids) \
1326 : gid = gids[i] - min; \
1327 : else \
1328 : gid = (oid) i; \
1329 : } \
1330 : if (is_##TYPE1##_nil(vals[i])) { \
1331 : if (!skip_nils) { \
1332 : prods[gid] = TYPE2##_nil; \
1333 : nils++; \
1334 : } \
1335 : } else { \
1336 : if (nil_if_empty && \
1337 : !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
1338 : seen[gid >> 5] |= 1U << (gid & 0x1F); \
1339 : prods[gid] = 1; \
1340 : } \
1341 : if (!is_##TYPE2##_nil(prods[gid])) { \
1342 : if (ABSOLUTE(vals[i]) > 1 && \
1343 : GDK_##TYPE2##_max / ABSOLUTE(vals[i]) < ABSOLUTE(prods[gid])) { \
1344 : goto overflow; \
1345 : } else { \
1346 : prods[gid] *= vals[i]; \
1347 : } \
1348 : } \
1349 : } \
1350 : } \
1351 : } \
1352 : TIMEOUT_CHECK(timeoffset, \
1353 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
1354 : } while (0)
1355 :
1356 : static BUN
1357 32 : doprod(const void *restrict values, oid seqb, struct canditer *restrict ci,
1358 : void *restrict results, BUN ngrp, int tp1, int tp2,
1359 : const oid *restrict gids, bool gidincr, oid min, oid max,
1360 : bool skip_nils, bool nil_if_empty, const char *func)
1361 : {
1362 32 : BUN nils = 0;
1363 32 : BUN i;
1364 32 : oid gid;
1365 32 : unsigned int *restrict seen; /* bitmask for groups that we've seen */
1366 :
1367 32 : lng timeoffset = 0;
1368 32 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
1369 32 : if (qry_ctx != NULL) {
1370 32 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
1371 : }
1372 :
1373 : /* allocate bitmap for seen group ids */
1374 32 : seen = GDKzalloc(((ngrp + 31) / 32) * sizeof(int));
1375 32 : if (seen == NULL) {
1376 : return BUN_NONE;
1377 : }
1378 :
1379 32 : switch (tp2) {
1380 0 : case TYPE_bte: {
1381 0 : bte *restrict prods = (bte *) results;
1382 0 : switch (tp1) {
1383 0 : case TYPE_bte:
1384 0 : AGGR_PROD(bte, bte, sht);
1385 : break;
1386 0 : default:
1387 0 : goto unsupported;
1388 : }
1389 : break;
1390 : }
1391 0 : case TYPE_sht: {
1392 0 : sht *restrict prods = (sht *) results;
1393 0 : switch (tp1) {
1394 0 : case TYPE_bte:
1395 0 : AGGR_PROD(bte, sht, int);
1396 : break;
1397 0 : case TYPE_sht:
1398 0 : AGGR_PROD(sht, sht, int);
1399 : break;
1400 0 : default:
1401 0 : goto unsupported;
1402 : }
1403 : break;
1404 : }
1405 0 : case TYPE_int: {
1406 0 : int *restrict prods = (int *) results;
1407 0 : switch (tp1) {
1408 0 : case TYPE_bte:
1409 0 : AGGR_PROD(bte, int, lng);
1410 : break;
1411 0 : case TYPE_sht:
1412 0 : AGGR_PROD(sht, int, lng);
1413 : break;
1414 0 : case TYPE_int:
1415 0 : AGGR_PROD(int, int, lng);
1416 : break;
1417 0 : default:
1418 0 : goto unsupported;
1419 : }
1420 : break;
1421 : }
1422 : #ifdef HAVE_HGE
1423 0 : case TYPE_lng: {
1424 0 : lng *prods = (lng *) results;
1425 0 : switch (tp1) {
1426 0 : case TYPE_bte:
1427 0 : AGGR_PROD(bte, lng, hge);
1428 : break;
1429 0 : case TYPE_sht:
1430 0 : AGGR_PROD(sht, lng, hge);
1431 : break;
1432 0 : case TYPE_int:
1433 0 : AGGR_PROD(int, lng, hge);
1434 : break;
1435 0 : case TYPE_lng:
1436 0 : AGGR_PROD(lng, lng, hge);
1437 : break;
1438 0 : default:
1439 0 : goto unsupported;
1440 : }
1441 : break;
1442 : }
1443 25 : case TYPE_hge: {
1444 25 : hge *prods = (hge *) results;
1445 25 : switch (tp1) {
1446 4 : case TYPE_bte:
1447 14 : AGGR_PROD_HGE(bte);
1448 : break;
1449 4 : case TYPE_sht:
1450 14 : AGGR_PROD_HGE(sht);
1451 : break;
1452 4 : case TYPE_int:
1453 14 : AGGR_PROD_HGE(int);
1454 : break;
1455 8 : case TYPE_lng:
1456 28 : AGGR_PROD_HGE(lng);
1457 : break;
1458 5 : case TYPE_hge:
1459 30 : AGGR_PROD_HGE(hge);
1460 : break;
1461 0 : default:
1462 0 : goto unsupported;
1463 : }
1464 : break;
1465 : }
1466 : #else
1467 : case TYPE_lng: {
1468 : lng *restrict prods = (lng *) results;
1469 : switch (tp1) {
1470 : case TYPE_bte:
1471 : AGGR_PROD_LNG(bte);
1472 : break;
1473 : case TYPE_sht:
1474 : AGGR_PROD_LNG(sht);
1475 : break;
1476 : case TYPE_int:
1477 : AGGR_PROD_LNG(int);
1478 : break;
1479 : case TYPE_lng:
1480 : AGGR_PROD_LNG(lng);
1481 : break;
1482 : default:
1483 : goto unsupported;
1484 : }
1485 : break;
1486 : }
1487 : #endif
1488 0 : case TYPE_flt: {
1489 0 : flt *restrict prods = (flt *) results;
1490 0 : switch (tp1) {
1491 0 : case TYPE_bte:
1492 0 : AGGR_PROD_FLOAT(bte, flt);
1493 : break;
1494 0 : case TYPE_sht:
1495 0 : AGGR_PROD_FLOAT(sht, flt);
1496 : break;
1497 0 : case TYPE_int:
1498 0 : AGGR_PROD_FLOAT(int, flt);
1499 : break;
1500 0 : case TYPE_lng:
1501 0 : AGGR_PROD_FLOAT(lng, flt);
1502 : break;
1503 : #ifdef HAVE_HGE
1504 0 : case TYPE_hge:
1505 0 : AGGR_PROD_FLOAT(hge, flt);
1506 : break;
1507 : #endif
1508 0 : case TYPE_flt:
1509 0 : AGGR_PROD_FLOAT(flt, flt);
1510 : break;
1511 0 : default:
1512 0 : goto unsupported;
1513 : }
1514 : break;
1515 : }
1516 7 : case TYPE_dbl: {
1517 7 : dbl *restrict prods = (dbl *) results;
1518 7 : switch (tp1) {
1519 0 : case TYPE_bte:
1520 0 : AGGR_PROD_FLOAT(bte, dbl);
1521 : break;
1522 0 : case TYPE_sht:
1523 0 : AGGR_PROD_FLOAT(sht, dbl);
1524 : break;
1525 0 : case TYPE_int:
1526 0 : AGGR_PROD_FLOAT(int, dbl);
1527 : break;
1528 0 : case TYPE_lng:
1529 0 : AGGR_PROD_FLOAT(lng, dbl);
1530 : break;
1531 : #ifdef HAVE_HGE
1532 0 : case TYPE_hge:
1533 0 : AGGR_PROD_FLOAT(hge, dbl);
1534 : break;
1535 : #endif
1536 0 : case TYPE_flt:
1537 0 : AGGR_PROD_FLOAT(flt, dbl);
1538 : break;
1539 7 : case TYPE_dbl:
1540 35 : AGGR_PROD_FLOAT(dbl, dbl);
1541 : break;
1542 0 : default:
1543 0 : goto unsupported;
1544 : }
1545 : break;
1546 : }
1547 0 : default:
1548 0 : goto unsupported;
1549 : }
1550 :
1551 32 : if (nils == 0 && nil_if_empty) {
1552 : /* figure out whether there were any empty groups
1553 : * (that result in a nil value) */
1554 32 : if (ngrp & 0x1F) {
1555 : /* fill last slot */
1556 32 : seen[ngrp >> 5] |= ~0U << (ngrp & 0x1F);
1557 : }
1558 64 : for (i = 0, ngrp = (ngrp + 31) / 32; i < ngrp; i++) {
1559 32 : if (seen[i] != ~0U) {
1560 : nils = 1;
1561 : break;
1562 : }
1563 : }
1564 : }
1565 32 : GDKfree(seen);
1566 :
1567 32 : return nils;
1568 :
1569 0 : unsupported:
1570 0 : GDKfree(seen);
1571 0 : GDKerror("%s: type combination (mul(%s)->%s) not supported.\n",
1572 : func, ATOMname(tp1), ATOMname(tp2));
1573 0 : return BUN_NONE;
1574 :
1575 0 : overflow:
1576 0 : GDKfree(seen);
1577 0 : GDKerror("22003!overflow in product aggregate.\n");
1578 0 : return BUN_NONE;
1579 :
1580 0 : bailout:
1581 0 : GDKfree(seen);
1582 0 : return BUN_NONE;
1583 : }
1584 :
1585 : /* calculate group products with optional candidates list */
1586 : BAT *
1587 17 : BATgroupprod(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
1588 : {
1589 17 : const oid *restrict gids;
1590 17 : oid min, max;
1591 17 : BUN ngrp;
1592 17 : BUN nils;
1593 17 : BAT *bn;
1594 17 : struct canditer ci;
1595 17 : const char *err;
1596 17 : lng t0 = 0;
1597 :
1598 17 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
1599 :
1600 17 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
1601 0 : GDKerror("%s\n", err);
1602 0 : return NULL;
1603 : }
1604 17 : if (g == NULL) {
1605 0 : GDKerror("b and g must be aligned\n");
1606 0 : return NULL;
1607 : }
1608 :
1609 17 : if (ci.ncand == 0 || ngrp == 0) {
1610 : /* trivial: no products, so return bat aligned with g
1611 : * with nil in the tail */
1612 6 : return BATconstant(ngrp == 0 ? 0 : min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
1613 : }
1614 :
1615 11 : if ((e == NULL ||
1616 11 : (BATcount(e) == ci.ncand && e->hseqbase == ci.hseq)) &&
1617 9 : (BATtdense(g) || (g->tkey && g->tnonil))) {
1618 : /* trivial: singleton groups, so all results are equal
1619 : * to the inputs (but possibly a different type) */
1620 9 : return BATconvert(b, s, tp, 0, 0, 0);
1621 : }
1622 :
1623 2 : bn = BATconstant(min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
1624 2 : if (bn == NULL) {
1625 : return NULL;
1626 : }
1627 :
1628 2 : if (BATtdense(g))
1629 : gids = NULL;
1630 : else
1631 2 : gids = (const oid *) Tloc(g, 0);
1632 :
1633 2 : BATiter bi = bat_iterator(b);
1634 4 : nils = doprod(bi.base, b->hseqbase, &ci, Tloc(bn, 0), ngrp,
1635 2 : bi.type, tp, gids, true, min, max, skip_nils,
1636 : true, __func__);
1637 2 : bat_iterator_end(&bi);
1638 :
1639 2 : if (nils < BUN_NONE) {
1640 2 : BATsetcount(bn, ngrp);
1641 2 : bn->tkey = BATcount(bn) <= 1;
1642 2 : bn->tsorted = BATcount(bn) <= 1;
1643 2 : bn->trevsorted = BATcount(bn) <= 1;
1644 2 : bn->tnil = nils != 0;
1645 2 : bn->tnonil = nils == 0;
1646 : } else {
1647 0 : BBPunfix(bn->batCacheid);
1648 0 : bn = NULL;
1649 : }
1650 :
1651 2 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
1652 : "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
1653 : "; start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
1654 : ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
1655 : ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
1656 : ci.seq, ci.ncand, GDKusec() - t0);
1657 :
1658 : return bn;
1659 : }
1660 :
1661 : gdk_return
1662 30 : BATprod(void *res, int tp, BAT *b, BAT *s, bool skip_nils, bool nil_if_empty)
1663 : {
1664 30 : oid min, max;
1665 30 : BUN ngrp;
1666 30 : BUN nils;
1667 30 : struct canditer ci;
1668 30 : const char *err;
1669 30 : lng t0 = 0;
1670 :
1671 30 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
1672 :
1673 30 : if ((err = BATgroupaggrinit(b, NULL, NULL, s, &min, &max, &ngrp, &ci)) != NULL) {
1674 0 : GDKerror("%s\n", err);
1675 0 : return GDK_FAIL;
1676 : }
1677 30 : switch (tp) {
1678 0 : case TYPE_bte:
1679 0 : * (bte *) res = nil_if_empty ? bte_nil : (bte) 1;
1680 0 : break;
1681 0 : case TYPE_sht:
1682 0 : * (sht *) res = nil_if_empty ? sht_nil : (sht) 1;
1683 0 : break;
1684 0 : case TYPE_int:
1685 0 : * (int *) res = nil_if_empty ? int_nil : (int) 1;
1686 0 : break;
1687 0 : case TYPE_lng:
1688 0 : * (lng *) res = nil_if_empty ? lng_nil : (lng) 1;
1689 0 : break;
1690 : #ifdef HAVE_HGE
1691 25 : case TYPE_hge:
1692 25 : * (hge *) res = nil_if_empty ? hge_nil : (hge) 1;
1693 25 : break;
1694 : #endif
1695 0 : case TYPE_flt:
1696 0 : * (flt *) res = nil_if_empty ? flt_nil : (flt) 1;
1697 0 : break;
1698 5 : case TYPE_dbl:
1699 5 : * (dbl *) res = nil_if_empty ? dbl_nil : (dbl) 1;
1700 5 : break;
1701 0 : default:
1702 0 : GDKerror("type combination (prod(%s)->%s) not supported.\n",
1703 : ATOMname(b->ttype), ATOMname(tp));
1704 0 : return GDK_FAIL;
1705 : }
1706 30 : if (ci.ncand == 0)
1707 : return GDK_SUCCEED;
1708 30 : BATiter bi = bat_iterator(b);
1709 60 : nils = doprod(bi.base, b->hseqbase, &ci, res, true,
1710 30 : bi.type, tp, &min, false, min, max,
1711 : skip_nils, nil_if_empty, __func__);
1712 30 : bat_iterator_end(&bi);
1713 30 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",s=" ALGOOPTBATFMT "; "
1714 : "start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
1715 : ALGOBATPAR(b), ALGOOPTBATPAR(s),
1716 : ci.seq, ci.ncand, GDKusec() - t0);
1717 30 : return nils < BUN_NONE ? GDK_SUCCEED : GDK_FAIL;
1718 : }
1719 :
1720 : /* ---------------------------------------------------------------------- */
1721 : /* average */
1722 :
1723 : #define GOTO_BAILOUT() \
1724 : do { \
1725 : GDKfree(avgs); \
1726 : GOTO_LABEL_TIMEOUT_HANDLER(bailout); \
1727 : } while (0)
1728 :
1729 : #define AGGR_AVG(TYPE) \
1730 : do { \
1731 : const TYPE *restrict vals = (const TYPE *) bi.base; \
1732 : TYPE *restrict avgs = GDKzalloc(ngrp * sizeof(TYPE)); \
1733 : if (avgs == NULL) \
1734 : goto bailout; \
1735 : TIMEOUT_LOOP(ci.ncand, timeoffset) { \
1736 : i = canditer_next(&ci) - b->hseqbase; \
1737 : if (gids == NULL || \
1738 : (gids[i] >= min && gids[i] <= max)) { \
1739 : if (gids) \
1740 : gid = gids[i] - min; \
1741 : else \
1742 : gid = (oid) i; \
1743 : if (is_##TYPE##_nil(vals[i])) { \
1744 : if (!skip_nils) \
1745 : cnts[gid] = lng_nil; \
1746 : } else if (!is_lng_nil(cnts[gid])) { \
1747 : AVERAGE_ITER(TYPE, vals[i], \
1748 : avgs[gid], \
1749 : rems[gid], \
1750 : cnts[gid]); \
1751 : } \
1752 : } \
1753 : } \
1754 : TIMEOUT_CHECK(timeoffset, GOTO_BAILOUT()); \
1755 : for (i = 0; i < ngrp; i++) { \
1756 : if (cnts[i] == 0 || is_lng_nil(cnts[i])) { \
1757 : dbls[i] = dbl_nil; \
1758 : cnts[i] = 0; \
1759 : nils++; \
1760 : } else { \
1761 : dbls[i] = avgs[i] + (dbl) rems[i] / cnts[i]; \
1762 : } \
1763 : } \
1764 : GDKfree(avgs); \
1765 : } while (0)
1766 :
1767 : #define AGGR_AVG_FLOAT(TYPE) \
1768 : do { \
1769 : const TYPE *restrict vals = (const TYPE *) bi.base; \
1770 : for (i = 0; i < ngrp; i++) \
1771 : dbls[i] = 0; \
1772 : TIMEOUT_LOOP(ci.ncand, timeoffset) { \
1773 : i = canditer_next(&ci) - b->hseqbase; \
1774 : if (gids == NULL || \
1775 : (gids[i] >= min && gids[i] <= max)) { \
1776 : if (gids) \
1777 : gid = gids[i] - min; \
1778 : else \
1779 : gid = (oid) i; \
1780 : if (is_##TYPE##_nil(vals[i])) { \
1781 : if (!skip_nils) \
1782 : cnts[gid] = lng_nil; \
1783 : } else if (!is_lng_nil(cnts[gid])) { \
1784 : AVERAGE_ITER_FLOAT(TYPE, vals[i], \
1785 : dbls[gid], \
1786 : cnts[gid]); \
1787 : } \
1788 : } \
1789 : } \
1790 : TIMEOUT_CHECK(timeoffset, \
1791 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
1792 : for (i = 0; i < ngrp; i++) { \
1793 : if (cnts[i] == 0 || is_lng_nil(cnts[i])) { \
1794 : dbls[i] = dbl_nil; \
1795 : cnts[i] = 0; \
1796 : nils++; \
1797 : } \
1798 : } \
1799 : } while (0)
1800 :
1801 : /* There are three functions that are used for calculating averages.
1802 : * The first one (BATgroupavg) returns averages as a floating point
1803 : * value, the other two (BATgroupavg3 and BATgroupavg3combine) work
1804 : * together to return averages in the domain type (which should be an
1805 : * integer type). */
1806 :
1807 : /* Calculate group averages with optional candidates list. The average
1808 : * that is calculated is returned in a dbl, independent of the type of
1809 : * the input. The average is calculated exactly, so not in a floating
1810 : * point which could potentially losse bits during processing
1811 : * (e.g. average of 2**62 and a billion 1's). */
1812 : gdk_return
1813 240 : BATgroupavg(BAT **bnp, BAT **cntsp, BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, int scale)
1814 : {
1815 240 : const oid *restrict gids;
1816 240 : oid gid;
1817 240 : oid min, max;
1818 240 : BUN i, ngrp;
1819 240 : BUN nils = 0;
1820 240 : lng *restrict rems = NULL;
1821 240 : lng *restrict cnts = NULL;
1822 240 : dbl *restrict dbls;
1823 240 : BAT *bn = NULL, *cn = NULL;
1824 240 : struct canditer ci;
1825 240 : const char *err;
1826 240 : lng t0 = 0;
1827 240 : BATiter bi = {0};
1828 :
1829 240 : lng timeoffset = 0;
1830 240 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
1831 240 : if (qry_ctx != NULL) {
1832 240 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
1833 : }
1834 :
1835 240 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
1836 :
1837 240 : assert(tp == TYPE_dbl);
1838 240 : (void) tp; /* compatibility (with other BATgroup*
1839 : * functions) argument */
1840 :
1841 240 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
1842 0 : GDKerror("%s\n", err);
1843 0 : return GDK_FAIL;
1844 : }
1845 240 : if (g == NULL) {
1846 0 : GDKerror("b and g must be aligned\n");
1847 0 : return GDK_FAIL;
1848 : }
1849 :
1850 240 : if (ci.ncand == 0 || ngrp == 0) {
1851 : /* trivial: no averages, so return bat aligned with g
1852 : * with nil in the tail */
1853 7 : bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
1854 7 : if (bn == NULL) {
1855 : return GDK_FAIL;
1856 : }
1857 7 : if (cntsp) {
1858 0 : lng zero = 0;
1859 0 : if ((cn = BATconstant(ngrp == 0 ? 0 : min, TYPE_lng, &zero, ngrp, TRANSIENT)) == NULL) {
1860 0 : BBPreclaim(bn);
1861 0 : return GDK_FAIL;
1862 : }
1863 0 : *cntsp = cn;
1864 : }
1865 7 : *bnp = bn;
1866 7 : return GDK_SUCCEED;
1867 : }
1868 :
1869 233 : if ((!skip_nils || cntsp == NULL || b->tnonil) &&
1870 164 : (e == NULL ||
1871 164 : (BATcount(e) == ci.ncand && e->hseqbase == b->hseqbase)) &&
1872 112 : (BATtdense(g) || (g->tkey && g->tnonil))) {
1873 : /* trivial: singleton groups, so all results are equal
1874 : * to the inputs (but possibly a different type) */
1875 112 : if ((bn = BATconvert(b, s, TYPE_dbl, 0, 0, 0)) == NULL)
1876 : return GDK_FAIL;
1877 112 : if (cntsp) {
1878 6 : lng one = 1;
1879 6 : if ((cn = BATconstant(ngrp == 0 ? 0 : min, TYPE_lng, &one, ngrp, TRANSIENT)) == NULL) {
1880 0 : BBPreclaim(bn);
1881 0 : return GDK_FAIL;
1882 : }
1883 6 : *cntsp = cn;
1884 : }
1885 112 : *bnp = bn;
1886 112 : return GDK_SUCCEED;
1887 : }
1888 :
1889 : /* allocate temporary space to do per group calculations */
1890 121 : switch (b->ttype) {
1891 106 : case TYPE_bte:
1892 : case TYPE_sht:
1893 : case TYPE_int:
1894 : case TYPE_lng:
1895 : #ifdef HAVE_HGE
1896 : case TYPE_hge:
1897 : #endif
1898 106 : rems = GDKzalloc(ngrp * sizeof(lng));
1899 106 : if (rems == NULL)
1900 0 : goto bailout1;
1901 : break;
1902 : default:
1903 : break;
1904 : }
1905 121 : if (cntsp) {
1906 99 : if ((cn = COLnew(min, TYPE_lng, ngrp, TRANSIENT)) == NULL)
1907 0 : goto bailout1;
1908 99 : cnts = (lng *) Tloc(cn, 0);
1909 99 : memset(cnts, 0, ngrp * sizeof(lng));
1910 : } else {
1911 22 : cnts = GDKzalloc(ngrp * sizeof(lng));
1912 22 : if (cnts == NULL)
1913 0 : goto bailout1;
1914 : }
1915 :
1916 121 : bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
1917 121 : if (bn == NULL)
1918 0 : goto bailout1;
1919 121 : dbls = (dbl *) Tloc(bn, 0);
1920 :
1921 121 : if (BATtdense(g))
1922 : gids = NULL;
1923 : else
1924 84 : gids = (const oid *) Tloc(g, 0);
1925 :
1926 121 : bi = bat_iterator(b);
1927 121 : switch (b->ttype) {
1928 0 : case TYPE_bte:
1929 0 : AGGR_AVG(bte);
1930 0 : break;
1931 4 : case TYPE_sht:
1932 17 : AGGR_AVG(sht);
1933 4 : break;
1934 88 : case TYPE_int:
1935 21847445 : AGGR_AVG(int);
1936 88 : break;
1937 14 : case TYPE_lng:
1938 107 : AGGR_AVG(lng);
1939 14 : break;
1940 : #ifdef HAVE_HGE
1941 0 : case TYPE_hge:
1942 0 : AGGR_AVG(hge);
1943 0 : break;
1944 : #endif
1945 6 : case TYPE_flt:
1946 2965 : AGGR_AVG_FLOAT(flt);
1947 : break;
1948 9 : case TYPE_dbl:
1949 298 : AGGR_AVG_FLOAT(dbl);
1950 : break;
1951 0 : default:
1952 0 : GDKerror("type (%s) not supported.\n", ATOMname(bi.type));
1953 0 : goto bailout;
1954 : }
1955 121 : bat_iterator_end(&bi);
1956 121 : GDKfree(rems);
1957 121 : if (cn == NULL)
1958 22 : GDKfree(cnts);
1959 : else {
1960 99 : BATsetcount(cn, ngrp);
1961 99 : cn->tkey = BATcount(cn) <= 1;
1962 99 : cn->tsorted = BATcount(cn) <= 1;
1963 99 : cn->trevsorted = BATcount(cn) <= 1;
1964 99 : cn->tnil = false;
1965 99 : cn->tnonil = true;
1966 99 : *cntsp = cn;
1967 : }
1968 121 : if (scale != 0) {
1969 0 : dbl fac = pow(10.0, (double) scale);
1970 0 : for (i = 0; i < ngrp; i++) {
1971 0 : if (!is_dbl_nil(dbls[i]))
1972 0 : dbls[i] /= fac;
1973 : }
1974 : }
1975 121 : BATsetcount(bn, ngrp);
1976 121 : bn->tkey = BATcount(bn) <= 1;
1977 121 : bn->tsorted = BATcount(bn) <= 1;
1978 121 : bn->trevsorted = BATcount(bn) <= 1;
1979 121 : bn->tnil = nils != 0;
1980 121 : bn->tnonil = nils == 0;
1981 121 : *bnp = bn;
1982 121 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
1983 : "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
1984 : "; start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
1985 : ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
1986 : ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
1987 : ci.seq, ci.ncand, GDKusec() - t0);
1988 : return GDK_SUCCEED;
1989 0 : bailout:
1990 0 : bat_iterator_end(&bi);
1991 0 : bailout1:
1992 0 : BBPreclaim(bn);
1993 0 : GDKfree(rems);
1994 0 : if (cntsp) {
1995 0 : BBPreclaim(*cntsp);
1996 0 : *cntsp = NULL;
1997 0 : } else if (cnts) {
1998 0 : GDKfree(cnts);
1999 : }
2000 : return GDK_FAIL;
2001 : }
2002 :
2003 : /* An exact numeric average of a bunch of values consists of three
2004 : * parts: the average rounded down (towards minus infinity), the number
2005 : * of values that participated in the calculation, and the remainder.
2006 : * The remainder is in the range 0 (inclusive) to count (not inclusive).
2007 : * BATgroupavg3 calculates these values for each given group. The
2008 : * function below, BATgroupavg3combine, combines averages calculated
2009 : * this way to correct averages by rounding or truncating towards zero
2010 : * (depending on the symbol TRUNCATE_NUMBERS). */
2011 : gdk_return
2012 270 : BATgroupavg3(BAT **avgp, BAT **remp, BAT **cntp, BAT *b, BAT *g, BAT *e, BAT *s, bool skip_nils)
2013 : {
2014 270 : const char *err;
2015 270 : oid min, max;
2016 270 : BUN ngrp;
2017 270 : struct canditer ci;
2018 270 : BAT *bn, *rn, *cn;
2019 270 : BUN i;
2020 270 : oid o;
2021 :
2022 270 : lng timeoffset = 0;
2023 270 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
2024 270 : if (qry_ctx != NULL) {
2025 270 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
2026 : }
2027 :
2028 270 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
2029 0 : GDKerror("%s\n", err);
2030 0 : return GDK_FAIL;
2031 : }
2032 270 : if (ci.ncand == 0 || ngrp == 0) {
2033 18 : if (ngrp == 0)
2034 12 : min = 0;
2035 18 : bn = BATconstant(min, b->ttype, ATOMnilptr(b->ttype),
2036 : ngrp, TRANSIENT);
2037 18 : rn = BATconstant(min, TYPE_lng, &lng_nil, ngrp, TRANSIENT);
2038 18 : cn = BATconstant(min, TYPE_lng, &(lng){0}, ngrp, TRANSIENT);
2039 18 : if (bn == NULL || rn == NULL || cn == NULL) {
2040 0 : BBPreclaim(bn);
2041 0 : BBPreclaim(rn);
2042 0 : BBPreclaim(cn);
2043 0 : return GDK_FAIL;
2044 : }
2045 18 : *avgp = bn;
2046 18 : *remp = rn;
2047 18 : *cntp = cn;
2048 18 : return GDK_SUCCEED;
2049 : }
2050 252 : ValRecord zero;
2051 252 : (void) VALinit(&zero, TYPE_bte, &(bte){0});
2052 252 : bn = BATconstant(min, b->ttype, VALconvert(b->ttype, &zero),
2053 : ngrp, TRANSIENT);
2054 252 : rn = BATconstant(min, TYPE_lng, &(lng){0}, ngrp, TRANSIENT);
2055 252 : cn = BATconstant(min, TYPE_lng, &(lng){0}, ngrp, TRANSIENT);
2056 252 : if (bn == NULL || rn == NULL || cn == NULL) {
2057 0 : BBPreclaim(bn);
2058 0 : BBPreclaim(rn);
2059 0 : BBPreclaim(cn);
2060 0 : return GDK_FAIL;
2061 : }
2062 252 : lng *rems = Tloc(rn, 0);
2063 252 : lng *cnts = Tloc(cn, 0);
2064 252 : const oid *gids = g && !BATtdense(g) ? Tloc(g, 0) : NULL;
2065 252 : oid gid = ngrp == 1 && gids ? gids[0] - min : 0;
2066 :
2067 252 : BATiter bi = bat_iterator(b);
2068 :
2069 504 : switch (ATOMbasetype(b->ttype)) {
2070 8 : case TYPE_bte: {
2071 8 : const bte *vals = (const bte *) bi.base;
2072 8 : bte *avgs = Tloc(bn, 0);
2073 257 : TIMEOUT_LOOP(ci.ncand, timeoffset) {
2074 241 : o = canditer_next(&ci) - b->hseqbase;
2075 241 : if (ngrp > 1)
2076 0 : gid = gids ? gids[o] - min : o;
2077 241 : if (is_bte_nil(vals[o])) {
2078 0 : if (!skip_nils) {
2079 0 : avgs[gid] = bte_nil;
2080 0 : rems[gid] = lng_nil;
2081 0 : cnts[gid] = lng_nil;
2082 0 : bn->tnil = true;
2083 0 : rn->tnil = true;
2084 0 : cn->tnil = true;
2085 : }
2086 241 : } else if (!is_lng_nil(cnts[gid])) {
2087 241 : AVERAGE_ITER(bte, vals[o], avgs[gid], rems[gid], cnts[gid]);
2088 : }
2089 : }
2090 32 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset) {
2091 8 : if (cnts[i] == 0) {
2092 0 : avgs[i] = bte_nil;
2093 0 : bn->tnil = true;
2094 : } else
2095 : #ifdef TRUNCATE_NUMBERS
2096 : if (!is_bte_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0) {
2097 : avgs[i]++;
2098 : rems[i] -= cnts[i];
2099 : }
2100 : #else
2101 8 : if (!is_bte_nil(avgs[i]) && rems[i] > 0) {
2102 6 : if (avgs[i] < 0) {
2103 0 : if (2*rems[i] > cnts[i]) {
2104 0 : avgs[i]++;
2105 0 : rems[i] -= cnts[i];
2106 : }
2107 : } else {
2108 6 : if (2*rems[i] >= cnts[i]) {
2109 2 : avgs[i]++;
2110 2 : rems[i] -= cnts[i];
2111 : }
2112 : }
2113 : }
2114 : #endif
2115 : }
2116 : break;
2117 : }
2118 1 : case TYPE_sht: {
2119 1 : const sht *vals = (const sht *) bi.base;
2120 1 : sht *avgs = Tloc(bn, 0);
2121 21 : TIMEOUT_LOOP(ci.ncand, timeoffset) {
2122 19 : o = canditer_next(&ci) - b->hseqbase;
2123 19 : if (ngrp > 1)
2124 0 : gid = gids ? gids[o] - min : o;
2125 19 : if (is_sht_nil(vals[o])) {
2126 0 : if (!skip_nils) {
2127 0 : avgs[gid] = sht_nil;
2128 0 : rems[gid] = lng_nil;
2129 0 : cnts[gid] = lng_nil;
2130 0 : bn->tnil = true;
2131 0 : rn->tnil = true;
2132 0 : cn->tnil = true;
2133 : }
2134 19 : } else if (!is_lng_nil(cnts[gid])) {
2135 19 : AVERAGE_ITER(sht, vals[o], avgs[gid], rems[gid], cnts[gid]);
2136 : }
2137 : }
2138 4 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset) {
2139 1 : if (cnts[i] == 0) {
2140 0 : avgs[i] = sht_nil;
2141 0 : bn->tnil = true;
2142 : } else
2143 : #ifdef TRUNCATE_NUMBERS
2144 : if (!is_sht_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0) {
2145 : avgs[i]++;
2146 : rems[i] -= cnts[i];
2147 : }
2148 : #else
2149 1 : if (!is_sht_nil(avgs[i]) && rems[i] > 0) {
2150 1 : if (avgs[i] < 0) {
2151 0 : if (2*rems[i] > cnts[i]) {
2152 0 : avgs[i]++;
2153 0 : rems[i] -= cnts[i];
2154 : }
2155 : } else {
2156 1 : if (2*rems[i] >= cnts[i]) {
2157 0 : avgs[i]++;
2158 0 : rems[i] -= cnts[i];
2159 : }
2160 : }
2161 : }
2162 : #endif
2163 : }
2164 : break;
2165 : }
2166 140 : case TYPE_int: {
2167 140 : const int *vals = (const int *) bi.base;
2168 140 : int *avgs = Tloc(bn, 0);
2169 2540981 : TIMEOUT_LOOP(ci.ncand, timeoffset) {
2170 2540355 : o = canditer_next(&ci) - b->hseqbase;
2171 2540354 : if (ngrp > 1)
2172 289609 : gid = gids ? gids[o] - min : o;
2173 2540354 : if (is_int_nil(vals[o])) {
2174 129045 : if (!skip_nils) {
2175 0 : avgs[gid] = int_nil;
2176 0 : rems[gid] = lng_nil;
2177 0 : cnts[gid] = lng_nil;
2178 0 : bn->tnil = true;
2179 0 : rn->tnil = true;
2180 0 : cn->tnil = true;
2181 : }
2182 2411309 : } else if (!is_lng_nil(cnts[gid])) {
2183 2568961 : AVERAGE_ITER(int, vals[o], avgs[gid], rems[gid], cnts[gid]);
2184 : }
2185 : }
2186 145663 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset) {
2187 145241 : if (cnts[i] == 0) {
2188 469 : avgs[i] = int_nil;
2189 469 : bn->tnil = true;
2190 : } else
2191 : #ifdef TRUNCATE_NUMBERS
2192 : if (!is_int_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0) {
2193 : avgs[i]++;
2194 : rems[i] -= cnts[i];
2195 : }
2196 : #else
2197 144772 : if (!is_int_nil(avgs[i]) && rems[i] > 0) {
2198 57982 : if (avgs[i] < 0) {
2199 43758 : if (2*rems[i] > cnts[i]) {
2200 19501 : avgs[i]++;
2201 19501 : rems[i] -= cnts[i];
2202 : }
2203 : } else {
2204 14224 : if (2*rems[i] >= cnts[i]) {
2205 12362 : avgs[i]++;
2206 12362 : rems[i] -= cnts[i];
2207 : }
2208 : }
2209 : }
2210 : #endif
2211 : }
2212 : break;
2213 : }
2214 92 : case TYPE_lng: {
2215 92 : const lng *vals = (const lng *) bi.base;
2216 92 : lng *avgs = Tloc(bn, 0);
2217 186759 : TIMEOUT_LOOP(ci.ncand, timeoffset) {
2218 186575 : o = canditer_next(&ci) - b->hseqbase;
2219 186576 : if (ngrp > 1)
2220 154227 : gid = gids ? gids[o] - min : o;
2221 186576 : if (is_lng_nil(vals[o])) {
2222 210 : if (!skip_nils) {
2223 0 : avgs[gid] = lng_nil;
2224 0 : rems[gid] = lng_nil;
2225 0 : cnts[gid] = lng_nil;
2226 0 : bn->tnil = true;
2227 0 : rn->tnil = true;
2228 0 : cn->tnil = true;
2229 : }
2230 186366 : } else if (!is_lng_nil(cnts[gid])) {
2231 187106 : AVERAGE_ITER(lng, vals[o], avgs[gid], rems[gid], cnts[gid]);
2232 : }
2233 : }
2234 69453 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset) {
2235 69177 : if (cnts[i] == 0) {
2236 152 : avgs[i] = lng_nil;
2237 152 : bn->tnil = true;
2238 : } else
2239 : #ifdef TRUNCATE_NUMBERS
2240 : if (!is_lng_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0) {
2241 : avgs[i]++;
2242 : rems[i] -= cnts[i];
2243 : }
2244 : #else
2245 69025 : if (!is_lng_nil(avgs[i]) && rems[i] > 0) {
2246 922 : if (avgs[i] < 0) {
2247 129 : if (2*rems[i] > cnts[i]) {
2248 0 : avgs[i]++;
2249 0 : rems[i] -= cnts[i];
2250 : }
2251 : } else {
2252 793 : if (2*rems[i] >= cnts[i]) {
2253 694 : avgs[i]++;
2254 694 : rems[i] -= cnts[i];
2255 : }
2256 : }
2257 : }
2258 : #endif
2259 : }
2260 : break;
2261 : }
2262 : #ifdef HAVE_HGE
2263 11 : case TYPE_hge: {
2264 11 : const hge *vals = (const hge *) bi.base;
2265 11 : hge *avgs = Tloc(bn, 0);
2266 11896810 : TIMEOUT_LOOP(ci.ncand, timeoffset) {
2267 11896065 : o = canditer_next(&ci) - b->hseqbase;
2268 11896065 : if (ngrp > 1)
2269 150707 : gid = gids ? gids[o] - min : o;
2270 11896065 : if (is_hge_nil(vals[o])) {
2271 248594 : if (!skip_nils) {
2272 0 : avgs[gid] = hge_nil;
2273 0 : rems[gid] = lng_nil;
2274 0 : cnts[gid] = lng_nil;
2275 0 : bn->tnil = true;
2276 0 : rn->tnil = true;
2277 0 : cn->tnil = true;
2278 : }
2279 11647471 : } else if (!is_lng_nil(cnts[gid])) {
2280 11896065 : AVERAGE_ITER(hge, vals[o], avgs[gid], rems[gid], cnts[gid]);
2281 : }
2282 : }
2283 157 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset) {
2284 124 : if (cnts[i] == 0) {
2285 0 : avgs[i] = hge_nil;
2286 0 : bn->tnil = true;
2287 : } else
2288 : #ifdef TRUNCATE_NUMBERS
2289 : if (!is_hge_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0) {
2290 : avgs[i]++;
2291 : rems[i] -= cnts[i];
2292 : }
2293 : #else
2294 124 : if (!is_hge_nil(avgs[i]) && rems[i] > 0) {
2295 119 : if (avgs[i] < 0) {
2296 0 : if (2*rems[i] > cnts[i]) {
2297 0 : avgs[i]++;
2298 0 : rems[i] -= cnts[i];
2299 : }
2300 : } else {
2301 119 : if (2*rems[i] >= cnts[i]) {
2302 52 : avgs[i]++;
2303 52 : rems[i] -= cnts[i];
2304 : }
2305 : }
2306 : }
2307 : #endif
2308 : }
2309 : break;
2310 : }
2311 : #endif
2312 : }
2313 252 : bat_iterator_end(&bi);
2314 252 : TIMEOUT_CHECK(timeoffset, GOTO_LABEL_TIMEOUT_HANDLER(bailout));
2315 252 : BATsetcount(bn, ngrp);
2316 252 : BATsetcount(rn, ngrp);
2317 252 : BATsetcount(cn, ngrp);
2318 251 : bn->tnonil = !bn->tnil;
2319 251 : rn->tnonil = !rn->tnil;
2320 251 : cn->tnonil = !cn->tnil;
2321 251 : bn->tkey = rn->tkey = cn->tkey = ngrp == 1;
2322 251 : bn->tsorted = rn->tsorted = cn->tsorted = ngrp == 1;
2323 251 : bn->trevsorted = rn->trevsorted = cn->trevsorted = ngrp == 1;
2324 251 : *avgp = bn;
2325 251 : *remp = rn;
2326 251 : *cntp = cn;
2327 251 : return GDK_SUCCEED;
2328 :
2329 0 : bailout:
2330 0 : BBPreclaim(bn);
2331 0 : BBPreclaim(rn);
2332 0 : BBPreclaim(cn);
2333 0 : return GDK_FAIL;
2334 : }
2335 :
2336 : #ifdef HAVE_HGE
2337 : #define BIGINT uhge
2338 : #else
2339 : #define BIGINT uint64_t
2340 : #endif
2341 : /* calculate the values (a*b)/c and (a*b)%c but without causing overflow
2342 : *
2343 : * this function is from https://stackoverflow.com/a/8757419, only with the
2344 : * slight addition of returning the remainder and changing the first argument
2345 : * and return value to BIGINT */
2346 : static inline BIGINT
2347 0 : multdiv(BIGINT a, uint64_t b, uint64_t c, uint64_t *rem)
2348 : {
2349 0 : static BIGINT const base = ((BIGINT)1)<<(sizeof(BIGINT)*4);
2350 : // static BIGINT const maxdiv = (base-1)*base + (base-1);
2351 0 : static BIGINT const maxdiv = ~(BIGINT)0;
2352 :
2353 : // First get the easy thing
2354 0 : BIGINT res = (a/c) * b + (a%c) * (b/c);
2355 0 : a %= c;
2356 0 : b %= c;
2357 : // Are we done?
2358 0 : if (a == 0 || b == 0) {
2359 0 : *rem = 0;
2360 0 : return res;
2361 : }
2362 : // Is it easy to compute what remain to be added?
2363 0 : if (c < base) {
2364 0 : *rem = (uint64_t) ((a*b)%c);
2365 0 : return res + (a*b/c);
2366 : }
2367 : // Now 0 < a < c, 0 < b < c, c >= 1ULL
2368 : // Normalize
2369 : BIGINT norm = maxdiv/c;
2370 : BIGINT cbig = c * norm; // orig: c *= norm; and below cbig was plain c
2371 : a *= norm;
2372 : // split into 2 digits
2373 : BIGINT ah = a / base, al = a % base;
2374 : BIGINT bh = b / base, bl = b % base;
2375 : BIGINT ch = cbig / base, cl = cbig % base;
2376 : // compute the product
2377 : BIGINT p0 = al*bl;
2378 : BIGINT p1 = p0 / base + al*bh;
2379 : p0 %= base;
2380 : BIGINT p2 = p1 / base + ah*bh;
2381 : p1 = (p1 % base) + ah * bl;
2382 : p2 += p1 / base;
2383 : p1 %= base;
2384 : // p2 holds 2 digits, p1 and p0 one
2385 :
2386 : // first digit is easy, not null only in case of overflow
2387 : //BIGINT q2 = p2 / cbig;
2388 : p2 = p2 % cbig;
2389 :
2390 : // second digit, estimate
2391 : BIGINT q1 = p2 / ch;
2392 : // and now adjust
2393 : BIGINT rhat = p2 % ch;
2394 : // the loop can be unrolled, it will be executed at most twice for
2395 : // even bases -- three times for odd one -- due to the normalisation above
2396 : while (q1 >= base || (rhat < base && q1*cl > rhat*base+p1)) {
2397 : q1--;
2398 : rhat += ch;
2399 : }
2400 : // subtract
2401 : p1 = ((p2 % base) * base + p1) - q1 * cl;
2402 : p2 = (p2 / base * base + p1 / base) - q1 * ch;
2403 : p1 = p1 % base + (p2 % base) * base;
2404 :
2405 : // now p1 hold 2 digits, p0 one and p2 is to be ignored
2406 : BIGINT q0 = p1 / ch;
2407 : rhat = p1 % ch;
2408 : while (q0 >= base || (rhat < base && q0*cl > rhat*base+p0)) {
2409 : q0--;
2410 : rhat += ch;
2411 : }
2412 : // subtract
2413 : p0 = ((p1 % base) * base + p0) - q0 * cl;
2414 : p1 = (p1 / base * base + p0 / base) - q0 * ch;
2415 : p0 = p0 % base + (p1 % base) * base;
2416 :
2417 : *rem = p0 / norm;
2418 : return res + q0 + q1 * base; // + q2 *base*base
2419 : }
2420 :
2421 : static inline void
2422 8 : combine_averages_bte(bte *avgp, lng *remp, lng *cntp,
2423 : bte avg1, lng rem1, lng cnt1,
2424 : bte avg2, lng rem2, lng cnt2)
2425 : {
2426 8 : lng cnt = cnt1 + cnt2;
2427 :
2428 8 : if (rem2 < 0) {
2429 2 : avg2--;
2430 2 : rem2 += cnt2;
2431 : }
2432 8 : *cntp = cnt;
2433 8 : lng v = avg1 * cnt1 + rem1 + avg2 * cnt2 + rem2;
2434 8 : bte a = (bte) (v / cnt);
2435 8 : v %= cnt;
2436 8 : if (v < 0) {
2437 0 : a--;
2438 0 : v += cnt;
2439 : }
2440 8 : *avgp = a;
2441 8 : *remp = v;
2442 8 : }
2443 :
2444 : static inline void
2445 1 : combine_averages_sht(sht *avgp, lng *remp, lng *cntp,
2446 : sht avg1, lng rem1, lng cnt1,
2447 : sht avg2, lng rem2, lng cnt2)
2448 : {
2449 1 : lng cnt = cnt1 + cnt2;
2450 :
2451 1 : if (rem2 < 0) {
2452 0 : avg2--;
2453 0 : rem2 += cnt2;
2454 : }
2455 1 : *cntp = cnt;
2456 1 : lng v = avg1 * cnt1 + rem1 + avg2 * cnt2 + rem2;
2457 1 : sht a = (sht) (v / cnt);
2458 1 : v %= cnt;
2459 1 : if (v < 0) {
2460 0 : a--;
2461 0 : v += cnt;
2462 : }
2463 1 : *avgp = a;
2464 1 : *remp = v;
2465 1 : }
2466 :
2467 :
2468 : static inline void
2469 81629 : combine_averages_int(int *avgp, lng *remp, lng *cntp,
2470 : int avg1, lng rem1, lng cnt1,
2471 : int avg2, lng rem2, lng cnt2)
2472 : {
2473 81629 : lng cnt = cnt1 + cnt2;
2474 :
2475 81629 : if (rem2 < 0) {
2476 23952 : avg2--;
2477 23952 : rem2 += cnt2;
2478 : }
2479 81629 : *cntp = cnt;
2480 : #ifdef HAVE_HGE
2481 81629 : hge v = avg1 * cnt1 + rem1 + avg2 * cnt2 + rem2;
2482 81629 : int a = (int) (v / cnt);
2483 81629 : v %= cnt;
2484 81629 : if (v < 0) {
2485 58824 : a--;
2486 58824 : v += cnt;
2487 : }
2488 81629 : *avgp = a;
2489 81629 : *remp = (lng) v;
2490 : #else
2491 : if (cnt1 == 0) {
2492 : avg1 = 0;
2493 : rem1 = 0;
2494 : }
2495 : lng rem = rem1 + rem2;
2496 : lng v;
2497 : uint64_t r;
2498 : if (avg1 < 0) {
2499 : avg1 = (int) multdiv((uint64_t) -avg1, cnt1, cnt, &r);
2500 : if (r > 0) {
2501 : avg1 = -avg1 - 1;
2502 : r = cnt - r;
2503 : } else {
2504 : avg1 = -avg1;
2505 : }
2506 : } else {
2507 : avg1 = (int) multdiv((uint64_t) avg1, cnt1, cnt, &r);
2508 : }
2509 : v = avg1;
2510 : rem += r;
2511 : if (avg2 < 0) {
2512 : avg2 = (int) multdiv((uint64_t) -avg2, cnt2, cnt, &r);
2513 : if (r > 0) {
2514 : avg2 = -avg2 - 1;
2515 : r = cnt - r;
2516 : } else {
2517 : avg2 = -avg2;
2518 : }
2519 : } else {
2520 : avg2 = (int) multdiv((uint64_t) avg2, cnt2, cnt, &r);
2521 : }
2522 : v += avg2;
2523 : rem += r;
2524 : while (rem >= cnt) { /* max twice */
2525 : v++;
2526 : rem -= cnt;
2527 : }
2528 : *avgp = (int) v;
2529 : *remp = rem;
2530 : #endif
2531 81629 : }
2532 :
2533 : static inline void
2534 52 : combine_averages_lng(lng *avgp, lng *remp, lng *cntp,
2535 : lng avg1, lng rem1, lng cnt1,
2536 : lng avg2, lng rem2, lng cnt2)
2537 : {
2538 52 : lng cnt = cnt1 + cnt2;
2539 :
2540 52 : if (rem2 < 0) {
2541 24 : avg2--;
2542 24 : rem2 += cnt2;
2543 : }
2544 52 : *cntp = cnt;
2545 : #ifdef HAVE_HGE
2546 52 : hge v = avg1 * cnt1 + rem1 + avg2 * cnt2 + rem2;
2547 52 : lng a = (lng) (v / cnt);
2548 52 : v %= cnt;
2549 52 : if (v < 0) {
2550 0 : a--;
2551 0 : v += cnt;
2552 : }
2553 52 : *avgp = a;
2554 52 : *remp = (lng) v;
2555 : #else
2556 : if (cnt1 == 0) {
2557 : avg1 = 0;
2558 : rem1 = 0;
2559 : }
2560 : lng rem = rem1 + rem2;
2561 : lng v;
2562 : uint64_t r;
2563 : if (avg1 < 0) {
2564 : avg1 = (lng) multdiv((uint64_t) -avg1, cnt1, cnt, &r);
2565 : if (r > 0) {
2566 : avg1 = -avg1 - 1;
2567 : r = cnt - r;
2568 : } else {
2569 : avg1 = -avg1;
2570 : }
2571 : } else {
2572 : avg1 = (lng) multdiv((uint64_t) avg1, cnt1, cnt, &r);
2573 : }
2574 : v = avg1;
2575 : rem += r;
2576 : if (avg2 < 0) {
2577 : avg2 = (lng) multdiv((uint64_t) -avg2, cnt2, cnt, &r);
2578 : if (r > 0) {
2579 : avg2 = -avg2 - 1;
2580 : r = cnt - r;
2581 : } else {
2582 : avg2 = -avg2;
2583 : }
2584 : } else {
2585 : avg2 = (lng) multdiv((uint64_t) avg2, cnt2, cnt, &r);
2586 : }
2587 : v += avg2;
2588 : rem += r;
2589 : while (rem >= cnt) { /* max twice */
2590 : v++;
2591 : rem -= cnt;
2592 : }
2593 : *avgp = v;
2594 : *remp = rem;
2595 : #endif
2596 52 : }
2597 :
2598 : #ifdef HAVE_HGE
2599 : static inline void
2600 0 : combine_averages_hge(hge *avgp, lng *remp, lng *cntp,
2601 : hge avg1, lng rem1, lng cnt1,
2602 : hge avg2, lng rem2, lng cnt2)
2603 : {
2604 0 : if (cnt1 == 0) {
2605 0 : avg1 = 0;
2606 0 : rem1 = 0;
2607 : }
2608 0 : if (rem2 < 0) {
2609 0 : avg2--;
2610 0 : rem2 += cnt2;
2611 : }
2612 0 : lng cnt = cnt1 + cnt2;
2613 0 : lng rem = rem1 + rem2;
2614 0 : hge v;
2615 0 : uint64_t r;
2616 :
2617 0 : *cntp = cnt;
2618 0 : if (avg1 < 0) {
2619 0 : avg1 = (hge) multdiv((uhge) -avg1, cnt1, cnt, &r);
2620 0 : if (r > 0) {
2621 0 : avg1 = -avg1 - 1;
2622 0 : r = cnt - r;
2623 : } else {
2624 0 : avg1 = -avg1;
2625 : }
2626 : } else {
2627 0 : avg1 = (hge) multdiv((uhge) avg1, cnt1, cnt, &r);
2628 : }
2629 0 : v = avg1;
2630 0 : rem += r;
2631 0 : if (avg2 < 0) {
2632 0 : avg2 = (hge) multdiv((uhge) -avg2, cnt2, cnt, &r);
2633 0 : if (r > 0) {
2634 0 : avg2 = -avg2 - 1;
2635 0 : r = cnt - r;
2636 : } else {
2637 0 : avg2 = -avg2;
2638 : }
2639 : } else {
2640 0 : avg2 = (hge) multdiv((uhge) avg2, cnt2, cnt, &r);
2641 : }
2642 0 : v += avg2;
2643 0 : rem += r;
2644 0 : while (rem >= cnt) { /* max twice */
2645 0 : v++;
2646 0 : rem -= cnt;
2647 : }
2648 0 : *avgp = v;
2649 0 : *remp = rem;
2650 0 : }
2651 : #endif
2652 :
2653 : BAT *
2654 39 : BATgroupavg3combine(BAT *avg, BAT *rem, BAT *cnt, BAT *g, BAT *e, bool skip_nils)
2655 : {
2656 39 : const char *err;
2657 39 : oid min, max;
2658 39 : BUN ngrp;
2659 39 : struct canditer ci;
2660 39 : BUN i;
2661 39 : BAT *bn, *rn, *cn;
2662 :
2663 39 : lng timeoffset = 0;
2664 39 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
2665 40 : if (qry_ctx != NULL) {
2666 40 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
2667 : }
2668 :
2669 40 : if ((err = BATgroupaggrinit(avg, g, e, NULL, &min, &max, &ngrp, &ci)) != NULL) {
2670 0 : GDKerror("%s\n", err);
2671 0 : return NULL;
2672 : }
2673 40 : assert(ci.tpe == cand_dense);
2674 40 : if (ci.ncand != BATcount(rem) || ci.ncand != BATcount(cnt)) {
2675 0 : GDKerror("input bats not aligned");
2676 0 : return NULL;
2677 : }
2678 40 : if (ci.ncand == 0 || ngrp == 0) {
2679 1 : return BATconstant(ngrp == 0 ? 0 : min, avg->ttype,
2680 1 : ATOMnilptr(avg->ttype), ngrp, TRANSIENT);
2681 : }
2682 39 : ValRecord zero;
2683 39 : (void) VALinit(&zero, TYPE_bte, &(bte){0});
2684 39 : bn = BATconstant(min, avg->ttype, VALconvert(avg->ttype, &zero),
2685 : ngrp, TRANSIENT);
2686 : /* rn and cn are temporary storage of intermediates */
2687 39 : rn = BATconstant(min, TYPE_lng, &(lng){0}, ngrp, TRANSIENT);
2688 39 : cn = BATconstant(min, TYPE_lng, &(lng){0}, ngrp, TRANSIENT);
2689 39 : if (bn == NULL || rn == NULL || cn == NULL) {
2690 0 : BBPreclaim(bn);
2691 0 : BBPreclaim(rn);
2692 0 : BBPreclaim(cn);
2693 0 : return NULL;
2694 : }
2695 39 : lng *rems = Tloc(rn, 0);
2696 39 : lng *cnts = Tloc(cn, 0);
2697 39 : const lng *orems = Tloc(rem, 0);
2698 39 : const lng *ocnts = Tloc(cnt, 0);
2699 39 : const oid *gids = g && !BATtdense(g) ? Tloc(g, 0) : NULL;
2700 39 : oid gid = ngrp == 1 && gids ? gids[0] - min : 0;
2701 :
2702 39 : BATiter bi = bat_iterator(avg);
2703 :
2704 76 : switch (ATOMbasetype(avg->ttype)) {
2705 2 : case TYPE_bte: {
2706 2 : const bte *vals = (const bte *) bi.base;
2707 2 : bte *avgs = Tloc(bn, 0);
2708 14 : TIMEOUT_LOOP_IDX(i, ci.ncand, timeoffset) {
2709 8 : if (ngrp > 1)
2710 0 : gid = gids ? gids[i] - min : i;
2711 8 : if (is_bte_nil(vals[i])) {
2712 0 : if (!skip_nils) {
2713 0 : avgs[gid] = bte_nil;
2714 0 : bn->tnil = true;
2715 : }
2716 8 : } else if (!is_bte_nil(avgs[gid])) {
2717 8 : combine_averages_bte(&avgs[gid], &rems[gid],
2718 : &cnts[gid], avgs[gid],
2719 8 : rems[gid], cnts[gid],
2720 8 : vals[i], orems[i],
2721 8 : ocnts[i]);
2722 : }
2723 : }
2724 8 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset) {
2725 2 : if (cnts[i] == 0) {
2726 0 : avgs[i] = bte_nil;
2727 0 : bn->tnil = true;
2728 : } else
2729 : #ifdef TRUNCATE_NUMBERS
2730 : if (!is_bte_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0)
2731 : avgs[i]++;
2732 : #else
2733 2 : if (!is_bte_nil(avgs[i]) && rems[i] > 0) {
2734 2 : if (avgs[i] < 0) {
2735 0 : if (2*rems[i] > cnts[i])
2736 0 : avgs[i]++;
2737 : } else {
2738 2 : if (2*rems[i] >= cnts[i])
2739 0 : avgs[i]++;
2740 : }
2741 : }
2742 : #endif
2743 : }
2744 : break;
2745 : }
2746 1 : case TYPE_sht: {
2747 1 : const sht *vals = (const sht *) bi.base;
2748 1 : sht *avgs = Tloc(bn, 0);
2749 4 : TIMEOUT_LOOP_IDX(i, ci.ncand, timeoffset) {
2750 1 : if (ngrp > 1)
2751 0 : gid = gids ? gids[i] - min : i;
2752 1 : if (is_sht_nil(vals[i])) {
2753 0 : if (!skip_nils) {
2754 0 : avgs[gid] = sht_nil;
2755 0 : bn->tnil = true;
2756 : }
2757 1 : } else if (!is_sht_nil(avgs[gid])) {
2758 1 : combine_averages_sht(&avgs[gid], &rems[gid],
2759 : &cnts[gid], avgs[gid],
2760 1 : rems[gid], cnts[gid],
2761 1 : vals[i], orems[i],
2762 1 : ocnts[i]);
2763 : }
2764 : }
2765 4 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset) {
2766 1 : if (cnts[i] == 0) {
2767 0 : avgs[i] = sht_nil;
2768 0 : bn->tnil = true;
2769 : } else
2770 : #ifdef TRUNCATE_NUMBERS
2771 : if (!is_sht_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0)
2772 : avgs[i]++;
2773 : #else
2774 1 : if (!is_sht_nil(avgs[i]) && rems[i] > 0) {
2775 1 : if (avgs[i] < 0) {
2776 0 : if (2*rems[i] > cnts[i])
2777 0 : avgs[i]++;
2778 : } else {
2779 1 : if (2*rems[i] >= cnts[i])
2780 0 : avgs[i]++;
2781 : }
2782 : }
2783 : #endif
2784 : }
2785 : break;
2786 : }
2787 31 : case TYPE_int: {
2788 31 : const int *vals = (const int *) bi.base;
2789 31 : int *avgs = Tloc(bn, 0);
2790 81866 : TIMEOUT_LOOP_IDX(i, ci.ncand, timeoffset) {
2791 81767 : if (ngrp > 1)
2792 81692 : gid = gids ? gids[i] - min : i;
2793 81767 : if (is_int_nil(vals[i])) {
2794 138 : if (!skip_nils) {
2795 0 : avgs[gid] = int_nil;
2796 0 : bn->tnil = true;
2797 : }
2798 81629 : } else if (!is_int_nil(avgs[gid])) {
2799 81629 : combine_averages_int(&avgs[gid], &rems[gid],
2800 : &cnts[gid], avgs[gid],
2801 81629 : rems[gid], cnts[gid],
2802 81629 : vals[i], orems[i],
2803 81629 : ocnts[i]);
2804 : }
2805 : }
2806 34852 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset) {
2807 34756 : if (cnts[i] == 0) {
2808 42 : avgs[i] = int_nil;
2809 42 : bn->tnil = true;
2810 : } else
2811 : #ifdef TRUNCATE_NUMBERS
2812 : if (!is_int_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0)
2813 : avgs[i]++;
2814 : #else
2815 34714 : if (!is_int_nil(avgs[i]) && rems[i] > 0) {
2816 19886 : if (avgs[i] < 0) {
2817 16717 : if (2*rems[i] > cnts[i])
2818 8094 : avgs[i]++;
2819 : } else {
2820 3169 : if (2*rems[i] >= cnts[i])
2821 2904 : avgs[i]++;
2822 : }
2823 : }
2824 : #endif
2825 : }
2826 : break;
2827 : }
2828 4 : case TYPE_lng: {
2829 4 : const lng *vals = (const lng *) bi.base;
2830 4 : lng *avgs = Tloc(bn, 0);
2831 64 : TIMEOUT_LOOP_IDX(i, ci.ncand, timeoffset) {
2832 52 : if (ngrp > 1)
2833 48 : gid = gids ? gids[i] - min : i;
2834 52 : if (is_lng_nil(vals[i])) {
2835 0 : if (!skip_nils) {
2836 0 : avgs[gid] = lng_nil;
2837 0 : bn->tnil = true;
2838 : }
2839 52 : } else if (!is_lng_nil(avgs[gid])) {
2840 52 : combine_averages_lng(&avgs[gid], &rems[gid],
2841 : &cnts[gid], avgs[gid],
2842 52 : rems[gid], cnts[gid],
2843 52 : vals[i], orems[i],
2844 52 : ocnts[i]);
2845 : }
2846 : }
2847 25 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset) {
2848 13 : if (cnts[i] == 0) {
2849 0 : avgs[i] = lng_nil;
2850 0 : bn->tnil = true;
2851 : } else
2852 : #ifdef TRUNCATE_NUMBERS
2853 : if (!is_lng_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0)
2854 : avgs[i]++;
2855 : #else
2856 13 : if (!is_lng_nil(avgs[i]) && rems[i] > 0) {
2857 12 : if (avgs[i] < 0) {
2858 0 : if (2*rems[i] > cnts[i])
2859 0 : avgs[i]++;
2860 : } else {
2861 12 : if (2*rems[i] >= cnts[i])
2862 10 : avgs[i]++;
2863 : }
2864 : }
2865 : #endif
2866 : }
2867 : break;
2868 : }
2869 : #ifdef HAVE_HGE
2870 0 : case TYPE_hge: {
2871 0 : const hge *vals = (const hge *) bi.base;
2872 0 : hge *avgs = Tloc(bn, 0);
2873 0 : TIMEOUT_LOOP_IDX(i, ci.ncand, timeoffset) {
2874 0 : if (ngrp > 1)
2875 0 : gid = gids ? gids[i] - min : i;
2876 0 : if (is_hge_nil(vals[i])) {
2877 0 : if (!skip_nils) {
2878 0 : avgs[gid] = hge_nil;
2879 0 : bn->tnil = true;
2880 : }
2881 0 : } else if (!is_hge_nil(avgs[gid])) {
2882 0 : combine_averages_hge(&avgs[gid], &rems[gid],
2883 : &cnts[gid], avgs[gid],
2884 0 : rems[gid], cnts[gid],
2885 0 : vals[i], orems[i],
2886 0 : ocnts[i]);
2887 : }
2888 : }
2889 0 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset) {
2890 0 : if (cnts[i] == 0) {
2891 0 : avgs[i] = hge_nil;
2892 0 : bn->tnil = true;
2893 : } else
2894 : #ifdef TRUNCATE_NUMBERS
2895 : if (!is_hge_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0)
2896 : avgs[i]++;
2897 : #else
2898 0 : if (!is_hge_nil(avgs[i]) && rems[i] > 0) {
2899 0 : if (avgs[i] < 0) {
2900 0 : if (2*rems[i] > cnts[i])
2901 0 : avgs[i]++;
2902 : } else {
2903 0 : if (2*rems[i] >= cnts[i])
2904 0 : avgs[i]++;
2905 : }
2906 : }
2907 : #endif
2908 : }
2909 : break;
2910 : }
2911 : #endif
2912 : }
2913 39 : bat_iterator_end(&bi);
2914 39 : BBPreclaim(rn);
2915 39 : BBPreclaim(cn);
2916 39 : TIMEOUT_CHECK(timeoffset, GOTO_LABEL_TIMEOUT_HANDLER(bailout));
2917 39 : BATsetcount(bn, ngrp);
2918 39 : bn->tnonil = !bn->tnil;
2919 39 : bn->tkey = ngrp == 1;
2920 39 : bn->tsorted = ngrp == 1;
2921 39 : bn->trevsorted = ngrp == 1;
2922 39 : return bn;
2923 :
2924 0 : bailout:
2925 0 : BBPreclaim(bn);
2926 0 : return NULL;
2927 : }
2928 :
2929 : #define AVERAGE_TYPE_LNG_HGE(TYPE,lng_hge) \
2930 : do { \
2931 : TYPE x, a; \
2932 : \
2933 : /* first try to calculate the sum of all values into a */ \
2934 : /* lng/hge */ \
2935 : TIMEOUT_LOOP(ci.ncand, timeoffset) { \
2936 : i = canditer_next(&ci) - b->hseqbase; \
2937 : x = ((const TYPE *) src)[i]; \
2938 : if (is_##TYPE##_nil(x)) \
2939 : continue; \
2940 : ADD_WITH_CHECK(x, sum, \
2941 : lng_hge, sum, \
2942 : GDK_##lng_hge##_max, \
2943 : goto overflow##TYPE); \
2944 : /* don't count value until after overflow check */ \
2945 : n++; \
2946 : } \
2947 : /* the sum fit, so now we can calculate the average */ \
2948 : *avg = n > 0 ? (dbl) sum / n : dbl_nil; \
2949 : if (0) { \
2950 : overflow##TYPE: \
2951 : /* we get here if sum(x[0],...,x[i]) doesn't */ \
2952 : /* fit in a lng/hge but sum(x[0],...,x[i-1]) did */ \
2953 : /* the variable sum contains that sum */ \
2954 : /* the rest of the calculation is done */ \
2955 : /* according to the loop invariant described */ \
2956 : /* in the below loop */ \
2957 : /* note that n necessarily is > 0 (else no */ \
2958 : /* overflow possible) */ \
2959 : assert(n > 0); \
2960 : if (sum >= 0) { \
2961 : a = (TYPE) (sum / n); /* this fits */ \
2962 : r = (lng) (sum % n); \
2963 : } else { \
2964 : sum = -sum; \
2965 : a = - (TYPE) (sum / n); /* this fits */ \
2966 : r = (lng) (sum % n); \
2967 : if (r) { \
2968 : a--; \
2969 : r = n - r; \
2970 : } \
2971 : } \
2972 : TIMEOUT_LOOP(ci.ncand, timeoffset) { \
2973 : /* loop invariant: */ \
2974 : /* a + r/n == average(x[0],...,x[n]); */ \
2975 : /* 0 <= r < n */ \
2976 : i = canditer_next(&ci) - b->hseqbase; \
2977 : x = ((const TYPE *) src)[i]; \
2978 : if (is_##TYPE##_nil(x)) \
2979 : continue; \
2980 : AVERAGE_ITER(TYPE, x, a, r, n); \
2981 : } \
2982 : *avg = a + (dbl) r / n; \
2983 : } \
2984 : TIMEOUT_CHECK(timeoffset, \
2985 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
2986 : } while (0)
2987 :
2988 : #ifdef HAVE_HGE
2989 : #define AVERAGE_TYPE(TYPE) AVERAGE_TYPE_LNG_HGE(TYPE,hge)
2990 : #else
2991 : #define AVERAGE_TYPE(TYPE) AVERAGE_TYPE_LNG_HGE(TYPE,lng)
2992 : #endif
2993 :
2994 : #define AVERAGE_FLOATTYPE(TYPE) \
2995 : do { \
2996 : double a = 0; \
2997 : TYPE x; \
2998 : TIMEOUT_LOOP(ci.ncand, timeoffset) { \
2999 : i = canditer_next(&ci) - b->hseqbase; \
3000 : x = ((const TYPE *) src)[i]; \
3001 : if (is_##TYPE##_nil(x)) \
3002 : continue; \
3003 : AVERAGE_ITER_FLOAT(TYPE, x, a, n); \
3004 : } \
3005 : TIMEOUT_CHECK(timeoffset, \
3006 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
3007 : *avg = n > 0 ? a : dbl_nil; \
3008 : } while (0)
3009 :
3010 : gdk_return
3011 4042 : BATcalcavg(BAT *b, BAT *s, dbl *avg, BUN *vals, int scale)
3012 : {
3013 4042 : lng n = 0, r = 0;
3014 4042 : BUN i = 0;
3015 : #ifdef HAVE_HGE
3016 4042 : hge sum = 0;
3017 : #else
3018 : lng sum = 0;
3019 : #endif
3020 4042 : struct canditer ci;
3021 4042 : const void *restrict src;
3022 :
3023 4042 : lng timeoffset = 0;
3024 4042 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
3025 4042 : if (qry_ctx != NULL) {
3026 4042 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
3027 : }
3028 :
3029 4042 : canditer_init(&ci, b, s);
3030 :
3031 4041 : BATiter bi = bat_iterator(b);
3032 4042 : src = bi.base;
3033 :
3034 4042 : switch (b->ttype) {
3035 8 : case TYPE_bte:
3036 24 : AVERAGE_TYPE(bte);
3037 : break;
3038 0 : case TYPE_sht:
3039 0 : AVERAGE_TYPE(sht);
3040 : break;
3041 3982 : case TYPE_int:
3042 4742280 : AVERAGE_TYPE(int);
3043 : break;
3044 9 : case TYPE_lng:
3045 10000637 : AVERAGE_TYPE(lng);
3046 : break;
3047 : #ifdef HAVE_HGE
3048 0 : case TYPE_hge:
3049 0 : AVERAGE_TYPE(hge);
3050 : break;
3051 : #endif
3052 1 : case TYPE_flt:
3053 100006105 : AVERAGE_FLOATTYPE(flt);
3054 1 : break;
3055 42 : case TYPE_dbl:
3056 155 : AVERAGE_FLOATTYPE(dbl);
3057 42 : break;
3058 0 : default:
3059 0 : GDKerror("average of type %s unsupported.\n",
3060 : ATOMname(bi.type));
3061 0 : goto bailout;
3062 : }
3063 4042 : bat_iterator_end(&bi);
3064 4042 : if (scale != 0 && !is_dbl_nil(*avg))
3065 0 : *avg /= pow(10.0, (double) scale);
3066 4042 : if (vals)
3067 4042 : *vals = (BUN) n;
3068 : return GDK_SUCCEED;
3069 0 : bailout:
3070 0 : bat_iterator_end(&bi);
3071 0 : return GDK_FAIL;
3072 : }
3073 :
3074 : /* ---------------------------------------------------------------------- */
3075 : /* count */
3076 :
3077 : #define AGGR_COUNT(TYPE) \
3078 : do { \
3079 : const TYPE *restrict vals = (const TYPE *) bi.base; \
3080 : TIMEOUT_LOOP(ci.ncand, timeoffset) { \
3081 : i = canditer_next(&ci) - b->hseqbase; \
3082 : if (gids == NULL || \
3083 : (gids[i] >= min && gids[i] <= max)) { \
3084 : if (gids) \
3085 : gid = gids[i] - min; \
3086 : else \
3087 : gid = (oid) i; \
3088 : if (!is_##TYPE##_nil(vals[i])) { \
3089 : cnts[gid]++; \
3090 : } \
3091 : } \
3092 : } \
3093 : } while (0)
3094 :
3095 : /* calculate group counts with optional candidates list */
3096 : BAT *
3097 11472 : BATgroupcount(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
3098 : {
3099 11472 : const oid *restrict gids;
3100 11472 : oid gid;
3101 11472 : oid min, max;
3102 11472 : BUN i, ngrp;
3103 11472 : lng *restrict cnts;
3104 11472 : BAT *bn = NULL;
3105 11472 : int t;
3106 11472 : const void *nil;
3107 11472 : int (*atomcmp)(const void *, const void *);
3108 11472 : struct canditer ci;
3109 11472 : const char *err;
3110 11472 : lng t0 = 0;
3111 11472 : BATiter bi = {0};
3112 :
3113 11472 : lng timeoffset = 0;
3114 11472 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
3115 11472 : if (qry_ctx != NULL) {
3116 11472 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
3117 : }
3118 :
3119 11472 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
3120 :
3121 11472 : assert(tp == TYPE_lng);
3122 11472 : (void) tp; /* compatibility (with other BATgroup* */
3123 :
3124 11472 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
3125 0 : GDKerror("%s\n", err);
3126 0 : return NULL;
3127 : }
3128 11472 : if (g == NULL) {
3129 0 : GDKerror("b and g must be aligned\n");
3130 0 : return NULL;
3131 : }
3132 :
3133 11472 : if (ci.ncand == 0 || ngrp == 0) {
3134 : /* trivial: no counts, so return bat aligned with g
3135 : * with zero in the tail */
3136 5050 : lng zero = 0;
3137 5050 : return BATconstant(ngrp == 0 ? 0 : min, TYPE_lng, &zero, ngrp, TRANSIENT);
3138 : }
3139 :
3140 6422 : bn = COLnew(min, TYPE_lng, ngrp, TRANSIENT);
3141 6422 : if (bn == NULL)
3142 : return NULL;
3143 6422 : cnts = (lng *) Tloc(bn, 0);
3144 6422 : memset(cnts, 0, ngrp * sizeof(lng));
3145 :
3146 6422 : if (BATtdense(g))
3147 : gids = NULL;
3148 : else
3149 4874 : gids = (const oid *) Tloc(g, 0);
3150 :
3151 6422 : if (!skip_nils || b->tnonil) {
3152 : /* if nils are nothing special, or if there are no
3153 : * nils, we don't need to look at the values at all */
3154 6186 : if (gids) {
3155 17342638 : TIMEOUT_LOOP(ci.ncand, timeoffset) {
3156 17331920 : i = canditer_next(&ci) - b->hseqbase;
3157 17331921 : if (gids[i] >= min && gids[i] <= max)
3158 17331820 : cnts[gids[i] - min]++;
3159 : }
3160 : } else {
3161 92425 : TIMEOUT_LOOP(ci.ncand, timeoffset) {
3162 89709 : i = canditer_next(&ci) - b->hseqbase;
3163 89710 : cnts[i] = 1;
3164 : }
3165 : }
3166 : } else {
3167 236 : t = b->ttype;
3168 236 : nil = ATOMnilptr(t);
3169 236 : atomcmp = ATOMcompare(t);
3170 236 : t = ATOMbasetype(t);
3171 :
3172 236 : bi = bat_iterator(b);
3173 :
3174 236 : switch (t) {
3175 6 : case TYPE_bte:
3176 20 : AGGR_COUNT(bte);
3177 : break;
3178 0 : case TYPE_sht:
3179 0 : AGGR_COUNT(sht);
3180 : break;
3181 210 : case TYPE_int:
3182 16180 : AGGR_COUNT(int);
3183 : break;
3184 4 : case TYPE_lng:
3185 12 : AGGR_COUNT(lng);
3186 : break;
3187 : #ifdef HAVE_HGE
3188 0 : case TYPE_hge:
3189 0 : AGGR_COUNT(hge);
3190 : break;
3191 : #endif
3192 0 : case TYPE_flt:
3193 0 : AGGR_COUNT(flt);
3194 : break;
3195 5 : case TYPE_dbl:
3196 155 : AGGR_COUNT(dbl);
3197 : break;
3198 11 : default:
3199 199 : TIMEOUT_LOOP(ci.ncand, timeoffset) {
3200 177 : i = canditer_next(&ci) - b->hseqbase;
3201 178 : if (gids == NULL ||
3202 176 : (gids[i] >= min && gids[i] <= max)) {
3203 176 : if (gids)
3204 176 : gid = gids[i] - min;
3205 : else
3206 : gid = (oid) i;
3207 178 : if ((*atomcmp)(BUNtail(bi, i), nil) != 0) {
3208 92 : cnts[gid]++;
3209 : }
3210 : }
3211 : }
3212 : break;
3213 : }
3214 236 : bat_iterator_end(&bi);
3215 : }
3216 6422 : TIMEOUT_CHECK(timeoffset, GOTO_LABEL_TIMEOUT_HANDLER(bailout));
3217 6422 : BATsetcount(bn, ngrp);
3218 6422 : bn->tkey = BATcount(bn) <= 1;
3219 6422 : bn->tsorted = BATcount(bn) <= 1;
3220 6422 : bn->trevsorted = BATcount(bn) <= 1;
3221 6422 : bn->tnil = false;
3222 6422 : bn->tnonil = true;
3223 6422 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
3224 : "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
3225 : "; start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
3226 : ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
3227 : ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
3228 : ci.seq, ci.ncand, GDKusec() - t0);
3229 : return bn;
3230 :
3231 0 : bailout:
3232 0 : BBPreclaim(bn);
3233 0 : return NULL;
3234 : }
3235 :
3236 : /* ---------------------------------------------------------------------- */
3237 : /* min and max */
3238 :
3239 : #define AGGR_CMP(TYPE, OP) \
3240 : do { \
3241 : const TYPE *restrict vals = (const TYPE *) bi->base; \
3242 : if (ngrp == ci->ncand) { \
3243 : /* single element groups */ \
3244 : TIMEOUT_LOOP(ci->ncand, timeoffset) { \
3245 : i = canditer_next(ci) - hseq; \
3246 : if (!skip_nils || \
3247 : !is_##TYPE##_nil(vals[i])) { \
3248 : oids[gid] = i + hseq; \
3249 : nils--; \
3250 : } \
3251 : gid++; \
3252 : } \
3253 : } else { \
3254 : TIMEOUT_LOOP(ci->ncand, timeoffset) { \
3255 : i = canditer_next(ci) - hseq; \
3256 : if (gids == NULL || \
3257 : (gids[i] >= min && gids[i] <= max)) { \
3258 : if (gids) \
3259 : gid = gids[i] - min; \
3260 : if (!skip_nils || !is_##TYPE##_nil(vals[i])) { \
3261 : if (is_oid_nil(oids[gid])) { \
3262 : oids[gid] = i + hseq; \
3263 : nils--; \
3264 : } else if (!is_##TYPE##_nil(vals[oids[gid] - hseq]) && \
3265 : (is_##TYPE##_nil(vals[i]) || \
3266 : OP(vals[i], vals[oids[gid] - hseq]))) \
3267 : oids[gid] = i + hseq; \
3268 : } \
3269 : } \
3270 : } \
3271 : } \
3272 : } while (0)
3273 :
3274 : /* calculate group minimums with optional candidates list
3275 : *
3276 : * note that this functions returns *positions* of where the minimum
3277 : * values occur */
3278 : static BUN
3279 544 : do_groupmin(oid *restrict oids, BATiter *bi, const oid *restrict gids, BUN ngrp,
3280 : oid min, oid max, struct canditer *restrict ci,
3281 : bool skip_nils, bool gdense)
3282 : {
3283 544 : oid gid = 0;
3284 544 : BUN i, nils;
3285 544 : int t;
3286 544 : const void *nil;
3287 544 : int (*atomcmp)(const void *, const void *);
3288 :
3289 544 : lng timeoffset = 0;
3290 544 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
3291 544 : if (qry_ctx != NULL) {
3292 544 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
3293 : }
3294 :
3295 1088 : nils = ngrp;
3296 61701 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset)
3297 60069 : oids[i] = oid_nil;
3298 544 : if (ci->ncand == 0)
3299 : return nils;
3300 :
3301 544 : t = bi->b->ttype;
3302 544 : nil = ATOMnilptr(t);
3303 544 : atomcmp = ATOMcompare(t);
3304 544 : t = ATOMbasetype(t);
3305 544 : oid hseq = bi->b->hseqbase;
3306 :
3307 544 : switch (t) {
3308 16 : case TYPE_bte:
3309 88 : AGGR_CMP(bte, LT);
3310 : break;
3311 1 : case TYPE_sht:
3312 1760 : AGGR_CMP(sht, LT);
3313 : break;
3314 383 : case TYPE_int:
3315 78907 : AGGR_CMP(int, LT);
3316 : break;
3317 67 : case TYPE_lng:
3318 102185 : AGGR_CMP(lng, LT);
3319 : break;
3320 : #ifdef HAVE_HGE
3321 8 : case TYPE_hge:
3322 30 : AGGR_CMP(hge, LT);
3323 : break;
3324 : #endif
3325 0 : case TYPE_flt:
3326 0 : AGGR_CMP(flt, LT);
3327 : break;
3328 27 : case TYPE_dbl:
3329 318 : AGGR_CMP(dbl, LT);
3330 : break;
3331 0 : case TYPE_void:
3332 0 : if (!skip_nils || !is_oid_nil(bi->tseq)) {
3333 0 : if (!gdense && gids == NULL) {
3334 0 : oids[0] = canditer_next(ci);
3335 0 : nils--;
3336 0 : } else if (gdense) {
3337 : /* single element groups */
3338 0 : TIMEOUT_LOOP(ci->ncand, timeoffset) {
3339 0 : i = canditer_next(ci);
3340 0 : oids[gid++] = i;
3341 0 : nils--;
3342 : }
3343 : } else {
3344 0 : TIMEOUT_LOOP(ci->ncand, timeoffset) {
3345 0 : i = canditer_next(ci);
3346 0 : gid = gids[i - hseq] - min;
3347 0 : if (is_oid_nil(oids[gid])) {
3348 0 : oids[gid] = i;
3349 0 : nils--;
3350 : }
3351 : }
3352 : }
3353 : }
3354 : break;
3355 42 : default:
3356 42 : assert(bi->b->ttype != TYPE_oid);
3357 :
3358 42 : if (gdense) {
3359 : /* single element groups */
3360 3 : TIMEOUT_LOOP(ci->ncand, timeoffset) {
3361 1 : i = canditer_next(ci) - hseq;
3362 2 : if (!skip_nils ||
3363 1 : (*atomcmp)(BUNtail(*bi, i), nil) != 0) {
3364 0 : oids[gid] = i + hseq;
3365 0 : nils--;
3366 : }
3367 1 : gid++;
3368 : }
3369 : } else {
3370 22262 : TIMEOUT_LOOP(ci->ncand, timeoffset) {
3371 22180 : i = canditer_next(ci) - hseq;
3372 22180 : if (gids == NULL ||
3373 51 : (gids[i] >= min && gids[i] <= max)) {
3374 22180 : const void *v = BUNtail(*bi, i);
3375 22180 : if (gids)
3376 51 : gid = gids[i] - min;
3377 44360 : if (!skip_nils ||
3378 22180 : (*atomcmp)(v, nil) != 0) {
3379 19795 : if (is_oid_nil(oids[gid])) {
3380 60 : oids[gid] = i + hseq;
3381 60 : nils--;
3382 19735 : } else if (t != TYPE_void) {
3383 19735 : const void *g = BUNtail(*bi, (BUN) (oids[gid] - hseq));
3384 39470 : if ((*atomcmp)(g, nil) != 0 &&
3385 39470 : ((*atomcmp)(v, nil) == 0 ||
3386 19735 : LT((*atomcmp)(v, g), 0)))
3387 81 : oids[gid] = i + hseq;
3388 : }
3389 : }
3390 : }
3391 : }
3392 : }
3393 : break;
3394 : }
3395 544 : TIMEOUT_CHECK(timeoffset, TIMEOUT_HANDLER(BUN_NONE));
3396 :
3397 : return nils;
3398 : }
3399 :
3400 : /* calculate group maximums with optional candidates list
3401 : *
3402 : * note that this functions returns *positions* of where the maximum
3403 : * values occur */
3404 : static BUN
3405 535 : do_groupmax(oid *restrict oids, BATiter *bi, const oid *restrict gids, BUN ngrp,
3406 : oid min, oid max, struct canditer *restrict ci,
3407 : bool skip_nils, bool gdense)
3408 : {
3409 535 : oid gid = 0;
3410 535 : BUN i, nils;
3411 535 : int t;
3412 535 : const void *nil;
3413 535 : int (*atomcmp)(const void *, const void *);
3414 :
3415 535 : lng timeoffset = 0;
3416 535 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
3417 535 : if (qry_ctx != NULL) {
3418 535 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
3419 : }
3420 :
3421 1070 : nils = ngrp;
3422 703574 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset)
3423 701895 : oids[i] = oid_nil;
3424 535 : if (ci->ncand == 0)
3425 : return nils;
3426 :
3427 535 : t = bi->b->ttype;
3428 535 : nil = ATOMnilptr(t);
3429 535 : atomcmp = ATOMcompare(t);
3430 535 : t = ATOMbasetype(t);
3431 535 : oid hseq = bi->b->hseqbase;
3432 :
3433 535 : switch (t) {
3434 22 : case TYPE_bte:
3435 122 : AGGR_CMP(bte, GT);
3436 : break;
3437 1 : case TYPE_sht:
3438 168 : AGGR_CMP(sht, GT);
3439 : break;
3440 376 : case TYPE_int:
3441 153236 : AGGR_CMP(int, GT);
3442 : break;
3443 74 : case TYPE_lng:
3444 103898 : AGGR_CMP(lng, GT);
3445 : break;
3446 : #ifdef HAVE_HGE
3447 9 : case TYPE_hge:
3448 33860205 : AGGR_CMP(hge, GT);
3449 : break;
3450 : #endif
3451 1 : case TYPE_flt:
3452 8 : AGGR_CMP(flt, GT);
3453 : break;
3454 22 : case TYPE_dbl:
3455 253 : AGGR_CMP(dbl, GT);
3456 : break;
3457 0 : case TYPE_void:
3458 0 : if (!skip_nils || !is_oid_nil(bi->tseq)) {
3459 0 : if (!gdense && gids == NULL) {
3460 0 : oids[0] = canditer_last(ci);
3461 0 : nils--;
3462 0 : } else if (gdense) {
3463 : /* single element groups */
3464 0 : TIMEOUT_LOOP(ci->ncand, timeoffset) {
3465 0 : i = canditer_next(ci);
3466 0 : oids[gid++] = i;
3467 0 : nils--;
3468 : }
3469 : } else {
3470 0 : TIMEOUT_LOOP(ci->ncand, timeoffset) {
3471 0 : i = canditer_next(ci);
3472 0 : gid = gids[i - hseq] - min;
3473 0 : if (is_oid_nil(oids[gid]))
3474 0 : nils--;
3475 0 : oids[gid] = i;
3476 : }
3477 : }
3478 : }
3479 : break;
3480 30 : default:
3481 30 : assert(bi->b->ttype != TYPE_oid);
3482 :
3483 30 : if (gdense) {
3484 : /* single element groups */
3485 9 : TIMEOUT_LOOP(ci->ncand, timeoffset) {
3486 3 : i = canditer_next(ci) - hseq;
3487 6 : if (!skip_nils ||
3488 3 : (*atomcmp)(BUNtail(*bi, i), nil) != 0) {
3489 3 : oids[gid] = i + hseq;
3490 3 : nils--;
3491 : }
3492 3 : gid++;
3493 : }
3494 : } else {
3495 21300 : TIMEOUT_LOOP(ci->ncand, timeoffset) {
3496 21246 : i = canditer_next(ci) - hseq;
3497 21246 : if (gids == NULL ||
3498 26 : (gids[i] >= min && gids[i] <= max)) {
3499 21246 : const void *v = BUNtail(*bi, i);
3500 21246 : if (gids)
3501 26 : gid = gids[i] - min;
3502 42492 : if (!skip_nils ||
3503 21246 : (*atomcmp)(v, nil) != 0) {
3504 18863 : if (is_oid_nil(oids[gid])) {
3505 34 : oids[gid] = i + hseq;
3506 34 : nils--;
3507 : } else {
3508 18829 : const void *g = BUNtail(*bi, (BUN) (oids[gid] - hseq));
3509 18829 : if (t == TYPE_void ||
3510 37658 : ((*atomcmp)(g, nil) != 0 &&
3511 37658 : ((*atomcmp)(v, nil) == 0 ||
3512 18829 : GT((*atomcmp)(v, g), 0))))
3513 204 : oids[gid] = i + hseq;
3514 : }
3515 : }
3516 : }
3517 : }
3518 : }
3519 : break;
3520 : }
3521 535 : TIMEOUT_CHECK(timeoffset, TIMEOUT_HANDLER(BUN_NONE));
3522 :
3523 : return nils;
3524 : }
3525 :
3526 : static BAT *
3527 1027 : BATgroupminmax(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils,
3528 : BUN (*minmax)(oid *restrict, BATiter *, const oid *restrict, BUN,
3529 : oid, oid, struct canditer *restrict,
3530 : bool, bool),
3531 : const char *name)
3532 : {
3533 1027 : const oid *restrict gids;
3534 1027 : oid min, max;
3535 1027 : BUN ngrp;
3536 1027 : oid *restrict oids;
3537 1027 : BAT *bn = NULL;
3538 1027 : BUN nils;
3539 1027 : struct canditer ci;
3540 1027 : const char *err;
3541 1027 : lng t0 = 0;
3542 :
3543 1027 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
3544 :
3545 1027 : assert(tp == TYPE_oid);
3546 1027 : (void) tp; /* compatibility (with other BATgroup* */
3547 :
3548 1027 : if (!ATOMlinear(b->ttype)) {
3549 0 : GDKerror("%s: cannot determine minimum on "
3550 : "non-linear type %s\n", name, ATOMname(b->ttype));
3551 0 : return NULL;
3552 : }
3553 :
3554 1027 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
3555 0 : GDKerror("%s: %s\n", name, err);
3556 0 : return NULL;
3557 : }
3558 :
3559 1027 : if (ci.ncand == 0 || ngrp == 0) {
3560 : /* trivial: no minimums, so return bat aligned with g
3561 : * with nil in the tail */
3562 112 : return BATconstant(ngrp == 0 ? 0 : min, TYPE_oid, &oid_nil, ngrp, TRANSIENT);
3563 : }
3564 :
3565 915 : bn = COLnew(min, TYPE_oid, ngrp, TRANSIENT);
3566 915 : if (bn == NULL)
3567 : return NULL;
3568 915 : oids = (oid *) Tloc(bn, 0);
3569 :
3570 915 : if (g == NULL || BATtdense(g))
3571 : gids = NULL;
3572 : else
3573 461 : gids = (const oid *) Tloc(g, 0);
3574 :
3575 915 : BATiter bi = bat_iterator(b);
3576 915 : nils = (*minmax)(oids, &bi, gids, ngrp, min, max, &ci,
3577 915 : skip_nils, g && BATtdense(g));
3578 915 : bat_iterator_end(&bi);
3579 915 : if (nils == BUN_NONE) {
3580 0 : BBPreclaim(bn);
3581 0 : return NULL;
3582 : }
3583 :
3584 915 : BATsetcount(bn, ngrp);
3585 :
3586 915 : bn->tkey = BATcount(bn) <= 1;
3587 915 : bn->tsorted = BATcount(bn) <= 1;
3588 915 : bn->trevsorted = BATcount(bn) <= 1;
3589 915 : bn->tnil = nils != 0;
3590 915 : bn->tnonil = nils == 0;
3591 915 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
3592 : "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
3593 : "; start " OIDFMT ", count " BUNFMT " (%s -- " LLFMT " usec)\n",
3594 : ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
3595 : ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
3596 : ci.seq, ci.ncand, name, GDKusec() - t0);
3597 : return bn;
3598 : }
3599 :
3600 : BAT *
3601 494 : BATgroupmin(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
3602 : {
3603 494 : return BATgroupminmax(b, g, e, s, tp, skip_nils,
3604 : do_groupmin, __func__);
3605 : }
3606 :
3607 : /* return pointer to smallest non-nil value in b, or pointer to nil if
3608 : * there is no such value (no values at all, or only nil) */
3609 : void *
3610 1234 : BATmin_skipnil(BAT *b, void *aggr, bit skipnil)
3611 : {
3612 1234 : const void *res = NULL;
3613 1234 : size_t s;
3614 1234 : lng t0 = 0;
3615 :
3616 1234 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
3617 :
3618 1234 : if (!ATOMlinear(b->ttype)) {
3619 : /* there is no such thing as a smallest value if you
3620 : * can't compare values */
3621 0 : GDKerror("non-linear type");
3622 0 : return NULL;
3623 : }
3624 1234 : BATiter bi = bat_iterator(b);
3625 1234 : if (bi.count == 0) {
3626 287 : res = ATOMnilptr(bi.type);
3627 947 : } else if (bi.minpos != BUN_NONE) {
3628 446 : res = BUNtail(bi, bi.minpos);
3629 : } else {
3630 501 : oid pos;
3631 657 : BAT *pb = BATdescriptor(VIEWtparent(b));
3632 501 : Heap *oidxh = NULL;
3633 501 : bool usepoidx = false;
3634 :
3635 501 : if (BATordered(b)) {
3636 415 : if (skipnil && !bi.nonil) {
3637 355 : pos = binsearch(NULL, 0, bi.type, bi.base,
3638 355 : bi.vh ? bi.vh->base : NULL,
3639 355 : bi.width, 0, bi.count,
3640 355 : ATOMnilptr(bi.type), 1, 1);
3641 355 : if (pos == bi.count)
3642 48 : pos = oid_nil;
3643 : else
3644 307 : pos += b->hseqbase;
3645 : } else {
3646 60 : pos = b->hseqbase;
3647 : }
3648 86 : } else if (BATordered_rev(b)) {
3649 4 : if (skipnil && !bi.nonil) {
3650 1 : pos = binsearch(NULL, 0, bi.type, bi.base,
3651 1 : bi.vh ? bi.vh->base : NULL,
3652 1 : bi.width, 0, bi.count,
3653 1 : ATOMnilptr(bi.type), -1, 1);
3654 1 : if (pos == 0)
3655 0 : pos = oid_nil;
3656 : else
3657 1 : pos = pos + b->hseqbase - 1;
3658 : } else {
3659 3 : pos = bi.count + b->hseqbase - 1;
3660 : }
3661 : } else {
3662 82 : if (BATcheckorderidx(b)) {
3663 0 : MT_lock_set(&b->batIdxLock);
3664 0 : oidxh = b->torderidx;
3665 0 : if (oidxh != NULL)
3666 0 : HEAPincref(oidxh);
3667 0 : MT_lock_unset(&b->batIdxLock);
3668 : }
3669 82 : if (oidxh == NULL &&
3670 121 : pb != NULL &&
3671 39 : BATcheckorderidx(pb)) {
3672 : /* no lock on b needed since it's a view */
3673 0 : MT_lock_set(&pb->batIdxLock);
3674 0 : MT_lock_set(&pb->theaplock);
3675 0 : if (pb->tbaseoff == bi.baseoff &&
3676 0 : BATcount(pb) == bi.count &&
3677 0 : pb->hseqbase == b->hseqbase &&
3678 0 : (oidxh = pb->torderidx) != NULL) {
3679 0 : HEAPincref(oidxh);
3680 0 : usepoidx = true;
3681 : }
3682 0 : MT_lock_unset(&pb->batIdxLock);
3683 0 : MT_lock_unset(&pb->theaplock);
3684 : }
3685 82 : if (oidxh != NULL) {
3686 0 : const oid *ords = (const oid *) oidxh->base + ORDERIDXOFF;
3687 0 : BUN r;
3688 0 : if (skipnil && !bi.nonil) {
3689 0 : MT_thread_setalgorithm(usepoidx ? "binsearch on parent oidx" : "binsearch on oidx");
3690 0 : r = binsearch(ords, 0, bi.type, bi.base,
3691 0 : bi.vh ? bi.vh->base : NULL,
3692 0 : bi.width, 0, bi.count,
3693 0 : ATOMnilptr(bi.type), 1, 1);
3694 0 : if (r == 0) {
3695 : /* there are no nils, record that */
3696 0 : MT_lock_set(&b->theaplock);
3697 0 : if (b->batCount == bi.count) {
3698 0 : b->tnonil = true;
3699 : }
3700 0 : MT_lock_unset(&b->theaplock);
3701 : }
3702 : } else {
3703 : r = 0;
3704 : }
3705 0 : if (r == bi.count) {
3706 : /* no non-nil values */
3707 0 : pos = oid_nil;
3708 : } else {
3709 0 : MT_thread_setalgorithm(usepoidx ? "using parent oidx" : "using oidx");
3710 0 : pos = ords[r];
3711 : }
3712 0 : HEAPdecref(oidxh, false);
3713 : } else {
3714 82 : Imprints *imprints = NULL;
3715 149 : if ((pb == NULL || bi.count == BATcount(pb)) &&
3716 67 : BATcheckimprints(b)) {
3717 0 : if (pb != NULL) {
3718 0 : BAT *pb = BBP_cache(VIEWtparent(b));
3719 0 : MT_lock_set(&pb->batIdxLock);
3720 0 : imprints = pb->timprints;
3721 0 : if (imprints != NULL)
3722 0 : IMPSincref(imprints);
3723 : else
3724 : imprints = NULL;
3725 0 : MT_lock_unset(&pb->batIdxLock);
3726 : } else {
3727 0 : MT_lock_set(&b->batIdxLock);
3728 0 : imprints = b->timprints;
3729 0 : if (imprints != NULL)
3730 0 : IMPSincref(imprints);
3731 : else
3732 : imprints = NULL;
3733 0 : MT_lock_unset(&b->batIdxLock);
3734 : }
3735 : }
3736 0 : if (imprints) {
3737 0 : int i;
3738 :
3739 0 : MT_thread_setalgorithm(VIEWtparent(b) ? "using parent imprints" : "using imprints");
3740 0 : pos = oid_nil;
3741 : /* find first non-empty bin */
3742 0 : for (i = 0; i < imprints->bits; i++) {
3743 0 : if (imprints->stats[i + 128]) {
3744 0 : pos = imprints->stats[i] + b->hseqbase;
3745 0 : break;
3746 : }
3747 : }
3748 0 : IMPSdecref(imprints, false);
3749 : } else {
3750 82 : struct canditer ci;
3751 82 : canditer_init(&ci, b, NULL);
3752 82 : (void) do_groupmin(&pos, &bi, NULL, 1, 0, 0, &ci,
3753 : skipnil, false);
3754 : }
3755 : }
3756 : }
3757 501 : if (is_oid_nil(pos)) {
3758 48 : res = ATOMnilptr(bi.type);
3759 : } else {
3760 453 : bi.minpos = pos - b->hseqbase;
3761 453 : res = BUNtail(bi, bi.minpos);
3762 453 : MT_lock_set(&b->theaplock);
3763 452 : if (bi.count == BATcount(b) && bi.h == b->theap)
3764 453 : b->tminpos = bi.minpos;
3765 452 : MT_lock_unset(&b->theaplock);
3766 453 : if (pb) {
3767 331 : MT_lock_set(&pb->theaplock);
3768 331 : if (bi.count == BATcount(pb) &&
3769 96 : bi.h == pb->theap)
3770 96 : pb->tminpos = bi.minpos;
3771 331 : MT_lock_unset(&pb->theaplock);
3772 : }
3773 : }
3774 846 : BBPreclaim(pb);
3775 : }
3776 1234 : if (aggr == NULL) {
3777 526 : s = ATOMlen(bi.type, res);
3778 526 : aggr = GDKmalloc(s);
3779 : } else {
3780 1416 : s = ATOMsize(ATOMtype(bi.type));
3781 : }
3782 1234 : if (aggr != NULL) /* else: malloc error */
3783 1234 : memcpy(aggr, res, s);
3784 1234 : bat_iterator_end(&bi);
3785 1234 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",skipnil=%d; (" LLFMT " usec)\n",
3786 : ALGOBATPAR(b), skipnil, GDKusec() - t0);
3787 : return aggr;
3788 : }
3789 :
3790 : void *
3791 474 : BATmin(BAT *b, void *aggr)
3792 : {
3793 474 : return BATmin_skipnil(b, aggr, 1);
3794 : }
3795 :
3796 : BAT *
3797 533 : BATgroupmax(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
3798 : {
3799 533 : return BATgroupminmax(b, g, e, s, tp, skip_nils,
3800 : do_groupmax, __func__);
3801 : }
3802 :
3803 : void *
3804 1378 : BATmax_skipnil(BAT *b, void *aggr, bit skipnil)
3805 : {
3806 1378 : const void *res = NULL;
3807 1378 : size_t s;
3808 1378 : lng t0 = 0;
3809 :
3810 1378 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
3811 :
3812 1378 : if (!ATOMlinear(b->ttype)) {
3813 : /* there is no such thing as a largest value if you
3814 : * can't compare values */
3815 0 : GDKerror("non-linear type");
3816 0 : return NULL;
3817 : }
3818 1378 : BATiter bi = bat_iterator(b);
3819 1378 : if (bi.count == 0) {
3820 202 : res = ATOMnilptr(bi.type);
3821 1176 : } else if (bi.maxpos != BUN_NONE) {
3822 734 : res = BUNtail(bi, bi.maxpos);
3823 : } else {
3824 442 : oid pos;
3825 575 : BAT *pb = BATdescriptor(VIEWtparent(b));
3826 442 : Heap *oidxh = NULL;
3827 442 : bool usepoidx = false;
3828 :
3829 442 : if (BATordered(b)) {
3830 351 : pos = bi.count - 1 + b->hseqbase;
3831 636 : if (skipnil && !bi.nonil &&
3832 285 : ATOMcmp(bi.type, BUNtail(bi, bi.count - 1),
3833 : ATOMnilptr(bi.type)) == 0)
3834 38 : pos = oid_nil; /* no non-nil values */
3835 91 : } else if (BATordered_rev(b)) {
3836 9 : pos = b->hseqbase;
3837 15 : if (skipnil && !bi.nonil &&
3838 6 : ATOMcmp(bi.type, BUNtail(bi, 0),
3839 : ATOMnilptr(bi.type)) == 0)
3840 0 : pos = oid_nil; /* no non-nil values */
3841 : } else {
3842 82 : if (BATcheckorderidx(b)) {
3843 0 : MT_lock_set(&b->batIdxLock);
3844 0 : oidxh = b->torderidx;
3845 0 : if (oidxh != NULL)
3846 0 : HEAPincref(oidxh);
3847 0 : MT_lock_unset(&b->batIdxLock);
3848 : }
3849 82 : if (oidxh == NULL &&
3850 110 : pb != NULL &&
3851 28 : BATcheckorderidx(pb)) {
3852 : /* no lock on b needed since it's a view */
3853 0 : MT_lock_set(&pb->batIdxLock);
3854 0 : MT_lock_set(&pb->theaplock);
3855 0 : if (pb->tbaseoff == bi.baseoff &&
3856 0 : BATcount(pb) == bi.count &&
3857 0 : pb->hseqbase == b->hseqbase &&
3858 0 : (oidxh = pb->torderidx) != NULL) {
3859 0 : HEAPincref(oidxh);
3860 0 : usepoidx = true;
3861 : }
3862 0 : MT_lock_unset(&pb->batIdxLock);
3863 0 : MT_lock_unset(&pb->theaplock);
3864 : }
3865 82 : if (oidxh != NULL) {
3866 0 : const oid *ords = (const oid *) oidxh->base + ORDERIDXOFF;
3867 :
3868 0 : MT_thread_setalgorithm(usepoidx ? "using parent oidx" : "using oids");
3869 0 : pos = ords[bi.count - 1];
3870 : /* nils are first, ie !skipnil, check for nils */
3871 0 : if (!skipnil) {
3872 0 : BUN z = ords[0];
3873 :
3874 0 : res = BUNtail(bi, z - b->hseqbase);
3875 :
3876 0 : if (ATOMcmp(bi.type, res, ATOMnilptr(bi.type)) == 0)
3877 0 : pos = z;
3878 : }
3879 0 : HEAPdecref(oidxh, false);
3880 : } else {
3881 82 : Imprints *imprints = NULL;
3882 147 : if ((pb == NULL || BATcount(b) == BATcount(pb)) &&
3883 65 : BATcheckimprints(b)) {
3884 0 : if (pb != NULL) {
3885 0 : BAT *pb = BBP_cache(VIEWtparent(b));
3886 0 : MT_lock_set(&pb->batIdxLock);
3887 0 : imprints = pb->timprints;
3888 0 : if (imprints != NULL)
3889 0 : IMPSincref(imprints);
3890 : else
3891 : imprints = NULL;
3892 0 : MT_lock_unset(&pb->batIdxLock);
3893 : } else {
3894 0 : MT_lock_set(&b->batIdxLock);
3895 0 : imprints = b->timprints;
3896 0 : if (imprints != NULL)
3897 0 : IMPSincref(imprints);
3898 : else
3899 : imprints = NULL;
3900 0 : MT_lock_unset(&b->batIdxLock);
3901 : }
3902 : }
3903 0 : if (imprints) {
3904 0 : int i;
3905 :
3906 0 : MT_thread_setalgorithm(VIEWtparent(b) ? "using parent imprints" : "using imprints");
3907 0 : pos = oid_nil;
3908 : /* find last non-empty bin */
3909 0 : for (i = imprints->bits - 1; i >= 0; i--) {
3910 0 : if (imprints->stats[i + 128]) {
3911 0 : pos = imprints->stats[i + 64] + b->hseqbase;
3912 0 : break;
3913 : }
3914 : }
3915 0 : IMPSdecref(imprints, false);
3916 : } else {
3917 82 : struct canditer ci;
3918 82 : canditer_init(&ci, b, NULL);
3919 82 : (void) do_groupmax(&pos, &bi, NULL, 1, 0, 0, &ci,
3920 : skipnil, false);
3921 : }
3922 : }
3923 : }
3924 442 : if (is_oid_nil(pos)) {
3925 38 : res = ATOMnilptr(bi.type);
3926 : } else {
3927 404 : bi.maxpos = pos - b->hseqbase;
3928 404 : res = BUNtail(bi, bi.maxpos);
3929 404 : MT_lock_set(&b->theaplock);
3930 404 : if (bi.count == BATcount(b) && bi.h == b->theap)
3931 404 : b->tmaxpos = bi.maxpos;
3932 404 : MT_lock_unset(&b->theaplock);
3933 404 : if (pb) {
3934 301 : MT_lock_set(&pb->theaplock);
3935 301 : if (bi.count == BATcount(pb) &&
3936 112 : bi.h == pb->theap)
3937 112 : pb->tmaxpos = bi.maxpos;
3938 301 : MT_lock_unset(&pb->theaplock);
3939 : }
3940 : }
3941 751 : BBPreclaim(pb);
3942 : }
3943 1378 : if (aggr == NULL) {
3944 515 : s = ATOMlen(bi.type, res);
3945 515 : aggr = GDKmalloc(s);
3946 : } else {
3947 1724 : s = ATOMsize(ATOMtype(bi.type));
3948 : }
3949 1378 : if (aggr != NULL) /* else: malloc error */
3950 1378 : memcpy(aggr, res, s);
3951 1378 : bat_iterator_end(&bi);
3952 1378 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",skipnil=%d; (" LLFMT " usec)\n",
3953 : ALGOBATPAR(b), skipnil, GDKusec() - t0);
3954 : return aggr;
3955 : }
3956 :
3957 : void *
3958 788 : BATmax(BAT *b, void *aggr)
3959 : {
3960 788 : return BATmax_skipnil(b, aggr, 1);
3961 : }
3962 :
3963 :
3964 : /* ---------------------------------------------------------------------- */
3965 : /* quantiles/median */
3966 :
3967 : #if SIZEOF_OID == SIZEOF_INT
3968 : #define binsearch_oid(indir, offset, vals, lo, hi, v, ordering, last) binsearch_int(indir, offset, (const int *) vals, lo, hi, (int) (v), ordering, last)
3969 : #endif
3970 : #if SIZEOF_OID == SIZEOF_LNG
3971 : #define binsearch_oid(indir, offset, vals, lo, hi, v, ordering, last) binsearch_lng(indir, offset, (const lng *) vals, lo, hi, (lng) (v), ordering, last)
3972 : #endif
3973 :
3974 : #define DO_QUANTILE_AVG(TPE) \
3975 : do { \
3976 : BUN idxlo, idxhi; \
3977 : if (ords) { \
3978 : idxlo = ords[r + (BUN) lo] - b->hseqbase; \
3979 : idxhi = ords[r + (BUN) hi] - b->hseqbase; \
3980 : } else { \
3981 : idxlo = r + (BUN) lo; \
3982 : idxhi = r + (BUN) hi; \
3983 : } \
3984 : TPE low = *(TPE*) BUNtloc(bi, idxhi); \
3985 : TPE high = *(TPE*) BUNtloc(bi, idxlo); \
3986 : if (is_##TPE##_nil(low) || is_##TPE##_nil(high)) { \
3987 : val = dbl_nil; \
3988 : nils++; \
3989 : } else { \
3990 : val = (f - lo) * low + (lo + 1 - f) * high; \
3991 : } \
3992 : } while (0)
3993 :
3994 : static BAT *
3995 66 : doBATgroupquantile(BAT *b, BAT *g, BAT *e, BAT *s, int tp, double quantile,
3996 : bool skip_nils, bool average)
3997 : {
3998 66 : BAT *origb = b;
3999 66 : BAT *origg = g;
4000 66 : oid min, max;
4001 66 : BUN ngrp;
4002 66 : BUN nils = 0;
4003 66 : BAT *bn = NULL;
4004 66 : struct canditer ci;
4005 66 : BAT *t1, *t2;
4006 66 : BATiter bi = {0};
4007 66 : const void *v;
4008 66 : const void *nil = ATOMnilptr(tp);
4009 66 : const void *dnil = nil;
4010 66 : dbl val; /* only used for average */
4011 66 : int (*atomcmp)(const void *, const void *) = ATOMcompare(tp);
4012 66 : const char *err;
4013 66 : lng t0 = 0;
4014 :
4015 66 : size_t counter = 0;
4016 66 : lng timeoffset = 0;
4017 66 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
4018 66 : if (qry_ctx != NULL) {
4019 66 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
4020 : }
4021 :
4022 66 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4023 :
4024 66 : if (average) {
4025 12 : switch (ATOMbasetype(b->ttype)) {
4026 : case TYPE_bte:
4027 : case TYPE_sht:
4028 : case TYPE_int:
4029 : case TYPE_lng:
4030 : #ifdef HAVE_HGE
4031 : case TYPE_hge:
4032 : #endif
4033 : case TYPE_flt:
4034 : case TYPE_dbl:
4035 : break;
4036 0 : default:
4037 0 : GDKerror("incompatible type\n");
4038 0 : return NULL;
4039 : }
4040 : dnil = &dbl_nil;
4041 : }
4042 66 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
4043 0 : GDKerror("%s\n", err);
4044 0 : return NULL;
4045 : }
4046 66 : assert(tp == b->ttype);
4047 66 : if (!ATOMlinear(tp)) {
4048 0 : GDKerror("cannot determine quantile on "
4049 : "non-linear type %s\n", ATOMname(tp));
4050 0 : return NULL;
4051 : }
4052 66 : if (quantile < 0 || quantile > 1) {
4053 2 : GDKerror("cannot determine quantile for "
4054 : "p=%f (p has to be in [0,1])\n", quantile);
4055 2 : return NULL;
4056 : }
4057 :
4058 64 : if (BATcount(b) == 0 || ngrp == 0 || is_dbl_nil(quantile)) {
4059 : /* trivial: no values, thus also no quantiles,
4060 : * so return bat aligned with e with nil in the tail
4061 : * The same happens for a NULL quantile */
4062 24 : return BATconstant(ngrp == 0 ? 0 : min, average ? TYPE_dbl : tp, dnil, ngrp, TRANSIENT);
4063 : }
4064 :
4065 52 : if (s) {
4066 : /* there is a candidate list, replace b (and g, if
4067 : * given) with just the values we're interested in */
4068 0 : b = BATproject(s, b);
4069 0 : if (b == NULL)
4070 : return NULL;
4071 0 : if (g) {
4072 0 : g = BATproject(s, g);
4073 0 : if (g == NULL)
4074 0 : goto bailout;
4075 : }
4076 : }
4077 :
4078 : /* we want to sort b so that we can figure out the quantile, but
4079 : * if g is given, sort g and subsort b so that we can get the
4080 : * quantile for each group */
4081 52 : if (g) {
4082 19 : const oid *restrict grps;
4083 19 : oid prev;
4084 19 : BUN p, q, r;
4085 :
4086 19 : if (BATtdense(g)) {
4087 : /* singleton groups, so calculating quantile is
4088 : * easy */
4089 1 : if (average)
4090 0 : bn = BATconvert(b, NULL, TYPE_dbl, 0, 0, 0);
4091 : else
4092 1 : bn = COLcopy(b, tp, false, TRANSIENT);
4093 1 : BAThseqbase(bn, g->tseqbase); /* deals with NULL */
4094 1 : if (b != origb)
4095 0 : BBPunfix(b->batCacheid);
4096 1 : if (g != origg)
4097 0 : BBPunfix(g->batCacheid);
4098 1 : return bn;
4099 : }
4100 18 : if (BATsort(&t1, &t2, NULL, g, NULL, NULL, false, false, false) != GDK_SUCCEED)
4101 0 : goto bailout;
4102 18 : if (g != origg)
4103 0 : BBPunfix(g->batCacheid);
4104 18 : g = t1;
4105 :
4106 18 : if (BATsort(&t1, NULL, NULL, b, t2, g, false, false, false) != GDK_SUCCEED) {
4107 0 : BBPunfix(t2->batCacheid);
4108 0 : goto bailout;
4109 : }
4110 18 : if (b != origb)
4111 0 : BBPunfix(b->batCacheid);
4112 18 : b = t1;
4113 18 : BBPunfix(t2->batCacheid);
4114 :
4115 18 : if (average)
4116 2 : bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
4117 : else
4118 16 : bn = COLnew(min, tp, ngrp, TRANSIENT);
4119 18 : if (bn == NULL)
4120 0 : goto bailout;
4121 :
4122 18 : bi = bat_iterator(b);
4123 :
4124 18 : grps = (const oid *) Tloc(g, 0);
4125 : /* for each group (r and p are the beginning and end
4126 : * of the current group, respectively) */
4127 68 : for (r = 0, q = BATcount(g); r < q; r = p) {
4128 50 : GDK_CHECK_TIMEOUT(timeoffset, counter,
4129 : GOTO_LABEL_TIMEOUT_HANDLER(bunins_failed));
4130 50 : BUN qindex;
4131 50 : prev = grps[r];
4132 : /* search for end of current group (grps is
4133 : * sorted so we can use binary search) */
4134 50 : p = binsearch_oid(NULL, 0, grps, r, q - 1, prev, 1, 1);
4135 50 : if (skip_nils && !bi.nonil) {
4136 : /* within group, locate start of non-nils */
4137 23 : r = binsearch(NULL, 0, tp, bi.base,
4138 23 : bi.vh ? bi.vh->base : NULL,
4139 23 : bi.width, r, p, nil,
4140 : 1, 1);
4141 : }
4142 50 : if (r == p) {
4143 : /* no non-nils found */
4144 8 : v = dnil;
4145 8 : nils++;
4146 42 : } else if (average) {
4147 4 : double f = (p - r - 1) * quantile;
4148 4 : double lo = floor(f);
4149 4 : double hi = ceil(f);
4150 4 : const oid *const ords = NULL;
4151 8 : switch (ATOMbasetype(tp)) {
4152 : case TYPE_bte:
4153 4 : DO_QUANTILE_AVG(bte);
4154 : break;
4155 : case TYPE_sht:
4156 0 : DO_QUANTILE_AVG(sht);
4157 : break;
4158 : case TYPE_int:
4159 0 : DO_QUANTILE_AVG(int);
4160 : break;
4161 : case TYPE_lng:
4162 0 : DO_QUANTILE_AVG(lng);
4163 : break;
4164 : #ifdef HAVE_HGE
4165 : case TYPE_hge:
4166 0 : DO_QUANTILE_AVG(hge);
4167 : break;
4168 : #endif
4169 : case TYPE_flt:
4170 0 : DO_QUANTILE_AVG(flt);
4171 : break;
4172 : case TYPE_dbl:
4173 0 : DO_QUANTILE_AVG(dbl);
4174 : break;
4175 : }
4176 : v = &val;
4177 : } else {
4178 : /* round *down* to nearest integer */
4179 38 : double f = (p - r - 1) * quantile;
4180 38 : qindex = r + p - (BUN) (p + 0.5 - f);
4181 : /* be a little paranoid about the index */
4182 38 : assert(qindex >= r && qindex < p);
4183 38 : v = BUNtail(bi, qindex);
4184 38 : if (!skip_nils && !bi.nonil)
4185 0 : nils += (*atomcmp)(v, dnil) == 0;
4186 : }
4187 50 : if (bunfastapp_nocheck(bn, v) != GDK_SUCCEED)
4188 0 : goto bunins_failed;
4189 : }
4190 18 : bat_iterator_end(&bi);
4191 18 : nils += ngrp - BATcount(bn);
4192 18 : while (BATcount(bn) < ngrp) {
4193 0 : if (bunfastapp_nocheck(bn, dnil) != GDK_SUCCEED)
4194 0 : goto bailout;
4195 : }
4196 18 : bn->theap->dirty = true;
4197 18 : BBPunfix(g->batCacheid);
4198 : } else {
4199 33 : BUN index, r, p = BATcount(b);
4200 33 : BAT *pb = NULL;
4201 33 : const oid *ords;
4202 33 : Heap *oidxh = NULL;
4203 :
4204 62 : bn = COLnew(0, average ? TYPE_dbl : tp, 1, TRANSIENT);
4205 33 : if (bn == NULL)
4206 0 : goto bailout;
4207 :
4208 33 : t1 = NULL;
4209 :
4210 33 : if (BATcheckorderidx(b)) {
4211 0 : MT_lock_set(&b->batIdxLock);
4212 0 : oidxh = b->torderidx;
4213 0 : if (oidxh != NULL)
4214 0 : HEAPincref(oidxh);
4215 0 : MT_lock_unset(&b->batIdxLock);
4216 : }
4217 4 : if (oidxh == NULL &&
4218 33 : VIEWtparent(b) &&
4219 4 : (pb = BATdescriptor(VIEWtparent(b))) != NULL) {
4220 4 : if (BATcheckorderidx(pb)) {
4221 : /* no lock on b needed since it's a view */
4222 0 : MT_lock_set(&pb->batIdxLock);
4223 0 : MT_lock_set(&pb->theaplock);
4224 0 : if (pb->tbaseoff == b->tbaseoff &&
4225 0 : BATcount(pb) == BATcount(b) &&
4226 0 : pb->hseqbase == b->hseqbase &&
4227 0 : (oidxh = pb->torderidx) != NULL) {
4228 0 : HEAPincref(oidxh);
4229 : }
4230 0 : MT_lock_unset(&pb->batIdxLock);
4231 0 : MT_lock_unset(&pb->theaplock);
4232 : }
4233 4 : BBPunfix(pb->batCacheid);
4234 : }
4235 33 : if (oidxh != NULL) {
4236 0 : MT_thread_setalgorithm(pb ? "using parent oidx" : "using oids");
4237 0 : ords = (const oid *) oidxh->base + ORDERIDXOFF;
4238 : } else {
4239 33 : if (BATsort(NULL, &t1, NULL, b, NULL, g, false, false, false) != GDK_SUCCEED)
4240 0 : goto bailout;
4241 33 : if (BATtdense(t1))
4242 : ords = NULL;
4243 : else
4244 11 : ords = (const oid *) Tloc(t1, 0);
4245 : }
4246 :
4247 33 : bi = bat_iterator(b);
4248 :
4249 33 : if (skip_nils && !bi.nonil)
4250 14 : r = binsearch(ords, 0, tp, bi.base,
4251 14 : bi.vh ? bi.vh->base : NULL,
4252 14 : bi.width, 0, p,
4253 : nil, 1, 1);
4254 : else
4255 : r = 0;
4256 :
4257 33 : if (r == p) {
4258 : /* no non-nil values, so quantile is nil */
4259 : v = dnil;
4260 : nils++;
4261 30 : } else if (average) {
4262 4 : double f = (p - r - 1) * quantile;
4263 4 : double lo = floor(f);
4264 4 : double hi = ceil(f);
4265 8 : switch (ATOMbasetype(tp)) {
4266 4 : case TYPE_bte:
4267 4 : DO_QUANTILE_AVG(bte);
4268 : break;
4269 0 : case TYPE_sht:
4270 0 : DO_QUANTILE_AVG(sht);
4271 : break;
4272 0 : case TYPE_int:
4273 0 : DO_QUANTILE_AVG(int);
4274 : break;
4275 0 : case TYPE_lng:
4276 0 : DO_QUANTILE_AVG(lng);
4277 : break;
4278 : #ifdef HAVE_HGE
4279 0 : case TYPE_hge:
4280 0 : DO_QUANTILE_AVG(hge);
4281 : break;
4282 : #endif
4283 0 : case TYPE_flt:
4284 0 : DO_QUANTILE_AVG(flt);
4285 : break;
4286 0 : case TYPE_dbl:
4287 0 : DO_QUANTILE_AVG(dbl);
4288 : break;
4289 : }
4290 : v = &val;
4291 : } else {
4292 26 : double f;
4293 : /* round (p-r-1)*quantile *down* to nearest
4294 : * integer (i.e., 1.49 and 1.5 are rounded to
4295 : * 1, 1.51 is rounded to 2) */
4296 26 : f = (p - r - 1) * quantile;
4297 26 : index = r + p - (BUN) (p + 0.5 - f);
4298 26 : if (ords)
4299 10 : index = ords[index] - b->hseqbase;
4300 : else
4301 16 : index = index + t1->tseqbase;
4302 26 : v = BUNtail(bi, index);
4303 26 : nils += (*atomcmp)(v, dnil) == 0;
4304 : }
4305 33 : if (oidxh != NULL)
4306 0 : HEAPdecref(oidxh, false);
4307 33 : BBPreclaim(t1);
4308 33 : gdk_return rc = BUNappend(bn, v, false);
4309 33 : bat_iterator_end(&bi);
4310 33 : if (rc != GDK_SUCCEED)
4311 0 : goto bailout;
4312 : }
4313 :
4314 51 : if (b != origb)
4315 18 : BBPunfix(b->batCacheid);
4316 :
4317 51 : bn->tkey = BATcount(bn) <= 1;
4318 51 : bn->tsorted = BATcount(bn) <= 1;
4319 51 : bn->trevsorted = BATcount(bn) <= 1;
4320 51 : bn->tnil = nils != 0;
4321 51 : bn->tnonil = nils == 0;
4322 51 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
4323 : "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT
4324 : ",quantile=%g,average=%s -> " ALGOOPTBATFMT
4325 : "; start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
4326 : ALGOBATPAR(origb), ALGOOPTBATPAR(origg), ALGOOPTBATPAR(e),
4327 : ALGOOPTBATPAR(s), quantile, average ? "true" : "false",
4328 : ALGOOPTBATPAR(bn), ci.seq, ci.ncand, GDKusec() - t0);
4329 : return bn;
4330 :
4331 0 : bunins_failed:
4332 0 : bat_iterator_end(&bi);
4333 0 : bailout:
4334 0 : if (b && b != origb)
4335 0 : BBPunfix(b->batCacheid);
4336 0 : if (g && g != origg)
4337 0 : BBPunfix(g->batCacheid);
4338 0 : BBPreclaim(bn);
4339 : return NULL;
4340 : }
4341 :
4342 : BAT *
4343 29 : BATgroupmedian(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
4344 : bool skip_nils)
4345 : {
4346 29 : return doBATgroupquantile(b, g, e, s, tp, 0.5,
4347 : skip_nils, false);
4348 : }
4349 :
4350 : BAT *
4351 31 : BATgroupquantile(BAT *b, BAT *g, BAT *e, BAT *s, int tp, double quantile,
4352 : bool skip_nils)
4353 : {
4354 31 : return doBATgroupquantile(b, g, e, s, tp, quantile,
4355 : skip_nils, false);
4356 : }
4357 :
4358 : BAT *
4359 3 : BATgroupmedian_avg(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
4360 : {
4361 3 : return doBATgroupquantile(b, g, e, s, tp, 0.5, skip_nils, true);
4362 : }
4363 :
4364 : BAT *
4365 3 : BATgroupquantile_avg(BAT *b, BAT *g, BAT *e, BAT *s, int tp, double quantile,
4366 : bool skip_nils)
4367 : {
4368 3 : return doBATgroupquantile(b, g, e, s, tp, quantile,
4369 : skip_nils, true);
4370 : }
4371 :
4372 : /* ---------------------------------------------------------------------- */
4373 : /* standard deviation (both biased and non-biased) */
4374 :
4375 : #define AGGR_STDEV_SINGLE(TYPE) \
4376 : do { \
4377 : TYPE x; \
4378 : TIMEOUT_LOOP_IDX(i, cnt, timeoffset) { \
4379 : x = ((const TYPE *) values)[i]; \
4380 : if (is_##TYPE##_nil(x)) \
4381 : continue; \
4382 : n++; \
4383 : delta = (dbl) x - mean; \
4384 : mean += delta / n; \
4385 : m2 += delta * ((dbl) x - mean); \
4386 : if (isinf(m2)) \
4387 : goto overflow; \
4388 : } \
4389 : TIMEOUT_CHECK(timeoffset, \
4390 : TIMEOUT_HANDLER(dbl_nil)); \
4391 : } while (0)
4392 :
4393 : static dbl
4394 25 : calcvariance(dbl *restrict avgp, const void *restrict values, BUN cnt, int tp, bool issample)
4395 : {
4396 25 : BUN n = 0, i;
4397 25 : dbl mean = 0;
4398 25 : dbl m2 = 0;
4399 25 : dbl delta;
4400 :
4401 25 : lng timeoffset = 0;
4402 25 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
4403 25 : if (qry_ctx != NULL) {
4404 25 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
4405 : }
4406 :
4407 25 : switch (tp) {
4408 : case TYPE_bte:
4409 18 : AGGR_STDEV_SINGLE(bte);
4410 : break;
4411 : case TYPE_sht:
4412 18 : AGGR_STDEV_SINGLE(sht);
4413 : break;
4414 : case TYPE_int:
4415 47 : AGGR_STDEV_SINGLE(int);
4416 : break;
4417 : case TYPE_lng:
4418 36 : AGGR_STDEV_SINGLE(lng);
4419 : break;
4420 : #ifdef HAVE_HGE
4421 : case TYPE_hge:
4422 0 : AGGR_STDEV_SINGLE(hge);
4423 : break;
4424 : #endif
4425 : case TYPE_flt:
4426 0 : AGGR_STDEV_SINGLE(flt);
4427 : break;
4428 : case TYPE_dbl:
4429 29 : AGGR_STDEV_SINGLE(dbl);
4430 : break;
4431 0 : default:
4432 0 : GDKerror("type (%s) not supported.\n", ATOMname(tp));
4433 0 : return dbl_nil;
4434 : }
4435 23 : if (n <= (BUN) issample) {
4436 7 : if (avgp)
4437 0 : *avgp = dbl_nil;
4438 7 : return dbl_nil;
4439 : }
4440 16 : if (avgp)
4441 0 : *avgp = mean;
4442 16 : return m2 / (n - issample);
4443 2 : overflow:
4444 2 : GDKerror("22003!overflow in calculation.\n");
4445 2 : return dbl_nil;
4446 : }
4447 :
4448 : dbl
4449 13 : BATcalcstdev_population(dbl *avgp, BAT *b)
4450 : {
4451 13 : lng t0 = 0;
4452 :
4453 13 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4454 13 : BATiter bi = bat_iterator(b);
4455 13 : dbl v = calcvariance(avgp, bi.base, bi.count, bi.type, false);
4456 13 : bat_iterator_end(&bi);
4457 13 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
4458 : ALGOBATPAR(b), GDKusec() - t0);
4459 13 : return is_dbl_nil(v) ? dbl_nil : sqrt(v);
4460 : }
4461 :
4462 : dbl
4463 7 : BATcalcstdev_sample(dbl *avgp, BAT *b)
4464 : {
4465 7 : lng t0 = 0;
4466 :
4467 7 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4468 7 : BATiter bi = bat_iterator(b);
4469 7 : dbl v = calcvariance(avgp, bi.base, bi.count, bi.type, true);
4470 7 : bat_iterator_end(&bi);
4471 7 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
4472 : ALGOBATPAR(b), GDKusec() - t0);
4473 7 : return is_dbl_nil(v) ? dbl_nil : sqrt(v);
4474 : }
4475 :
4476 : dbl
4477 3 : BATcalcvariance_population(dbl *avgp, BAT *b)
4478 : {
4479 3 : lng t0 = 0;
4480 :
4481 3 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4482 3 : BATiter bi = bat_iterator(b);
4483 3 : dbl v = calcvariance(avgp, bi.base, bi.count, bi.type, false);
4484 3 : bat_iterator_end(&bi);
4485 3 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
4486 : ALGOBATPAR(b), GDKusec() - t0);
4487 3 : return v;
4488 : }
4489 :
4490 : dbl
4491 2 : BATcalcvariance_sample(dbl *avgp, BAT *b)
4492 : {
4493 2 : lng t0 = 0;
4494 :
4495 2 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4496 2 : BATiter bi = bat_iterator(b);
4497 2 : dbl v = calcvariance(avgp, bi.base, bi.count, bi.type, true);
4498 2 : bat_iterator_end(&bi);
4499 2 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
4500 : ALGOBATPAR(b), GDKusec() - t0);
4501 2 : return v;
4502 : }
4503 :
4504 : #define AGGR_COVARIANCE_SINGLE(TYPE) \
4505 : do { \
4506 : TYPE x, y; \
4507 : TIMEOUT_LOOP_IDX(i, cnt, timeoffset) { \
4508 : x = ((const TYPE *) v1)[i]; \
4509 : y = ((const TYPE *) v2)[i]; \
4510 : if (is_##TYPE##_nil(x) || is_##TYPE##_nil(y)) \
4511 : continue; \
4512 : n++; \
4513 : delta1 = (dbl) x - mean1; \
4514 : mean1 += delta1 / n; \
4515 : delta2 = (dbl) y - mean2; \
4516 : mean2 += delta2 / n; \
4517 : m2 += delta1 * ((dbl) y - mean2); \
4518 : if (isinf(m2)) \
4519 : goto overflow; \
4520 : } \
4521 : TIMEOUT_CHECK(timeoffset, \
4522 : TIMEOUT_HANDLER(dbl_nil)); \
4523 : } while (0)
4524 :
4525 : static dbl
4526 13 : calccovariance(const void *v1, const void *v2, BUN cnt, int tp, bool issample)
4527 : {
4528 13 : BUN n = 0, i;
4529 13 : dbl mean1 = 0, mean2 = 0, m2 = 0, delta1, delta2;
4530 :
4531 13 : lng timeoffset = 0;
4532 13 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
4533 13 : if (qry_ctx != NULL) {
4534 13 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
4535 : }
4536 :
4537 :
4538 13 : switch (tp) {
4539 : case TYPE_bte:
4540 4 : AGGR_COVARIANCE_SINGLE(bte);
4541 : break;
4542 : case TYPE_sht:
4543 0 : AGGR_COVARIANCE_SINGLE(sht);
4544 : break;
4545 : case TYPE_int:
4546 111 : AGGR_COVARIANCE_SINGLE(int);
4547 : break;
4548 : case TYPE_lng:
4549 0 : AGGR_COVARIANCE_SINGLE(lng);
4550 : break;
4551 : #ifdef HAVE_HGE
4552 : case TYPE_hge:
4553 0 : AGGR_COVARIANCE_SINGLE(hge);
4554 : break;
4555 : #endif
4556 : case TYPE_flt:
4557 0 : AGGR_COVARIANCE_SINGLE(flt);
4558 : break;
4559 : case TYPE_dbl:
4560 10 : AGGR_COVARIANCE_SINGLE(dbl);
4561 : break;
4562 0 : default:
4563 0 : GDKerror("type (%s) not supported.\n", ATOMname(tp));
4564 0 : return dbl_nil;
4565 : }
4566 12 : if (n <= (BUN) issample)
4567 1 : return dbl_nil;
4568 11 : return m2 / (n - issample);
4569 1 : overflow:
4570 1 : GDKerror("22003!overflow in calculation.\n");
4571 1 : return dbl_nil;
4572 : }
4573 :
4574 : dbl
4575 8 : BATcalccovariance_population(BAT *b1, BAT *b2)
4576 : {
4577 8 : lng t0 = 0;
4578 :
4579 8 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4580 8 : BATiter b1i = bat_iterator(b1), b2i = bat_iterator(b2);
4581 8 : dbl v = calccovariance(b1i.base, b2i.base, b1i.count, b1i.type, false);
4582 8 : bat_iterator_end(&b1i);
4583 8 : bat_iterator_end(&b2i);
4584 8 : TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT " (" LLFMT " usec)\n",
4585 : ALGOBATPAR(b1), ALGOBATPAR(b2), GDKusec() - t0);
4586 8 : return v;
4587 : }
4588 :
4589 : dbl
4590 5 : BATcalccovariance_sample(BAT *b1, BAT *b2)
4591 : {
4592 5 : lng t0 = 0;
4593 :
4594 5 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4595 5 : BATiter b1i = bat_iterator(b1), b2i = bat_iterator(b2);
4596 5 : dbl v = calccovariance(b1i.base, b2i.base, b1i.count, b1i.type, true);
4597 5 : bat_iterator_end(&b1i);
4598 5 : bat_iterator_end(&b2i);
4599 5 : TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT " (" LLFMT " usec)\n",
4600 : ALGOBATPAR(b1), ALGOBATPAR(b2), GDKusec() - t0);
4601 5 : return v;
4602 : }
4603 :
4604 : #define AGGR_CORRELATION_SINGLE(TYPE) \
4605 : do { \
4606 : TYPE x, y; \
4607 : TIMEOUT_LOOP_IDX(i, cnt, timeoffset) { \
4608 : x = ((const TYPE *) v1)[i]; \
4609 : y = ((const TYPE *) v2)[i]; \
4610 : if (is_##TYPE##_nil(x) || is_##TYPE##_nil(y)) \
4611 : continue; \
4612 : n++; \
4613 : delta1 = (dbl) x - mean1; \
4614 : mean1 += delta1 / n; \
4615 : delta2 = (dbl) y - mean2; \
4616 : mean2 += delta2 / n; \
4617 : aux = (dbl) y - mean2; \
4618 : up += delta1 * aux; \
4619 : down1 += delta1 * ((dbl) x - mean1); \
4620 : down2 += delta2 * aux; \
4621 : if (isinf(up) || isinf(down1) || isinf(down2)) \
4622 : goto overflow; \
4623 : } \
4624 : TIMEOUT_CHECK(timeoffset, \
4625 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
4626 : } while (0)
4627 :
4628 : dbl
4629 19 : BATcalccorrelation(BAT *b1, BAT *b2)
4630 : {
4631 19 : BUN n = 0, i, cnt = BATcount(b1);
4632 19 : dbl mean1 = 0, mean2 = 0, up = 0, down1 = 0, down2 = 0, delta1, delta2, aux;
4633 19 : BATiter b1i = bat_iterator(b1), b2i = bat_iterator(b2);
4634 19 : const void *v1 = b1i.base, *v2 = b2i.base;
4635 19 : int tp = b1i.type;
4636 19 : lng t0 = 0;
4637 :
4638 19 : lng timeoffset = 0;
4639 19 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
4640 19 : if (qry_ctx != NULL) {
4641 19 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
4642 : }
4643 :
4644 19 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4645 :
4646 19 : switch (tp) {
4647 : case TYPE_bte:
4648 16 : AGGR_CORRELATION_SINGLE(bte);
4649 : break;
4650 : case TYPE_sht:
4651 0 : AGGR_CORRELATION_SINGLE(sht);
4652 : break;
4653 : case TYPE_int:
4654 95 : AGGR_CORRELATION_SINGLE(int);
4655 : break;
4656 : case TYPE_lng:
4657 0 : AGGR_CORRELATION_SINGLE(lng);
4658 : break;
4659 : #ifdef HAVE_HGE
4660 : case TYPE_hge:
4661 0 : AGGR_CORRELATION_SINGLE(hge);
4662 : break;
4663 : #endif
4664 : case TYPE_flt:
4665 0 : AGGR_CORRELATION_SINGLE(flt);
4666 : break;
4667 : case TYPE_dbl:
4668 12 : AGGR_CORRELATION_SINGLE(dbl);
4669 : break;
4670 0 : default:
4671 0 : GDKerror("type (%s) not supported.\n", ATOMname(tp));
4672 0 : goto bailout;
4673 : }
4674 18 : bat_iterator_end(&b1i);
4675 18 : bat_iterator_end(&b2i);
4676 18 : if (n != 0 && down1 != 0 && down2 != 0)
4677 9 : aux = (up / n) / (sqrt(down1 / n) * sqrt(down2 / n));
4678 : else
4679 9 : aux = dbl_nil;
4680 18 : TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT " (" LLFMT " usec)\n",
4681 : ALGOBATPAR(b1), ALGOBATPAR(b2), GDKusec() - t0);
4682 : return aux;
4683 1 : overflow:
4684 1 : GDKerror("22003!overflow in calculation.\n");
4685 1 : bailout:
4686 1 : bat_iterator_end(&b1i);
4687 1 : bat_iterator_end(&b2i);
4688 1 : return dbl_nil;
4689 : }
4690 :
4691 : #define AGGR_STDEV(TYPE) \
4692 : do { \
4693 : const TYPE *restrict vals = (const TYPE *) bi.base; \
4694 : TIMEOUT_LOOP(ci.ncand, timeoffset) { \
4695 : i = canditer_next(&ci) - b->hseqbase; \
4696 : if (gids == NULL || \
4697 : (gids[i] >= min && gids[i] <= max)) { \
4698 : if (gids) \
4699 : gid = gids[i] - min; \
4700 : else \
4701 : gid = (oid) i; \
4702 : if (is_##TYPE##_nil(vals[i])) { \
4703 : if (!skip_nils) \
4704 : cnts[gid] = BUN_NONE; \
4705 : } else if (cnts[gid] != BUN_NONE) { \
4706 : cnts[gid]++; \
4707 : delta[gid] = (dbl) vals[i] - mean[gid]; \
4708 : mean[gid] += delta[gid] / cnts[gid]; \
4709 : m2[gid] += delta[gid] * ((dbl) vals[i] - mean[gid]); \
4710 : } \
4711 : } \
4712 : } \
4713 : TIMEOUT_CHECK(timeoffset, \
4714 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
4715 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset) { \
4716 : if (cnts[i] == 0 || cnts[i] == BUN_NONE) { \
4717 : dbls[i] = dbl_nil; \
4718 : mean[i] = dbl_nil; \
4719 : nils++; \
4720 : } else if (cnts[i] == 1) { \
4721 : dbls[i] = issample ? dbl_nil : 0; \
4722 : nils2++; \
4723 : } else if (isinf(m2[i])) { \
4724 : goto overflow; \
4725 : } else if (variance) { \
4726 : dbls[i] = m2[i] / (cnts[i] - issample); \
4727 : } else { \
4728 : dbls[i] = sqrt(m2[i] / (cnts[i] - issample)); \
4729 : } \
4730 : } \
4731 : } while (0)
4732 :
4733 : /* Calculate group standard deviation (population (i.e. biased) or
4734 : * sample (i.e. non-biased)) with optional candidates list.
4735 : *
4736 : * Note that this helper function is prepared to return two BATs: one
4737 : * (as return value) with the standard deviation per group, and one
4738 : * (as return argument) with the average per group. This isn't
4739 : * currently used since it doesn't fit into the mold of grouped
4740 : * aggregates. */
4741 : static BAT *
4742 18 : dogroupstdev(BAT **avgb, BAT *b, BAT *g, BAT *e, BAT *s, int tp,
4743 : bool skip_nils, bool issample, bool variance, const char *func)
4744 : {
4745 18 : const oid *restrict gids;
4746 18 : oid gid;
4747 18 : oid min, max;
4748 18 : BUN i, ngrp;
4749 18 : BUN nils = 0, nils2 = 0;
4750 18 : BUN *restrict cnts = NULL;
4751 18 : dbl *restrict dbls, *restrict mean, *restrict delta, *restrict m2;
4752 18 : BAT *bn = NULL, *an = NULL;
4753 18 : struct canditer ci;
4754 18 : const char *err;
4755 18 : lng t0 = 0;
4756 18 : BATiter bi = {0};
4757 :
4758 18 : lng timeoffset = 0;
4759 18 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
4760 18 : if (qry_ctx != NULL) {
4761 18 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
4762 : }
4763 :
4764 18 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4765 :
4766 18 : assert(tp == TYPE_dbl);
4767 18 : (void) tp; /* compatibility (with other BATgroup*
4768 : * functions) argument */
4769 :
4770 18 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
4771 0 : GDKerror("%s: %s\n", func, err);
4772 0 : return NULL;
4773 : }
4774 18 : if (g == NULL) {
4775 0 : GDKerror("%s: b and g must be aligned\n", func);
4776 0 : return NULL;
4777 : }
4778 :
4779 18 : if (BATcount(b) == 0 || ngrp == 0) {
4780 : /* trivial: no products, so return bat aligned with g
4781 : * with nil in the tail */
4782 3 : bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
4783 3 : goto doreturn;
4784 : }
4785 :
4786 15 : if ((e == NULL ||
4787 15 : (BATcount(e) == ci.ncand && e->hseqbase == ci.hseq)) &&
4788 7 : (BATtdense(g) || (g->tkey && g->tnonil)) &&
4789 4 : (issample || b->tnonil)) {
4790 : /* trivial: singleton groups, so all results are equal
4791 : * to zero (population) or nil (sample) */
4792 3 : dbl v = issample ? dbl_nil : 0;
4793 4 : bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &v, ngrp, TRANSIENT);
4794 4 : goto doreturn;
4795 : }
4796 :
4797 11 : delta = GDKmalloc(ngrp * sizeof(dbl));
4798 11 : m2 = GDKmalloc(ngrp * sizeof(dbl));
4799 11 : cnts = GDKzalloc(ngrp * sizeof(BUN));
4800 11 : if (avgb) {
4801 0 : an = COLnew(0, TYPE_dbl, ngrp, TRANSIENT);
4802 0 : *avgb = an;
4803 0 : if (an == NULL) {
4804 0 : mean = NULL;
4805 0 : goto alloc_fail;
4806 : }
4807 0 : mean = (dbl *) Tloc(an, 0);
4808 : } else {
4809 11 : mean = GDKmalloc(ngrp * sizeof(dbl));
4810 : }
4811 11 : if (mean == NULL || delta == NULL || m2 == NULL || cnts == NULL)
4812 0 : goto alloc_fail;
4813 :
4814 11 : bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
4815 11 : if (bn == NULL)
4816 0 : goto alloc_fail;
4817 11 : dbls = (dbl *) Tloc(bn, 0);
4818 :
4819 1080135 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset) {
4820 1080038 : mean[i] = 0;
4821 1080038 : delta[i] = 0;
4822 1080038 : m2[i] = 0;
4823 : }
4824 :
4825 11 : bi = bat_iterator(b);
4826 11 : if (BATtdense(g))
4827 : gids = NULL;
4828 : else
4829 8 : gids = (const oid *) Tloc(g, 0);
4830 :
4831 11 : switch (b->ttype) {
4832 0 : case TYPE_bte:
4833 0 : AGGR_STDEV(bte);
4834 : break;
4835 0 : case TYPE_sht:
4836 0 : AGGR_STDEV(sht);
4837 : break;
4838 4 : case TYPE_int:
4839 5760376 : AGGR_STDEV(int);
4840 : break;
4841 0 : case TYPE_lng:
4842 0 : AGGR_STDEV(lng);
4843 : break;
4844 : #ifdef HAVE_HGE
4845 0 : case TYPE_hge:
4846 0 : AGGR_STDEV(hge);
4847 : break;
4848 : #endif
4849 0 : case TYPE_flt:
4850 0 : AGGR_STDEV(flt);
4851 : break;
4852 7 : case TYPE_dbl:
4853 104 : AGGR_STDEV(dbl);
4854 : break;
4855 0 : default:
4856 0 : GDKerror("%s: type (%s) not supported.\n",
4857 : func, ATOMname(b->ttype));
4858 0 : goto bailout;
4859 : }
4860 9 : bat_iterator_end(&bi);
4861 9 : if (an) {
4862 0 : BATsetcount(an, ngrp);
4863 0 : an->tkey = ngrp <= 1;
4864 0 : an->tsorted = ngrp <= 1;
4865 0 : an->trevsorted = ngrp <= 1;
4866 0 : an->tnil = nils != 0;
4867 0 : an->tnonil = nils == 0;
4868 : } else {
4869 9 : GDKfree(mean);
4870 : }
4871 9 : if (issample)
4872 3 : nils += nils2;
4873 9 : GDKfree(delta);
4874 9 : GDKfree(m2);
4875 9 : GDKfree(cnts);
4876 9 : BATsetcount(bn, ngrp);
4877 9 : bn->tkey = ngrp <= 1;
4878 9 : bn->tsorted = ngrp <= 1;
4879 9 : bn->trevsorted = ngrp <= 1;
4880 9 : bn->tnil = nils != 0;
4881 9 : bn->tnonil = nils == 0;
4882 16 : doreturn:
4883 16 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOBATFMT ",e=" ALGOOPTBATFMT
4884 : ",s=" ALGOOPTBATFMT
4885 : ",skip_nils=%s,issample=%s,variance=%s -> " ALGOOPTBATFMT
4886 : ",avgb=" ALGOOPTBATFMT " (%s -- " LLFMT " usec)\n",
4887 : ALGOBATPAR(b), ALGOBATPAR(g), ALGOOPTBATPAR(e),
4888 : ALGOOPTBATPAR(s),
4889 : skip_nils ? "true" : "false",
4890 : issample ? "true" : "false",
4891 : variance ? "true" : "false",
4892 : ALGOOPTBATPAR(bn), ALGOOPTBATPAR(an),
4893 : func, GDKusec() - t0);
4894 : return bn;
4895 2 : overflow:
4896 2 : GDKerror("22003!overflow in calculation.\n");
4897 2 : bailout:
4898 2 : bat_iterator_end(&bi);
4899 2 : alloc_fail:
4900 2 : if (an)
4901 0 : BBPreclaim(an);
4902 : else
4903 2 : GDKfree(mean);
4904 2 : BBPreclaim(bn);
4905 2 : GDKfree(delta);
4906 2 : GDKfree(m2);
4907 2 : GDKfree(cnts);
4908 2 : return NULL;
4909 : }
4910 :
4911 : BAT *
4912 7 : BATgroupstdev_sample(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
4913 : {
4914 7 : return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, true, false,
4915 : __func__);
4916 : }
4917 :
4918 : BAT *
4919 5 : BATgroupstdev_population(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
4920 : bool skip_nils)
4921 : {
4922 5 : return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, false, false,
4923 : __func__);
4924 : }
4925 :
4926 : BAT *
4927 1 : BATgroupvariance_sample(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
4928 : bool skip_nils)
4929 : {
4930 1 : return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, true, true,
4931 : __func__);
4932 : }
4933 :
4934 : BAT *
4935 5 : BATgroupvariance_population(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
4936 : bool skip_nils)
4937 : {
4938 5 : return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, false, true,
4939 : __func__);
4940 : }
4941 :
4942 : #define AGGR_COVARIANCE(TYPE) \
4943 : do { \
4944 : const TYPE *vals1 = (const TYPE *) b1i.base; \
4945 : const TYPE *vals2 = (const TYPE *) b2i.base; \
4946 : TIMEOUT_LOOP(ci.ncand, timeoffset) { \
4947 : i = canditer_next(&ci) - b1->hseqbase; \
4948 : if (gids == NULL || \
4949 : (gids[i] >= min && gids[i] <= max)) { \
4950 : if (gids) \
4951 : gid = gids[i] - min; \
4952 : else \
4953 : gid = (oid) i; \
4954 : if (is_##TYPE##_nil(vals1[i]) || is_##TYPE##_nil(vals2[i])) { \
4955 : if (!skip_nils) \
4956 : cnts[gid] = BUN_NONE; \
4957 : } else if (cnts[gid] != BUN_NONE) { \
4958 : cnts[gid]++; \
4959 : delta1[gid] = (dbl) vals1[i] - mean1[gid]; \
4960 : mean1[gid] += delta1[gid] / cnts[gid]; \
4961 : delta2[gid] = (dbl) vals2[i] - mean2[gid]; \
4962 : mean2[gid] += delta2[gid] / cnts[gid]; \
4963 : m2[gid] += delta1[gid] * ((dbl) vals2[i] - mean2[gid]); \
4964 : } \
4965 : } \
4966 : } \
4967 : TIMEOUT_CHECK(timeoffset, \
4968 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
4969 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset) { \
4970 : if (cnts[i] == 0 || cnts[i] == BUN_NONE) { \
4971 : dbls[i] = dbl_nil; \
4972 : nils++; \
4973 : } else if (cnts[i] == 1) { \
4974 : dbls[i] = issample ? dbl_nil : 0; \
4975 : nils2++; \
4976 : } else if (isinf(m2[i])) { \
4977 : goto overflow; \
4978 : } else { \
4979 : dbls[i] = m2[i] / (cnts[i] - issample); \
4980 : } \
4981 : } \
4982 : } while (0)
4983 :
4984 : static BAT *
4985 60 : dogroupcovariance(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp,
4986 : bool skip_nils, bool issample, const char *func)
4987 : {
4988 60 : const oid *restrict gids;
4989 60 : oid gid, min, max;
4990 60 : BUN i, ngrp, nils = 0, nils2 = 0;
4991 60 : BUN *restrict cnts = NULL;
4992 60 : dbl *restrict dbls, *restrict mean1, *restrict mean2, *restrict delta1, *restrict delta2, *restrict m2;
4993 60 : BAT *bn = NULL;
4994 60 : struct canditer ci;
4995 60 : const char *err;
4996 60 : lng t0 = 0;
4997 60 : BATiter b1i = {0}, b2i = {0};
4998 :
4999 60 : lng timeoffset = 0;
5000 60 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
5001 60 : if (qry_ctx != NULL) {
5002 60 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
5003 : }
5004 :
5005 :
5006 60 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
5007 :
5008 60 : assert(tp == TYPE_dbl && BATcount(b1) == BATcount(b2) && b1->ttype == b2->ttype && BATtdense(b1) == BATtdense(b2));
5009 60 : (void) tp;
5010 :
5011 60 : if ((err = BATgroupaggrinit(b1, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
5012 0 : GDKerror("%s: %s\n", func, err);
5013 0 : return NULL;
5014 : }
5015 60 : if (g == NULL) {
5016 0 : GDKerror("%s: b1, b2 and g must be aligned\n", func);
5017 0 : return NULL;
5018 : }
5019 :
5020 60 : if (BATcount(b1) == 0 || ngrp == 0) {
5021 0 : bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
5022 0 : goto doreturn;
5023 : }
5024 :
5025 60 : if ((e == NULL ||
5026 60 : (BATcount(e) == BATcount(b1) && (e->hseqbase == b1->hseqbase || e->hseqbase == b2->hseqbase))) &&
5027 20 : (BATtdense(g) || (g->tkey && g->tnonil)) &&
5028 10 : (issample || (b1->tnonil && b2->tnonil))) {
5029 : /* trivial: singleton groups, so all results are equal
5030 : * to zero (population) or nil (sample) */
5031 10 : dbl v = issample ? dbl_nil : 0;
5032 11 : bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &v, ngrp, TRANSIENT);
5033 11 : goto doreturn;
5034 : }
5035 :
5036 49 : delta1 = GDKmalloc(ngrp * sizeof(dbl));
5037 49 : delta2 = GDKmalloc(ngrp * sizeof(dbl));
5038 49 : m2 = GDKmalloc(ngrp * sizeof(dbl));
5039 49 : cnts = GDKzalloc(ngrp * sizeof(BUN));
5040 49 : mean1 = GDKmalloc(ngrp * sizeof(dbl));
5041 49 : mean2 = GDKmalloc(ngrp * sizeof(dbl));
5042 :
5043 49 : if (mean1 == NULL || mean2 == NULL || delta1 == NULL || delta2 == NULL || m2 == NULL || cnts == NULL)
5044 0 : goto alloc_fail;
5045 :
5046 477 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset) {
5047 330 : m2[i] = 0;
5048 330 : mean1[i] = 0;
5049 330 : mean2[i] = 0;
5050 : }
5051 :
5052 49 : bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
5053 49 : if (bn == NULL)
5054 0 : goto alloc_fail;
5055 49 : dbls = (dbl *) Tloc(bn, 0);
5056 :
5057 49 : if (!g || BATtdense(g))
5058 : gids = NULL;
5059 : else
5060 49 : gids = (const oid *) Tloc(g, 0);
5061 :
5062 49 : b1i = bat_iterator(b1);
5063 49 : b2i = bat_iterator(b2);
5064 49 : switch (b1->ttype) {
5065 19 : case TYPE_bte:
5066 392 : AGGR_COVARIANCE(bte);
5067 : break;
5068 0 : case TYPE_sht:
5069 0 : AGGR_COVARIANCE(sht);
5070 : break;
5071 30 : case TYPE_int:
5072 613 : AGGR_COVARIANCE(int);
5073 : break;
5074 0 : case TYPE_lng:
5075 0 : AGGR_COVARIANCE(lng);
5076 : break;
5077 : #ifdef HAVE_HGE
5078 0 : case TYPE_hge:
5079 0 : AGGR_COVARIANCE(hge);
5080 : break;
5081 : #endif
5082 0 : case TYPE_flt:
5083 0 : AGGR_COVARIANCE(flt);
5084 : break;
5085 0 : case TYPE_dbl:
5086 0 : AGGR_COVARIANCE(dbl);
5087 : break;
5088 0 : default:
5089 0 : GDKerror("%s: type (%s) not supported.\n", func, ATOMname(b1->ttype));
5090 0 : goto bailout;
5091 : }
5092 49 : bat_iterator_end(&b1i);
5093 49 : bat_iterator_end(&b2i);
5094 49 : GDKfree(mean1);
5095 49 : GDKfree(mean2);
5096 :
5097 49 : if (issample)
5098 20 : nils += nils2;
5099 49 : GDKfree(delta1);
5100 48 : GDKfree(delta2);
5101 49 : GDKfree(m2);
5102 49 : GDKfree(cnts);
5103 49 : BATsetcount(bn, ngrp);
5104 49 : bn->tkey = ngrp <= 1;
5105 49 : bn->tsorted = ngrp <= 1;
5106 49 : bn->trevsorted = ngrp <= 1;
5107 49 : bn->tnil = nils != 0;
5108 49 : bn->tnonil = nils == 0;
5109 60 : doreturn:
5110 60 : TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT ",g=" ALGOBATFMT
5111 : ",e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT
5112 : ",skip_nils=%s,issample=%s -> " ALGOOPTBATFMT
5113 : " (%s -- " LLFMT " usec)\n",
5114 : ALGOBATPAR(b1), ALGOBATPAR(b2), ALGOBATPAR(g),
5115 : ALGOOPTBATPAR(e), ALGOOPTBATPAR(s),
5116 : skip_nils ? "true" : "false",
5117 : issample ? "true" : "false",
5118 : ALGOOPTBATPAR(bn),
5119 : func, GDKusec() - t0);
5120 : return bn;
5121 0 : overflow:
5122 0 : GDKerror("22003!overflow in calculation.\n");
5123 0 : bailout:
5124 0 : bat_iterator_end(&b1i);
5125 0 : bat_iterator_end(&b2i);
5126 0 : alloc_fail:
5127 0 : BBPreclaim(bn);
5128 0 : GDKfree(mean1);
5129 0 : GDKfree(mean2);
5130 0 : GDKfree(delta1);
5131 0 : GDKfree(delta2);
5132 0 : GDKfree(m2);
5133 0 : GDKfree(cnts);
5134 0 : return NULL;
5135 : }
5136 :
5137 : BAT *
5138 30 : BATgroupcovariance_sample(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
5139 : {
5140 30 : return dogroupcovariance(b1, b2, g, e, s, tp, skip_nils, true,
5141 : __func__);
5142 : }
5143 :
5144 : BAT *
5145 30 : BATgroupcovariance_population(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
5146 : {
5147 30 : return dogroupcovariance(b1, b2, g, e, s, tp, skip_nils, false,
5148 : __func__);
5149 : }
5150 :
5151 : #define AGGR_CORRELATION(TYPE) \
5152 : do { \
5153 : const TYPE *vals1 = (const TYPE *) b1i.base; \
5154 : const TYPE *vals2 = (const TYPE *) b2i.base; \
5155 : TIMEOUT_LOOP(ci.ncand, timeoffset) { \
5156 : i = canditer_next(&ci) - b1->hseqbase; \
5157 : if (gids == NULL || \
5158 : (gids[i] >= min && gids[i] <= max)) { \
5159 : if (gids) \
5160 : gid = gids[i] - min; \
5161 : else \
5162 : gid = (oid) i; \
5163 : if (is_##TYPE##_nil(vals1[i]) || is_##TYPE##_nil(vals2[i])) { \
5164 : if (!skip_nils) \
5165 : cnts[gid] = BUN_NONE; \
5166 : } else if (cnts[gid] != BUN_NONE) { \
5167 : cnts[gid]++; \
5168 : delta1[gid] = (dbl) vals1[i] - mean1[gid]; \
5169 : mean1[gid] += delta1[gid] / cnts[gid]; \
5170 : delta2[gid] = (dbl) vals2[i] - mean2[gid]; \
5171 : mean2[gid] += delta2[gid] / cnts[gid]; \
5172 : aux = (dbl) vals2[i] - mean2[gid]; \
5173 : up[gid] += delta1[gid] * aux; \
5174 : down1[gid] += delta1[gid] * ((dbl) vals1[i] - mean1[gid]); \
5175 : down2[gid] += delta2[gid] * aux; \
5176 : } \
5177 : } \
5178 : } \
5179 : TIMEOUT_CHECK(timeoffset, \
5180 : GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
5181 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset) { \
5182 : if (cnts[i] <= 1 || cnts[i] == BUN_NONE || down1[i] == 0 || down2[i] == 0) { \
5183 : dbls[i] = dbl_nil; \
5184 : nils++; \
5185 : } else if (isinf(up[i]) || isinf(down1[i]) || isinf(down2[i])) { \
5186 : goto overflow; \
5187 : } else { \
5188 : dbls[i] = (up[i] / cnts[i]) / (sqrt(down1[i] / cnts[i]) * sqrt(down2[i] / cnts[i])); \
5189 : assert(!is_dbl_nil(dbls[i])); \
5190 : } \
5191 : } \
5192 : } while (0)
5193 :
5194 : BAT *
5195 49 : BATgroupcorrelation(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
5196 : {
5197 49 : const oid *restrict gids;
5198 49 : oid gid, min, max;
5199 49 : BUN i, ngrp, nils = 0;
5200 49 : BUN *restrict cnts = NULL;
5201 49 : dbl *restrict dbls, *restrict mean1, *restrict mean2, *restrict delta1, *restrict delta2, *restrict up, *restrict down1, *restrict down2, aux;
5202 49 : BAT *bn = NULL;
5203 49 : struct canditer ci;
5204 49 : const char *err;
5205 49 : lng t0 = 0;
5206 49 : BATiter b1i = {0}, b2i = {0};
5207 :
5208 49 : lng timeoffset = 0;
5209 49 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
5210 49 : if (qry_ctx != NULL) {
5211 49 : timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
5212 : }
5213 :
5214 49 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
5215 :
5216 49 : assert(tp == TYPE_dbl && BATcount(b1) == BATcount(b2) && b1->ttype == b2->ttype && BATtdense(b1) == BATtdense(b2));
5217 49 : (void) tp;
5218 :
5219 49 : if ((err = BATgroupaggrinit(b1, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
5220 0 : GDKerror("%s\n", err);
5221 0 : return NULL;
5222 : }
5223 49 : if (g == NULL) {
5224 0 : GDKerror("b1, b2 and g must be aligned\n");
5225 0 : return NULL;
5226 : }
5227 :
5228 49 : if (BATcount(b1) == 0 || ngrp == 0) {
5229 1 : bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
5230 1 : goto doreturn;
5231 : }
5232 :
5233 48 : if ((e == NULL ||
5234 48 : (BATcount(e) == BATcount(b1) && (e->hseqbase == b1->hseqbase || e->hseqbase == b2->hseqbase))) &&
5235 14 : (BATtdense(g) || (g->tkey && g->tnonil))) {
5236 14 : dbl v = dbl_nil;
5237 14 : bn = BATconstant(min, TYPE_dbl, &v, ngrp, TRANSIENT);
5238 14 : goto doreturn;
5239 : }
5240 :
5241 34 : delta1 = GDKmalloc(ngrp * sizeof(dbl));
5242 34 : delta2 = GDKmalloc(ngrp * sizeof(dbl));
5243 34 : up = GDKmalloc(ngrp * sizeof(dbl));
5244 34 : down1 = GDKmalloc(ngrp * sizeof(dbl));
5245 34 : down2 = GDKmalloc(ngrp * sizeof(dbl));
5246 34 : cnts = GDKzalloc(ngrp * sizeof(BUN));
5247 34 : mean1 = GDKmalloc(ngrp * sizeof(dbl));
5248 34 : mean2 = GDKmalloc(ngrp * sizeof(dbl));
5249 :
5250 34 : if (mean1 == NULL || mean2 == NULL || delta1 == NULL || delta2 == NULL || up == NULL || down1 == NULL || down2 == NULL || cnts == NULL)
5251 0 : goto alloc_fail;
5252 :
5253 307 : TIMEOUT_LOOP_IDX(i, ngrp, timeoffset) {
5254 205 : up[i] = 0;
5255 205 : down1[i] = 0;
5256 205 : down2[i] = 0;
5257 205 : mean1[i] = 0;
5258 205 : mean2[i] = 0;
5259 : }
5260 :
5261 34 : bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
5262 34 : if (bn == NULL)
5263 0 : goto alloc_fail;
5264 34 : dbls = (dbl *) Tloc(bn, 0);
5265 :
5266 34 : if (!g || BATtdense(g))
5267 : gids = NULL;
5268 : else
5269 34 : gids = (const oid *) Tloc(g, 0);
5270 :
5271 34 : b1i = bat_iterator(b1);
5272 34 : b2i = bat_iterator(b2);
5273 34 : switch (b1->ttype) {
5274 8 : case TYPE_bte:
5275 160 : AGGR_CORRELATION(bte);
5276 : break;
5277 0 : case TYPE_sht:
5278 0 : AGGR_CORRELATION(sht);
5279 : break;
5280 21 : case TYPE_int:
5281 791 : AGGR_CORRELATION(int);
5282 : break;
5283 4 : case TYPE_lng:
5284 324 : AGGR_CORRELATION(lng);
5285 : break;
5286 : #ifdef HAVE_HGE
5287 0 : case TYPE_hge:
5288 0 : AGGR_CORRELATION(hge);
5289 : break;
5290 : #endif
5291 0 : case TYPE_flt:
5292 0 : AGGR_CORRELATION(flt);
5293 : break;
5294 1 : case TYPE_dbl:
5295 4 : AGGR_CORRELATION(dbl);
5296 : break;
5297 0 : default:
5298 0 : GDKerror("type (%s) not supported.\n", ATOMname(b1->ttype));
5299 0 : goto bailout;
5300 : }
5301 33 : bat_iterator_end(&b1i);
5302 33 : bat_iterator_end(&b2i);
5303 33 : GDKfree(mean1);
5304 33 : GDKfree(mean2);
5305 33 : GDKfree(delta1);
5306 33 : GDKfree(delta2);
5307 33 : GDKfree(up);
5308 33 : GDKfree(down1);
5309 33 : GDKfree(down2);
5310 33 : GDKfree(cnts);
5311 33 : BATsetcount(bn, ngrp);
5312 33 : bn->tkey = ngrp <= 1;
5313 33 : bn->tsorted = ngrp <= 1;
5314 33 : bn->trevsorted = ngrp <= 1;
5315 33 : bn->tnil = nils != 0;
5316 33 : bn->tnonil = nils == 0;
5317 48 : doreturn:
5318 48 : TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT ",g=" ALGOBATFMT
5319 : ",e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT
5320 : ",skip_nils=%s -> " ALGOOPTBATFMT
5321 : " (" LLFMT " usec)\n",
5322 : ALGOBATPAR(b1), ALGOBATPAR(b2), ALGOBATPAR(g),
5323 : ALGOOPTBATPAR(e), ALGOOPTBATPAR(s),
5324 : skip_nils ? "true" : "false",
5325 : ALGOOPTBATPAR(bn),
5326 : GDKusec() - t0);
5327 : return bn;
5328 1 : overflow:
5329 1 : GDKerror("22003!overflow in calculation.\n");
5330 1 : bailout:
5331 1 : bat_iterator_end(&b1i);
5332 1 : bat_iterator_end(&b2i);
5333 1 : alloc_fail:
5334 1 : BBPreclaim(bn);
5335 1 : GDKfree(mean1);
5336 1 : GDKfree(mean2);
5337 1 : GDKfree(delta1);
5338 1 : GDKfree(delta2);
5339 1 : GDKfree(up);
5340 1 : GDKfree(down1);
5341 1 : GDKfree(down2);
5342 1 : GDKfree(cnts);
5343 1 : return NULL;
5344 : }
|
