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