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 45164 : BATgroupaggrinit(BAT *b, BAT *g, BAT *e, BAT *s,
66 : /* outputs: */
67 : oid *minp, oid *maxp, BUN *ngrpp,
68 : struct canditer *ci)
69 : {
70 45164 : oid min, max;
71 45164 : BUN i, ngrp;
72 45164 : const oid *restrict gids;
73 :
74 45164 : if (b == NULL)
75 : return "b must exist";
76 45164 : canditer_init(ci, b, s);
77 45194 : if (g) {
78 36208 : if (ci->ncand != BATcount(g) ||
79 20508 : (ci->ncand != 0 && ci->seq != g->hseqbase))
80 : return "b with s and g must be aligned";
81 36196 : assert(BATttype(g) == TYPE_oid);
82 : }
83 : if (g == NULL) {
84 : min = 0;
85 : max = 0;
86 : ngrp = 1;
87 36196 : } 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 36193 : ngrp = BATcount(e);
138 36193 : min = e->hseqbase;
139 36193 : max = e->hseqbase + ngrp - 1;
140 : }
141 45182 : *minp = min;
142 45182 : *maxp = max;
143 45182 : *ngrpp = ngrp;
144 :
145 45182 : 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 32775 : twosum(volatile double *hi, volatile double *lo, double x, double y)
163 : {
164 32775 : volatile double yr;
165 :
166 32775 : assert(fabs(x) >= fabs(y));
167 :
168 32775 : *hi = x + y;
169 32775 : yr = *hi - x;
170 32775 : *lo = y - yr;
171 32775 : }
172 :
173 : static inline void
174 14855 : exchange(double *x, double *y)
175 : {
176 14855 : double t = *x;
177 14855 : *x = *y;
178 14855 : *y = t;
179 14855 : }
180 :
181 : /* this function was adapted from https://bugs.python.org/file10357/msum4.py */
182 : BUN
183 1112 : 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 1112 : 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 1112 : BUN listi;
201 1112 : int parti;
202 1112 : int i;
203 1112 : BUN grp;
204 1112 : double x, y;
205 1112 : volatile double lo, hi;
206 1112 : double twopow = pow((double) FLT_RADIX, (double) (DBL_MAX_EXP - 1));
207 1112 : BUN nils = 0;
208 1112 : volatile flt f;
209 :
210 1112 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
211 :
212 : /* we only deal with the two floating point types */
213 1112 : assert(tp1 == TYPE_flt || tp1 == TYPE_dbl);
214 1112 : assert(tp2 == TYPE_flt || tp2 == TYPE_dbl);
215 : /* if input is dbl, then so it output */
216 1112 : assert(tp2 == TYPE_flt || tp1 == TYPE_dbl);
217 : /* if no gids, then we have a single group */
218 1112 : assert(ngrp == 1 || gids != NULL);
219 1112 : if (gids == NULL || ngrp == 1) {
220 1096 : min = max = 0;
221 1096 : ngrp = 1;
222 1096 : gids = NULL;
223 : }
224 1112 : pergroup = GDKmalloc(ngrp * sizeof(*pergroup));
225 1112 : if (pergroup == NULL)
226 : return BUN_NONE;
227 10754 : for (grp = 0; grp < ngrp; grp++) {
228 19284 : pergroup[grp] = (struct pergroup) {
229 : .maxpartials = 2,
230 9642 : .partials = GDKmalloc(2 * sizeof(double)),
231 : };
232 9642 : 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 29294 : TIMEOUT_LOOP(ci->ncand, qry_ctx) {
240 27070 : listi = canditer_next(ci) - seqb;
241 27070 : grp = gids ? gids[listi] : 0;
242 27070 : if (grp < min || grp > max)
243 0 : continue;
244 27070 : if (pergroup[grp].partials == NULL)
245 0 : continue;
246 27070 : if (tp1 == TYPE_flt && !is_flt_nil(((const flt *) values)[listi]))
247 16 : x = ((const flt *) values)[listi];
248 27054 : 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 26983 : pergroup[grp].valseen = true;
265 : #ifdef INFINITES_ALLOWED
266 : if (isinf(x)) {
267 : pergroup[grp].infs += x;
268 : continue;
269 : }
270 : #endif
271 26983 : i = 0;
272 58217 : for (parti = 0; parti < pergroup[grp].npartials; parti++) {
273 31234 : y = pergroup[grp].partials[parti];
274 31234 : if (fabs(x) < fabs(y))
275 14847 : exchange(&x, &y);
276 31234 : twosum(&hi, &lo, x, y);
277 31234 : 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 31234 : if (lo != 0)
287 16428 : pergroup[grp].partials[i++] = lo;
288 31234 : x = hi;
289 : }
290 26983 : if (x != 0) {
291 26912 : 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 26912 : pergroup[grp].partials[i++] = x;
300 : }
301 26983 : pergroup[grp].npartials = i;
302 : }
303 1112 : TIMEOUT_CHECK(qry_ctx, GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx));
304 10726 : for (grp = 0; grp < ngrp; grp++) {
305 9642 : if (pergroup[grp].partials == NULL)
306 0 : continue;
307 9642 : 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 9611 : 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 9605 : if (pergroup[grp].infs != 0)
373 28 : goto overflow;
374 :
375 9577 : 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 9563 : hi = pergroup[grp].partials[--pergroup[grp].npartials];
387 9563 : 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 9563 : 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 9563 : GDKfree(pergroup[grp].partials);
403 9563 : pergroup[grp].partials = NULL;
404 9563 : 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 9552 : } else if (is_dbl_nil(hi)) {
412 0 : goto overflow;
413 : } else {
414 9552 : ((dbl *) results)[grp] = hi;
415 : }
416 : }
417 1084 : GDKfree(pergroup);
418 1084 : 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 11151 : 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 11151 : BUN nils = 0;
716 11151 : BUN i;
717 11151 : oid gid;
718 11151 : unsigned int *restrict seen = NULL; /* bitmask for groups that we've seen */
719 :
720 11151 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
721 :
722 11155 : 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 1094 : if (tp1 != TYPE_flt && tp1 != TYPE_dbl)
729 0 : goto unsupported;
730 1106 : *algo = "sum: floating point";
731 1106 : 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 10049 : seen = GDKzalloc(((ngrp + 31) / 32) * sizeof(int));
738 10048 : if (seen == NULL) {
739 : return BUN_NONE;
740 : }
741 :
742 10048 : switch (tp2) {
743 5 : case TYPE_bte: {
744 5 : bte *restrict sums = (bte *) results;
745 5 : switch (tp1) {
746 5 : case TYPE_bte:
747 20 : 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 783 : case TYPE_int: {
772 783 : int *restrict sums = (int *) results;
773 783 : 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 783 : case TYPE_int:
787 5139 : AGGR_SUM(int, int);
788 : break;
789 0 : default:
790 0 : goto unsupported;
791 : }
792 : break;
793 : }
794 7581 : case TYPE_lng: {
795 7581 : lng *restrict sums = (lng *) results;
796 7581 : 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 110 : case TYPE_int:
821 110 : if (ci->ncand < ((BUN) 1 << ((sizeof(lng) - sizeof(int)) << 3)))
822 481017 : AGGR_SUM_NOOVL(int, lng);
823 : else
824 0 : AGGR_SUM(int, lng);
825 : break;
826 : #endif
827 7471 : case TYPE_lng:
828 8562666 : AGGR_SUM(lng, lng);
829 : break;
830 0 : default:
831 0 : goto unsupported;
832 : }
833 : break;
834 : }
835 : #ifdef HAVE_HGE
836 1654 : case TYPE_hge: {
837 1654 : hge *sums = (hge *) results;
838 1654 : switch (tp1) {
839 168 : case TYPE_bte:
840 4180137 : AGGR_SUM_NOOVL(bte, hge);
841 : break;
842 13 : case TYPE_sht:
843 197 : AGGR_SUM_NOOVL(sht, hge);
844 : break;
845 831 : case TYPE_int:
846 54203194 : AGGR_SUM_NOOVL(int, hge);
847 : break;
848 57 : case TYPE_lng:
849 4632683 : AGGR_SUM_NOOVL(lng, hge);
850 : break;
851 585 : case TYPE_hge:
852 35627468 : 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 10055 : if (nils == 0 && nil_if_empty) {
865 : /* figure out whether there were any empty groups
866 : * (that result in a nil value) */
867 10047 : if (ngrp & 0x1F) {
868 : /* fill last slot */
869 10036 : seen[ngrp >> 5] |= ~0U << (ngrp & 0x1F);
870 : }
871 232005 : for (i = 0, ngrp = (ngrp + 31) / 32; i < ngrp; i++) {
872 222465 : if (seen[i] != ~0U) {
873 : nils = 1;
874 : break;
875 : }
876 : }
877 : }
878 10055 : GDKfree(seen);
879 :
880 10055 : 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 6244 : BATgroupsum(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
901 : {
902 6244 : const oid *restrict gids;
903 6244 : oid min, max;
904 6244 : BUN ngrp;
905 6244 : BUN nils;
906 6244 : BAT *bn;
907 6244 : struct canditer ci;
908 6244 : const char *err;
909 6244 : const char *algo = NULL;
910 6244 : lng t0 = 0;
911 :
912 6244 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
913 :
914 6244 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
915 0 : GDKerror("%s\n", err);
916 0 : return NULL;
917 : }
918 6236 : if (g == NULL) {
919 0 : GDKerror("b and g must be aligned\n");
920 0 : return NULL;
921 : }
922 :
923 6236 : if (ci.ncand == 0 || ngrp == 0) {
924 : /* trivial: no sums, so return bat aligned with g with
925 : * nil in the tail */
926 1886 : return BATconstant(ngrp == 0 ? 0 : min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
927 : }
928 :
929 4350 : if ((e == NULL ||
930 4351 : (BATcount(e) == ci.ncand && e->hseqbase == ci.hseq)) &&
931 1412 : (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 1412 : return BATconvert(b, s, tp, 0, 0, 0);
935 : }
936 :
937 2938 : bn = BATconstant(min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
938 2929 : if (bn == NULL) {
939 : return NULL;
940 : }
941 :
942 2929 : if (BATtdense(g))
943 : gids = NULL;
944 : else
945 2616 : gids = (const oid *) Tloc(g, 0);
946 :
947 2929 : BATiter bi = bat_iterator(b);
948 5885 : nils = dosum(bi.base, bi.nonil, b->hseqbase, &ci,
949 2942 : Tloc(bn, 0), ngrp, bi.type, tp, gids, min, max,
950 : skip_nils, true, __func__, &algo);
951 2943 : bat_iterator_end(&bi);
952 :
953 2937 : if (nils < BUN_NONE) {
954 2937 : BATsetcount(bn, ngrp);
955 2943 : bn->tkey = BATcount(bn) <= 1;
956 2943 : bn->tsorted = BATcount(bn) <= 1;
957 2943 : bn->trevsorted = BATcount(bn) <= 1;
958 2943 : bn->tnil = nils != 0;
959 2943 : bn->tnonil = nils == 0;
960 : } else {
961 0 : BBPunfix(bn->batCacheid);
962 0 : bn = NULL;
963 : }
964 :
965 2943 : if (algo)
966 2943 : MT_thread_setalgorithm(algo);
967 2937 : 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 8646 : BATsum(void *res, int tp, BAT *b, BAT *s, bool skip_nils, bool nil_if_empty)
1019 : {
1020 8646 : oid min, max;
1021 8646 : BUN ngrp;
1022 8646 : struct canditer ci;
1023 8646 : const char *err;
1024 8646 : const char *algo = NULL;
1025 8646 : lng t0 = 0;
1026 :
1027 8646 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
1028 :
1029 8646 : 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 8646 : 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 8616 : switch (tp) {
1085 4 : case TYPE_bte:
1086 4 : * (bte *) res = nil_if_empty ? bte_nil : 0;
1087 4 : break;
1088 17 : case TYPE_sht:
1089 17 : * (sht *) res = nil_if_empty ? sht_nil : 0;
1090 17 : break;
1091 373 : case TYPE_int:
1092 373 : * (int *) res = nil_if_empty ? int_nil : 0;
1093 373 : break;
1094 6511 : case TYPE_lng:
1095 6511 : * (lng *) res = nil_if_empty ? lng_nil : 0;
1096 6511 : break;
1097 : #ifdef HAVE_HGE
1098 503 : case TYPE_hge:
1099 503 : * (hge *) res = nil_if_empty ? hge_nil : 0;
1100 503 : break;
1101 : #endif
1102 1208 : case TYPE_flt:
1103 : case TYPE_dbl:
1104 1208 : 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 1207 : break;
1158 : }
1159 1207 : if (b->ttype == TYPE_dbl)
1160 1194 : * (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 8615 : if (ci.ncand == 0)
1170 : return GDK_SUCCEED;
1171 8214 : BATiter bi = bat_iterator(b);
1172 16461 : BUN nils = dosum(bi.base, bi.nonil, b->hseqbase, &ci,
1173 8231 : res, true, bi.type, tp, &min, min, max,
1174 : skip_nils, nil_if_empty, __func__, &algo);
1175 8230 : bat_iterator_end(&bi);
1176 8227 : if (algo)
1177 8227 : MT_thread_setalgorithm(algo);
1178 8228 : 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 8228 : 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 35 : 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 44 : 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 88 : for (i = 0, ngrp = (ngrp + 31) / 32; i < ngrp; i++) {
1547 44 : if (seen[i] != ~0U) {
1548 : nils = 1;
1549 : break;
1550 : }
1551 : }
1552 : }
1553 44 : GDKfree(seen);
1554 :
1555 44 : 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 289 : BATgroupavg(BAT **bnp, BAT **cntsp, BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, int scale)
1802 : {
1803 289 : const oid *restrict gids;
1804 289 : oid gid;
1805 289 : oid min, max;
1806 289 : BUN i, ngrp;
1807 289 : BUN nils = 0;
1808 289 : lng *restrict rems = NULL;
1809 289 : lng *restrict cnts = NULL;
1810 289 : dbl *restrict dbls;
1811 289 : BAT *bn = NULL, *cn = NULL;
1812 289 : struct canditer ci;
1813 289 : const char *err;
1814 289 : lng t0 = 0;
1815 289 : BATiter bi = {0};
1816 :
1817 289 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
1818 :
1819 289 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
1820 :
1821 289 : assert(tp == TYPE_dbl);
1822 289 : (void) tp; /* compatibility (with other BATgroup*
1823 : * functions) argument */
1824 :
1825 289 : 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 289 : if (g == NULL) {
1830 0 : GDKerror("b and g must be aligned\n");
1831 0 : return GDK_FAIL;
1832 : }
1833 :
1834 289 : 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 280 : if ((!skip_nils || cntsp == NULL || b->tnonil) &&
1854 224 : (e == NULL ||
1855 224 : (BATcount(e) == ci.ncand && e->hseqbase == b->hseqbase)) &&
1856 116 : (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 116 : if ((bn = BATconvert(b, s, TYPE_dbl, 0, 0, 0)) == NULL)
1860 : return GDK_FAIL;
1861 116 : if (cntsp) {
1862 22 : lng one = 1;
1863 22 : if ((cn = BATconstant(ngrp == 0 ? 0 : min, TYPE_lng, &one, ngrp, TRANSIENT)) == NULL) {
1864 0 : BBPreclaim(bn);
1865 0 : return GDK_FAIL;
1866 : }
1867 22 : *cntsp = cn;
1868 : }
1869 116 : *bnp = bn;
1870 116 : return GDK_SUCCEED;
1871 : }
1872 :
1873 : /* allocate temporary space to do per group calculations */
1874 164 : switch (b->ttype) {
1875 138 : 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 138 : rems = GDKzalloc(ngrp * sizeof(lng));
1883 138 : if (rems == NULL)
1884 0 : goto bailout1;
1885 : break;
1886 : default:
1887 : break;
1888 : }
1889 164 : if (cntsp) {
1890 145 : if ((cn = COLnew(min, TYPE_lng, ngrp, TRANSIENT)) == NULL)
1891 0 : goto bailout1;
1892 145 : cnts = (lng *) Tloc(cn, 0);
1893 145 : memset(cnts, 0, ngrp * sizeof(lng));
1894 : } else {
1895 19 : cnts = GDKzalloc(ngrp * sizeof(lng));
1896 19 : if (cnts == NULL)
1897 0 : goto bailout1;
1898 : }
1899 :
1900 164 : bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
1901 164 : if (bn == NULL)
1902 0 : goto bailout1;
1903 164 : dbls = (dbl *) Tloc(bn, 0);
1904 :
1905 164 : if (BATtdense(g))
1906 : gids = NULL;
1907 : else
1908 106 : gids = (const oid *) Tloc(g, 0);
1909 :
1910 164 : bi = bat_iterator(b);
1911 164 : 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 116 : case TYPE_int:
1919 21398373 : AGGR_AVG(int);
1920 116 : 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 335 : 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 164 : bat_iterator_end(&bi);
1940 162 : GDKfree(rems);
1941 163 : if (cn == NULL)
1942 19 : GDKfree(cnts);
1943 : else {
1944 144 : BATsetcount(cn, ngrp);
1945 144 : cn->tkey = BATcount(cn) <= 1;
1946 144 : cn->tsorted = BATcount(cn) <= 1;
1947 144 : cn->trevsorted = BATcount(cn) <= 1;
1948 144 : cn->tnil = false;
1949 144 : cn->tnonil = true;
1950 144 : *cntsp = cn;
1951 : }
1952 163 : 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 163 : BATsetcount(bn, ngrp);
1960 163 : bn->tkey = BATcount(bn) <= 1;
1961 163 : bn->tsorted = BATcount(bn) <= 1;
1962 163 : bn->trevsorted = BATcount(bn) <= 1;
1963 163 : bn->tnil = nils != 0;
1964 163 : bn->tnonil = nils == 0;
1965 163 : *bnp = bn;
1966 163 : 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 417 : BATgroupavg3(BAT **avgp, BAT **remp, BAT **cntp, BAT *b, BAT *g, BAT *e, BAT *s, bool skip_nils)
1997 : {
1998 417 : const char *err;
1999 417 : oid min, max;
2000 417 : BUN ngrp;
2001 417 : struct canditer ci;
2002 417 : BAT *bn, *rn, *cn;
2003 417 : BUN i;
2004 417 : oid o;
2005 :
2006 417 : 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 416 : if (ci.ncand == 0 || ngrp == 0) {
2013 41 : if (ngrp == 0)
2014 29 : min = 0;
2015 41 : 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 375 : ValRecord zero;
2031 375 : (void) VALinit(&zero, TYPE_bte, &(bte){0});
2032 374 : bn = BATconstant(min, b->ttype, VALconvert(b->ttype, &zero),
2033 : ngrp, TRANSIENT);
2034 375 : rn = BATconstant(min, TYPE_lng, &(lng){0}, ngrp, TRANSIENT);
2035 375 : cn = BATconstant(min, TYPE_lng, &(lng){0}, ngrp, TRANSIENT);
2036 374 : 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 374 : lng *rems = Tloc(rn, 0);
2043 374 : lng *cnts = Tloc(cn, 0);
2044 374 : const oid *gids = g && !BATtdense(g) ? Tloc(g, 0) : NULL;
2045 374 : oid gid = ngrp == 1 && gids ? gids[0] - min : 0;
2046 :
2047 374 : BATiter bi = bat_iterator(b);
2048 :
2049 740 : 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 268 : TIMEOUT_LOOP(ci.ncand, qry_ctx) {
2054 238 : o = canditer_next(&ci) - b->hseqbase;
2055 237 : if (ngrp > 1)
2056 0 : gid = gids ? gids[o] - min : o;
2057 237 : 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 237 : } else if (!is_lng_nil(cnts[gid])) {
2067 237 : AVERAGE_ITER(bte, vals[o], avgs[gid], rems[gid], cnts[gid]);
2068 : }
2069 : }
2070 46 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) {
2071 15 : 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 15 : if (!is_bte_nil(avgs[i]) && rems[i] > 0) {
2082 14 : 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 14 : 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 241 : case TYPE_int: {
2147 241 : const int *vals = (const int *) bi.base;
2148 241 : int *avgs = Tloc(bn, 0);
2149 1833988 : TIMEOUT_LOOP(ci.ncand, qry_ctx) {
2150 1833164 : o = canditer_next(&ci) - b->hseqbase;
2151 1833163 : if (ngrp > 1)
2152 205041 : gid = gids ? gids[o] - min : o;
2153 1833163 : if (is_int_nil(vals[o])) {
2154 120025 : 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 1713138 : } else if (!is_lng_nil(cnts[gid])) {
2163 1853336 : AVERAGE_ITER(int, vals[o], avgs[gid], rems[gid], cnts[gid]);
2164 : }
2165 : }
2166 188593 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) {
2167 188105 : if (cnts[i] == 0) {
2168 848 : avgs[i] = int_nil;
2169 848 : 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 187257 : if (!is_int_nil(avgs[i]) && rems[i] > 0) {
2178 70724 : if (avgs[i] < 0) {
2179 54919 : if (2*rems[i] > cnts[i]) {
2180 20767 : avgs[i]++;
2181 20767 : rems[i] -= cnts[i];
2182 : }
2183 : } else {
2184 15805 : if (2*rems[i] >= cnts[i]) {
2185 14670 : avgs[i]++;
2186 14670 : rems[i] -= cnts[i];
2187 : }
2188 : }
2189 : }
2190 : #endif
2191 : }
2192 : break;
2193 : }
2194 102 : case TYPE_lng: {
2195 102 : const lng *vals = (const lng *) bi.base;
2196 102 : lng *avgs = Tloc(bn, 0);
2197 198886 : TIMEOUT_LOOP(ci.ncand, qry_ctx) {
2198 198681 : o = canditer_next(&ci) - b->hseqbase;
2199 198682 : if (ngrp > 1)
2200 176644 : gid = gids ? gids[o] - min : o;
2201 198682 : if (is_lng_nil(vals[o])) {
2202 209 : 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 198473 : } else if (!is_lng_nil(cnts[gid])) {
2211 199212 : AVERAGE_ITER(lng, vals[o], avgs[gid], rems[gid], cnts[gid]);
2212 : }
2213 : }
2214 61963 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) {
2215 61755 : if (cnts[i] == 0) {
2216 165 : avgs[i] = lng_nil;
2217 165 : 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 61590 : if (!is_lng_nil(avgs[i]) && rems[i] > 0) {
2226 846 : 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 717 : if (2*rems[i] >= cnts[i]) {
2233 716 : avgs[i]++;
2234 716 : 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 375 : bat_iterator_end(&bi);
2294 375 : TIMEOUT_CHECK(qry_ctx, GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx));
2295 375 : BATsetcount(bn, ngrp);
2296 375 : BATsetcount(rn, ngrp);
2297 375 : BATsetcount(cn, ngrp);
2298 376 : bn->tnonil = !bn->tnil;
2299 376 : rn->tnonil = !rn->tnil;
2300 376 : cn->tnonil = !cn->tnil;
2301 376 : bn->tkey = rn->tkey = cn->tkey = ngrp == 1;
2302 376 : bn->tsorted = rn->tsorted = cn->tsorted = ngrp == 1;
2303 376 : bn->trevsorted = rn->trevsorted = cn->trevsorted = ngrp == 1;
2304 376 : *avgp = bn;
2305 376 : *remp = rn;
2306 376 : *cntp = cn;
2307 376 : 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 190639 : 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 190639 : lng cnt = cnt1 + cnt2;
2454 :
2455 190639 : if (rem2 < 0) {
2456 39683 : avg2--;
2457 39683 : rem2 += cnt2;
2458 : }
2459 190639 : *cntp = cnt;
2460 : #ifdef HAVE_HGE
2461 190639 : hge v = avg1 * cnt1 + rem1 + avg2 * cnt2 + rem2;
2462 190639 : int a = (int) (v / cnt);
2463 190639 : v %= cnt;
2464 190639 : if (v < 0) {
2465 136668 : a--;
2466 136668 : v += cnt;
2467 : }
2468 190639 : *avgp = a;
2469 190639 : *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 190639 : }
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 40 : 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 82 : 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 34 : case TYPE_int: {
2764 34 : const int *vals = (const int *) bi.base;
2765 34 : int *avgs = Tloc(bn, 0);
2766 191741 : TIMEOUT_LOOP_IDX(i, ci.ncand, qry_ctx) {
2767 191625 : if (ngrp > 1)
2768 191458 : gid = gids ? gids[i] - min : i;
2769 191625 : if (is_int_nil(vals[i])) {
2770 986 : if (!skip_nils) {
2771 0 : avgs[gid] = int_nil;
2772 0 : bn->tnil = true;
2773 : }
2774 190639 : } else if (!is_int_nil(avgs[gid])) {
2775 190639 : combine_averages_int(&avgs[gid], &rems[gid],
2776 : &cnts[gid], avgs[gid],
2777 190639 : rems[gid], cnts[gid],
2778 190639 : vals[i], orems[i],
2779 190639 : ocnts[i]);
2780 : }
2781 : }
2782 57976 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) {
2783 57906 : if (cnts[i] == 0) {
2784 67 : avgs[i] = int_nil;
2785 67 : 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 57839 : if (!is_int_nil(avgs[i]) && rems[i] > 0) {
2792 37295 : if (avgs[i] < 0) {
2793 33434 : if (2*rems[i] > cnts[i])
2794 16188 : avgs[i]++;
2795 : } else {
2796 3861 : if (2*rems[i] >= cnts[i])
2797 3089 : 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 7971 : BATcalcavg(BAT *b, BAT *s, dbl *avg, BUN *vals, int scale)
2988 : {
2989 7971 : lng n = 0, r = 0;
2990 7971 : BUN i = 0;
2991 : #ifdef HAVE_HGE
2992 7971 : hge sum = 0;
2993 : #else
2994 : lng sum = 0;
2995 : #endif
2996 7971 : struct canditer ci;
2997 7971 : const void *restrict src;
2998 :
2999 7971 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
3000 :
3001 7975 : canditer_init(&ci, b, s);
3002 :
3003 7973 : BATiter bi = bat_iterator(b);
3004 7977 : src = bi.base;
3005 :
3006 7977 : switch (b->ttype) {
3007 8 : case TYPE_bte:
3008 24 : AVERAGE_TYPE(bte);
3009 : break;
3010 0 : case TYPE_sht:
3011 0 : AVERAGE_TYPE(sht);
3012 : break;
3013 7913 : case TYPE_int:
3014 4750123 : 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 7974 : bat_iterator_end(&bi);
3036 7976 : if (scale != 0 && !is_dbl_nil(*avg))
3037 0 : *avg /= pow(10.0, (double) scale);
3038 7976 : if (vals)
3039 7976 : *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 28171 : BATgroupcount(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
3070 : {
3071 28171 : const oid *restrict gids;
3072 28171 : oid gid;
3073 28171 : oid min, max;
3074 28171 : BUN i, ngrp;
3075 28171 : lng *restrict cnts;
3076 28171 : BAT *bn = NULL;
3077 28171 : int t;
3078 28171 : const void *nil;
3079 28171 : int (*atomcmp)(const void *, const void *);
3080 28171 : struct canditer ci;
3081 28171 : const char *err;
3082 28171 : lng t0 = 0;
3083 28171 : BATiter bi = {0};
3084 :
3085 28171 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
3086 :
3087 28173 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
3088 :
3089 28173 : assert(tp == TYPE_lng);
3090 28173 : (void) tp; /* compatibility (with other BATgroup* */
3091 :
3092 28173 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
3093 0 : GDKerror("%s\n", err);
3094 0 : return NULL;
3095 : }
3096 28177 : if (g == NULL) {
3097 0 : GDKerror("b and g must be aligned\n");
3098 0 : return NULL;
3099 : }
3100 :
3101 28177 : if (ci.ncand == 0 || ngrp == 0) {
3102 : /* trivial: no counts, so return bat aligned with g
3103 : * with zero in the tail */
3104 13615 : lng zero = 0;
3105 13615 : return BATconstant(ngrp == 0 ? 0 : min, TYPE_lng, &zero, ngrp, TRANSIENT);
3106 : }
3107 :
3108 14562 : bn = COLnew(min, TYPE_lng, ngrp, TRANSIENT);
3109 14558 : if (bn == NULL)
3110 : return NULL;
3111 14558 : cnts = (lng *) Tloc(bn, 0);
3112 14558 : memset(cnts, 0, ngrp * sizeof(lng));
3113 :
3114 14558 : if (BATtdense(g))
3115 : gids = NULL;
3116 : else
3117 11102 : gids = (const oid *) Tloc(g, 0);
3118 :
3119 14558 : 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 14469 : if (gids) {
3123 18095079 : TIMEOUT_LOOP(ci.ncand, qry_ctx) {
3124 18071872 : i = canditer_next(&ci) - b->hseqbase;
3125 18071874 : if (gids[i] >= min && gids[i] <= max)
3126 18101028 : cnts[gids[i] - min]++;
3127 : }
3128 : } else {
3129 101967 : TIMEOUT_LOOP(ci.ncand, qry_ctx) {
3130 95192 : i = canditer_next(&ci) - b->hseqbase;
3131 95189 : cnts[i] = 1;
3132 : }
3133 : }
3134 : } else {
3135 89 : t = b->ttype;
3136 89 : nil = ATOMnilptr(t);
3137 89 : atomcmp = ATOMcompare(t);
3138 89 : t = ATOMbasetype(t);
3139 :
3140 89 : bi = bat_iterator(b);
3141 :
3142 91 : 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 72 : case TYPE_int:
3150 15590 : 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 11 : default:
3167 129 : TIMEOUT_LOOP(ci.ncand, qry_ctx) {
3168 107 : i = canditer_next(&ci) - b->hseqbase;
3169 107 : if (gids == NULL ||
3170 105 : (gids[i] >= min && gids[i] <= max)) {
3171 105 : if (gids)
3172 105 : gid = gids[i] - min;
3173 : else
3174 : gid = (oid) i;
3175 107 : if ((*atomcmp)(BUNtail(bi, i), nil) != 0) {
3176 77 : cnts[gid]++;
3177 : }
3178 : }
3179 : }
3180 : break;
3181 : }
3182 91 : bat_iterator_end(&bi);
3183 : }
3184 14553 : TIMEOUT_CHECK(qry_ctx, GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx));
3185 14550 : BATsetcount(bn, ngrp);
3186 14549 : bn->tkey = BATcount(bn) <= 1;
3187 14549 : bn->tsorted = BATcount(bn) <= 1;
3188 14549 : bn->trevsorted = BATcount(bn) <= 1;
3189 14549 : bn->tnil = false;
3190 14549 : bn->tnonil = true;
3191 14549 : 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 488 : 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 488 : oid gid = 0;
3252 488 : BUN i, nils;
3253 488 : int t;
3254 488 : const void *nil;
3255 488 : int (*atomcmp)(const void *, const void *);
3256 :
3257 488 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
3258 :
3259 489 : nils = ngrp;
3260 55749 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx)
3261 54771 : oids[i] = oid_nil;
3262 490 : if (ci->ncand == 0)
3263 : return nils;
3264 :
3265 490 : t = bi->b->ttype;
3266 490 : nil = ATOMnilptr(t);
3267 490 : atomcmp = ATOMcompare(t);
3268 490 : t = ATOMbasetype(t);
3269 490 : oid hseq = bi->b->hseqbase;
3270 :
3271 490 : 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 328 : case TYPE_int:
3279 58849 : AGGR_CMP(int, LT);
3280 : break;
3281 57 : case TYPE_lng:
3282 91545 : 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 22202 : TIMEOUT_LOOP(ci->ncand, qry_ctx) {
3335 22122 : i = canditer_next(ci) - hseq;
3336 22121 : if (gids == NULL ||
3337 43 : (gids[i] >= min && gids[i] <= max)) {
3338 22121 : const void *v = BUNtail(*bi, i);
3339 22122 : if (gids)
3340 44 : gid = gids[i] - min;
3341 44243 : if (!skip_nils ||
3342 22121 : (*atomcmp)(v, nil) != 0) {
3343 19708 : if (is_oid_nil(oids[gid])) {
3344 58 : oids[gid] = i + hseq;
3345 58 : nils--;
3346 19650 : } else if (t != TYPE_void) {
3347 19650 : const void *g = BUNtail(*bi, (BUN) (oids[gid] - hseq));
3348 39300 : if ((*atomcmp)(g, nil) != 0 &&
3349 39300 : ((*atomcmp)(v, nil) == 0 ||
3350 19650 : LT((*atomcmp)(v, g), 0)))
3351 85 : oids[gid] = i + hseq;
3352 : }
3353 : }
3354 : }
3355 : }
3356 : }
3357 : break;
3358 : }
3359 489 : 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 540 : 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 540 : oid gid = 0;
3374 540 : BUN i, nils;
3375 540 : int t;
3376 540 : const void *nil;
3377 540 : int (*atomcmp)(const void *, const void *);
3378 :
3379 540 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
3380 :
3381 540 : nils = ngrp;
3382 815286 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx)
3383 814133 : oids[i] = oid_nil;
3384 540 : if (ci->ncand == 0)
3385 : return nils;
3386 :
3387 540 : t = bi->b->ttype;
3388 540 : nil = ATOMnilptr(t);
3389 540 : atomcmp = ATOMcompare(t);
3390 540 : t = ATOMbasetype(t);
3391 540 : oid hseq = bi->b->hseqbase;
3392 :
3393 540 : switch (t) {
3394 24 : case TYPE_bte:
3395 158 : AGGR_CMP(bte, GT);
3396 : break;
3397 1 : case TYPE_sht:
3398 169 : AGGR_CMP(sht, GT);
3399 : break;
3400 365 : case TYPE_int:
3401 126098 : AGGR_CMP(int, GT);
3402 : break;
3403 73 : case TYPE_lng:
3404 93011 : AGGR_CMP(lng, GT);
3405 : break;
3406 : #ifdef HAVE_HGE
3407 5 : case TYPE_hge:
3408 34770551 : 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 279 : 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 21260 : TIMEOUT_LOOP(ci->ncand, qry_ctx) {
3456 21208 : i = canditer_next(ci) - hseq;
3457 21208 : if (gids == NULL ||
3458 20 : (gids[i] >= min && gids[i] <= max)) {
3459 21208 : const void *v = BUNtail(*bi, i);
3460 21208 : if (gids)
3461 20 : gid = gids[i] - min;
3462 42416 : if (!skip_nils ||
3463 21208 : (*atomcmp)(v, nil) != 0) {
3464 18792 : if (is_oid_nil(oids[gid])) {
3465 31 : oids[gid] = i + hseq;
3466 31 : nils--;
3467 : } else {
3468 18761 : const void *g = BUNtail(*bi, (BUN) (oids[gid] - hseq));
3469 18761 : if (t == TYPE_void ||
3470 37522 : ((*atomcmp)(g, nil) != 0 &&
3471 37522 : ((*atomcmp)(v, nil) == 0 ||
3472 18761 : GT((*atomcmp)(v, g), 0))))
3473 206 : oids[gid] = i + hseq;
3474 : }
3475 : }
3476 : }
3477 : }
3478 : }
3479 : break;
3480 : }
3481 539 : TIMEOUT_CHECK(qry_ctx, TIMEOUT_HANDLER(BUN_NONE, qry_ctx));
3482 :
3483 : return nils;
3484 : }
3485 :
3486 : static BAT *
3487 987 : 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 987 : const oid *restrict gids;
3494 987 : oid min, max;
3495 987 : BUN ngrp;
3496 987 : oid *restrict oids;
3497 987 : BAT *bn = NULL;
3498 987 : BUN nils;
3499 987 : struct canditer ci;
3500 987 : const char *err;
3501 987 : lng t0 = 0;
3502 :
3503 987 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
3504 :
3505 987 : assert(tp == TYPE_oid);
3506 987 : (void) tp; /* compatibility (with other BATgroup* */
3507 :
3508 987 : 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 987 : 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 987 : if (ci.ncand == 0 || ngrp == 0) {
3520 : /* trivial: no minimums, so return bat aligned with g
3521 : * with nil in the tail */
3522 126 : return BATconstant(ngrp == 0 ? 0 : min, TYPE_oid, &oid_nil, ngrp, TRANSIENT);
3523 : }
3524 :
3525 861 : bn = COLnew(min, TYPE_oid, ngrp, TRANSIENT);
3526 859 : if (bn == NULL)
3527 : return NULL;
3528 859 : oids = (oid *) Tloc(bn, 0);
3529 :
3530 859 : if (g == NULL || BATtdense(g))
3531 : gids = NULL;
3532 : else
3533 467 : gids = (const oid *) Tloc(g, 0);
3534 :
3535 859 : BATiter bi = bat_iterator(b);
3536 862 : nils = (*minmax)(oids, &bi, gids, ngrp, min, max, &ci,
3537 862 : skip_nils, g && BATtdense(g));
3538 860 : bat_iterator_end(&bi);
3539 862 : if (nils == BUN_NONE) {
3540 0 : BBPreclaim(bn);
3541 0 : return NULL;
3542 : }
3543 :
3544 862 : BATsetcount(bn, ngrp);
3545 :
3546 862 : bn->tkey = BATcount(bn) <= 1;
3547 862 : bn->tsorted = BATcount(bn) <= 1;
3548 862 : bn->trevsorted = BATcount(bn) <= 1;
3549 862 : bn->tnil = nils != 0;
3550 862 : bn->tnonil = nils == 0;
3551 862 : 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 438 : BATgroupmin(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
3562 : {
3563 438 : 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 1140 : BATmin_skipnil(BAT *b, void *aggr, bit skipnil)
3571 : {
3572 1140 : const void *res = NULL;
3573 1140 : size_t s;
3574 1140 : lng t0 = 0;
3575 :
3576 1140 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
3577 :
3578 1140 : 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 1140 : BATiter bi = bat_iterator(b);
3585 1145 : if (bi.count == 0) {
3586 283 : res = ATOMnilptr(bi.type);
3587 862 : } else if (bi.minpos != BUN_NONE) {
3588 343 : res = BUNtail(bi, bi.minpos);
3589 : } else {
3590 519 : oid pos;
3591 770 : BAT *pb = BATdescriptor(VIEWtparent(b));
3592 519 : Heap *oidxh = NULL;
3593 519 : bool usepoidx = false;
3594 :
3595 519 : if (BATordered(b)) {
3596 318 : if (skipnil && !bi.nonil) {
3597 54 : pos = binsearch(NULL, 0, bi.type, bi.base,
3598 54 : bi.vh ? bi.vh->base : NULL,
3599 54 : bi.width, 0, bi.count,
3600 54 : ATOMnilptr(bi.type), 1, 1);
3601 54 : if (pos == bi.count)
3602 42 : pos = oid_nil;
3603 : else
3604 12 : pos += b->hseqbase;
3605 : } else {
3606 264 : pos = b->hseqbase;
3607 : }
3608 201 : } 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 77 : 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 77 : if (oidxh == NULL &&
3630 113 : pb != NULL &&
3631 36 : 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 77 : 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 77 : struct canditer ci;
3675 77 : canditer_init(&ci, b, NULL);
3676 77 : (void) do_groupmin(&pos, &bi, NULL, 1, 0, 0, &ci,
3677 : skipnil, false);
3678 : }
3679 : }
3680 518 : if (is_oid_nil(pos)) {
3681 42 : res = ATOMnilptr(bi.type);
3682 : } else {
3683 476 : bi.minpos = pos - b->hseqbase;
3684 476 : res = BUNtail(bi, bi.minpos);
3685 476 : MT_lock_set(&b->theaplock);
3686 477 : if (bi.count == BATcount(b) && bi.h == b->theap)
3687 477 : b->tminpos = bi.minpos;
3688 477 : MT_lock_unset(&b->theaplock);
3689 477 : if (pb) {
3690 265 : MT_lock_set(&pb->theaplock);
3691 265 : if (bi.count == BATcount(pb) &&
3692 45 : bi.h == pb->theap)
3693 45 : pb->tminpos = bi.minpos;
3694 265 : MT_lock_unset(&pb->theaplock);
3695 : }
3696 : }
3697 787 : BBPreclaim(pb);
3698 : }
3699 1145 : if (aggr == NULL) {
3700 527 : s = ATOMlen(bi.type, res);
3701 527 : aggr = GDKmalloc(s);
3702 : } else {
3703 1236 : s = ATOMsize(ATOMtype(bi.type));
3704 : }
3705 1145 : if (aggr != NULL) /* else: malloc error */
3706 1145 : memcpy(aggr, res, s);
3707 1145 : bat_iterator_end(&bi);
3708 1144 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",skipnil=%d; (" LLFMT " usec)\n",
3709 : ALGOBATPAR(b), skipnil, GDKusec() - t0);
3710 : return aggr;
3711 : }
3712 :
3713 : void *
3714 485 : BATmin(BAT *b, void *aggr)
3715 : {
3716 485 : return BATmin_skipnil(b, aggr, 1);
3717 : }
3718 :
3719 : BAT *
3720 551 : BATgroupmax(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
3721 : {
3722 551 : return BATgroupminmax(b, g, e, s, tp, skip_nils,
3723 : do_groupmax, __func__);
3724 : }
3725 :
3726 : void *
3727 1353 : BATmax_skipnil(BAT *b, void *aggr, bit skipnil)
3728 : {
3729 1353 : const void *res = NULL;
3730 1353 : size_t s;
3731 1353 : lng t0 = 0;
3732 :
3733 1353 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
3734 :
3735 1353 : if (!ATOMlinear(b->ttype)) {
3736 : /* there is no such thing as a largest value if you
3737 : * can't compare values */
3738 0 : GDKerror("non-linear type");
3739 0 : return NULL;
3740 : }
3741 1353 : BATiter bi = bat_iterator(b);
3742 1356 : if (bi.count == 0) {
3743 221 : res = ATOMnilptr(bi.type);
3744 1135 : } else if (bi.maxpos != BUN_NONE) {
3745 743 : res = BUNtail(bi, bi.maxpos);
3746 : } else {
3747 392 : oid pos;
3748 550 : BAT *pb = BATdescriptor(VIEWtparent(b));
3749 392 : Heap *oidxh = NULL;
3750 392 : bool usepoidx = false;
3751 :
3752 392 : if (BATordered(b)) {
3753 260 : pos = bi.count - 1 + b->hseqbase;
3754 304 : if (skipnil && !bi.nonil &&
3755 44 : ATOMcmp(bi.type, BUNtail(bi, bi.count - 1),
3756 : ATOMnilptr(bi.type)) == 0)
3757 34 : pos = oid_nil; /* no non-nil values */
3758 130 : } else if (BATordered_rev(b)) {
3759 39 : pos = b->hseqbase;
3760 73 : if (skipnil && !bi.nonil &&
3761 34 : ATOMcmp(bi.type, BUNtail(bi, 0),
3762 : ATOMnilptr(bi.type)) == 0)
3763 0 : pos = oid_nil; /* no non-nil values */
3764 : } else {
3765 91 : if (BATcheckorderidx(b)) {
3766 0 : MT_lock_set(&b->batIdxLock);
3767 0 : oidxh = b->torderidx;
3768 0 : if (oidxh != NULL)
3769 0 : HEAPincref(oidxh);
3770 0 : MT_lock_unset(&b->batIdxLock);
3771 : }
3772 91 : if (oidxh == NULL &&
3773 129 : pb != NULL &&
3774 38 : BATcheckorderidx(pb)) {
3775 : /* no lock on b needed since it's a view */
3776 0 : MT_lock_set(&pb->batIdxLock);
3777 0 : MT_lock_set(&pb->theaplock);
3778 0 : if (pb->tbaseoff == bi.baseoff &&
3779 0 : BATcount(pb) == bi.count &&
3780 0 : pb->hseqbase == b->hseqbase &&
3781 0 : (oidxh = pb->torderidx) != NULL) {
3782 0 : HEAPincref(oidxh);
3783 0 : usepoidx = true;
3784 : }
3785 0 : MT_lock_unset(&pb->batIdxLock);
3786 0 : MT_lock_unset(&pb->theaplock);
3787 : }
3788 91 : if (oidxh != NULL) {
3789 0 : const oid *ords = (const oid *) oidxh->base + ORDERIDXOFF;
3790 :
3791 0 : MT_thread_setalgorithm(usepoidx ? "using parent oidx" : "using oids");
3792 0 : pos = ords[bi.count - 1];
3793 : /* nils are first, ie !skipnil, check for nils */
3794 0 : if (!skipnil) {
3795 0 : BUN z = ords[0];
3796 :
3797 0 : res = BUNtail(bi, z - b->hseqbase);
3798 :
3799 0 : if (ATOMcmp(bi.type, res, ATOMnilptr(bi.type)) == 0)
3800 0 : pos = z;
3801 : }
3802 0 : HEAPdecref(oidxh, false);
3803 : } else {
3804 91 : struct canditer ci;
3805 91 : canditer_init(&ci, b, NULL);
3806 91 : (void) do_groupmax(&pos, &bi, NULL, 1, 0, 0, &ci,
3807 : skipnil, false);
3808 : }
3809 : }
3810 389 : if (is_oid_nil(pos)) {
3811 34 : res = ATOMnilptr(bi.type);
3812 : } else {
3813 355 : bi.maxpos = pos - b->hseqbase;
3814 355 : res = BUNtail(bi, bi.maxpos);
3815 355 : MT_lock_set(&b->theaplock);
3816 357 : if (bi.count == BATcount(b) && bi.h == b->theap)
3817 358 : b->tmaxpos = bi.maxpos;
3818 357 : MT_lock_unset(&b->theaplock);
3819 358 : if (pb) {
3820 232 : MT_lock_set(&pb->theaplock);
3821 232 : if (bi.count == BATcount(pb) &&
3822 64 : bi.h == pb->theap)
3823 64 : pb->tmaxpos = bi.maxpos;
3824 232 : MT_lock_unset(&pb->theaplock);
3825 : }
3826 : }
3827 626 : BBPreclaim(pb);
3828 : }
3829 1356 : if (aggr == NULL) {
3830 551 : s = ATOMlen(bi.type, res);
3831 552 : aggr = GDKmalloc(s);
3832 : } else {
3833 1609 : s = ATOMsize(ATOMtype(bi.type));
3834 : }
3835 1358 : if (aggr != NULL) /* else: malloc error */
3836 1357 : memcpy(aggr, res, s);
3837 1358 : bat_iterator_end(&bi);
3838 1356 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",skipnil=%d; (" LLFMT " usec)\n",
3839 : ALGOBATPAR(b), skipnil, GDKusec() - t0);
3840 : return aggr;
3841 : }
3842 :
3843 : void *
3844 799 : BATmax(BAT *b, void *aggr)
3845 : {
3846 799 : return BATmax_skipnil(b, aggr, 1);
3847 : }
3848 :
3849 :
3850 : /* ---------------------------------------------------------------------- */
3851 : /* quantiles/median */
3852 :
3853 : #if SIZEOF_OID == SIZEOF_INT
3854 : #define binsearch_oid(indir, offset, vals, lo, hi, v, ordering, last) binsearch_int(indir, offset, (const int *) vals, lo, hi, (int) (v), ordering, last)
3855 : #endif
3856 : #if SIZEOF_OID == SIZEOF_LNG
3857 : #define binsearch_oid(indir, offset, vals, lo, hi, v, ordering, last) binsearch_lng(indir, offset, (const lng *) vals, lo, hi, (lng) (v), ordering, last)
3858 : #endif
3859 :
3860 : #define DO_QUANTILE_AVG(TPE) \
3861 : do { \
3862 : BUN idxlo, idxhi; \
3863 : if (ords) { \
3864 : idxlo = ords[r + (BUN) lo] - b->hseqbase; \
3865 : idxhi = ords[r + (BUN) hi] - b->hseqbase; \
3866 : } else { \
3867 : idxlo = r + (BUN) lo; \
3868 : idxhi = r + (BUN) hi; \
3869 : } \
3870 : TPE low = *(TPE*) BUNtloc(bi, idxhi); \
3871 : TPE high = *(TPE*) BUNtloc(bi, idxlo); \
3872 : if (is_##TPE##_nil(low) || is_##TPE##_nil(high)) { \
3873 : val = dbl_nil; \
3874 : nils++; \
3875 : } else { \
3876 : val = (f - lo) * low + (lo + 1 - f) * high; \
3877 : } \
3878 : } while (0)
3879 :
3880 : static BAT *
3881 66 : doBATgroupquantile(BAT *b, BAT *g, BAT *e, BAT *s, int tp, double quantile,
3882 : bool skip_nils, bool average)
3883 : {
3884 66 : BAT *origb = b;
3885 66 : BAT *origg = g;
3886 66 : oid min, max;
3887 66 : BUN ngrp;
3888 66 : BUN nils = 0;
3889 66 : BAT *bn = NULL;
3890 66 : struct canditer ci;
3891 66 : BAT *t1, *t2;
3892 66 : BATiter bi = {0};
3893 66 : const void *v;
3894 66 : const void *nil = ATOMnilptr(tp);
3895 66 : const void *dnil = nil;
3896 66 : dbl val; /* only used for average */
3897 66 : int (*atomcmp)(const void *, const void *) = ATOMcompare(tp);
3898 66 : const char *err;
3899 66 : lng t0 = 0;
3900 :
3901 66 : size_t counter = 0;
3902 66 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
3903 :
3904 66 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
3905 :
3906 66 : if (average) {
3907 12 : switch (ATOMbasetype(b->ttype)) {
3908 : case TYPE_bte:
3909 : case TYPE_sht:
3910 : case TYPE_int:
3911 : case TYPE_lng:
3912 : #ifdef HAVE_HGE
3913 : case TYPE_hge:
3914 : #endif
3915 : case TYPE_flt:
3916 : case TYPE_dbl:
3917 : break;
3918 0 : default:
3919 0 : GDKerror("incompatible type\n");
3920 0 : return NULL;
3921 : }
3922 : dnil = &dbl_nil;
3923 : }
3924 66 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
3925 0 : GDKerror("%s\n", err);
3926 0 : return NULL;
3927 : }
3928 66 : assert(tp == b->ttype);
3929 66 : if (!ATOMlinear(tp)) {
3930 0 : GDKerror("cannot determine quantile on "
3931 : "non-linear type %s\n", ATOMname(tp));
3932 0 : return NULL;
3933 : }
3934 66 : if (quantile < 0 || quantile > 1) {
3935 2 : GDKerror("cannot determine quantile for "
3936 : "p=%f (p has to be in [0,1])\n", quantile);
3937 2 : return NULL;
3938 : }
3939 :
3940 64 : if (BATcount(b) == 0 || ngrp == 0 || is_dbl_nil(quantile)) {
3941 : /* trivial: no values, thus also no quantiles,
3942 : * so return bat aligned with e with nil in the tail
3943 : * The same happens for a NULL quantile */
3944 24 : return BATconstant(ngrp == 0 ? 0 : min, average ? TYPE_dbl : tp, dnil, ngrp, TRANSIENT);
3945 : }
3946 :
3947 52 : if (s) {
3948 : /* there is a candidate list, replace b (and g, if
3949 : * given) with just the values we're interested in */
3950 0 : b = BATproject(s, b);
3951 0 : if (b == NULL)
3952 : return NULL;
3953 0 : if (g) {
3954 0 : g = BATproject(s, g);
3955 0 : if (g == NULL)
3956 0 : goto bailout;
3957 : }
3958 : }
3959 :
3960 : /* we want to sort b so that we can figure out the quantile, but
3961 : * if g is given, sort g and subsort b so that we can get the
3962 : * quantile for each group */
3963 52 : if (g) {
3964 18 : const oid *restrict grps;
3965 18 : oid prev;
3966 18 : BUN p, q, r;
3967 :
3968 18 : if (BATtdense(g)) {
3969 : /* singleton groups, so calculating quantile is
3970 : * easy */
3971 1 : if (average)
3972 0 : bn = BATconvert(b, NULL, TYPE_dbl, 0, 0, 0);
3973 : else
3974 1 : bn = COLcopy(b, tp, false, TRANSIENT);
3975 1 : BAThseqbase(bn, g->tseqbase); /* deals with NULL */
3976 1 : if (b != origb)
3977 0 : BBPunfix(b->batCacheid);
3978 1 : if (g != origg)
3979 0 : BBPunfix(g->batCacheid);
3980 1 : return bn;
3981 : }
3982 17 : if (BATsort(&t1, &t2, NULL, g, NULL, NULL, false, false, false) != GDK_SUCCEED)
3983 0 : goto bailout;
3984 17 : if (g != origg)
3985 0 : BBPunfix(g->batCacheid);
3986 17 : g = t1;
3987 :
3988 17 : if (BATsort(&t1, NULL, NULL, b, t2, g, false, false, false) != GDK_SUCCEED) {
3989 0 : BBPunfix(t2->batCacheid);
3990 0 : goto bailout;
3991 : }
3992 17 : if (b != origb)
3993 0 : BBPunfix(b->batCacheid);
3994 17 : b = t1;
3995 17 : BBPunfix(t2->batCacheid);
3996 :
3997 17 : if (average)
3998 2 : bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
3999 : else
4000 15 : bn = COLnew(min, tp, ngrp, TRANSIENT);
4001 17 : if (bn == NULL)
4002 0 : goto bailout;
4003 :
4004 17 : bi = bat_iterator(b);
4005 :
4006 17 : grps = (const oid *) Tloc(g, 0);
4007 : /* for each group (r and p are the beginning and end
4008 : * of the current group, respectively) */
4009 70 : for (r = 0, q = BATcount(g); r < q; r = p) {
4010 53 : GDK_CHECK_TIMEOUT(qry_ctx, counter,
4011 : GOTO_LABEL_TIMEOUT_HANDLER(bunins_failed, qry_ctx));
4012 53 : BUN qindex;
4013 53 : prev = grps[r];
4014 : /* search for end of current group (grps is
4015 : * sorted so we can use binary search) */
4016 53 : p = binsearch_oid(NULL, 0, grps, r, q - 1, prev, 1, 1);
4017 53 : if (skip_nils && !bi.nonil) {
4018 : /* within group, locate start of non-nils */
4019 21 : r = binsearch(NULL, 0, tp, bi.base,
4020 21 : bi.vh ? bi.vh->base : NULL,
4021 21 : bi.width, r, p, nil,
4022 : 1, 1);
4023 : }
4024 53 : if (r == p) {
4025 : /* no non-nils found */
4026 11 : v = dnil;
4027 11 : nils++;
4028 42 : } else if (average) {
4029 4 : double f = (p - r - 1) * quantile;
4030 4 : double lo = floor(f);
4031 4 : double hi = ceil(f);
4032 4 : const oid *const ords = NULL;
4033 8 : switch (ATOMbasetype(tp)) {
4034 : case TYPE_bte:
4035 4 : DO_QUANTILE_AVG(bte);
4036 : break;
4037 : case TYPE_sht:
4038 0 : DO_QUANTILE_AVG(sht);
4039 : break;
4040 : case TYPE_int:
4041 0 : DO_QUANTILE_AVG(int);
4042 : break;
4043 : case TYPE_lng:
4044 0 : DO_QUANTILE_AVG(lng);
4045 : break;
4046 : #ifdef HAVE_HGE
4047 : case TYPE_hge:
4048 0 : DO_QUANTILE_AVG(hge);
4049 : break;
4050 : #endif
4051 : case TYPE_flt:
4052 0 : DO_QUANTILE_AVG(flt);
4053 : break;
4054 : case TYPE_dbl:
4055 0 : DO_QUANTILE_AVG(dbl);
4056 : break;
4057 : }
4058 : v = &val;
4059 : } else {
4060 : /* round *down* to nearest integer */
4061 38 : double f = (p - r - 1) * quantile;
4062 38 : qindex = r + p - (BUN) (p + 0.5 - f);
4063 : /* be a little paranoid about the index */
4064 38 : assert(qindex >= r && qindex < p);
4065 38 : v = BUNtail(bi, qindex);
4066 38 : if (!skip_nils && !bi.nonil)
4067 0 : nils += (*atomcmp)(v, dnil) == 0;
4068 : }
4069 53 : if (bunfastapp_nocheck(bn, v) != GDK_SUCCEED)
4070 0 : goto bunins_failed;
4071 : }
4072 17 : bat_iterator_end(&bi);
4073 17 : nils += ngrp - BATcount(bn);
4074 17 : while (BATcount(bn) < ngrp) {
4075 0 : if (bunfastapp_nocheck(bn, dnil) != GDK_SUCCEED)
4076 0 : goto bailout;
4077 : }
4078 17 : bn->theap->dirty = true;
4079 17 : BBPunfix(g->batCacheid);
4080 : } else {
4081 34 : BUN index, r, p = BATcount(b);
4082 34 : BAT *pb = NULL;
4083 34 : const oid *ords;
4084 34 : Heap *oidxh = NULL;
4085 :
4086 64 : bn = COLnew(0, average ? TYPE_dbl : tp, 1, TRANSIENT);
4087 34 : if (bn == NULL)
4088 0 : goto bailout;
4089 :
4090 34 : t1 = NULL;
4091 :
4092 34 : if (BATcheckorderidx(b)) {
4093 0 : MT_lock_set(&b->batIdxLock);
4094 0 : oidxh = b->torderidx;
4095 0 : if (oidxh != NULL)
4096 0 : HEAPincref(oidxh);
4097 0 : MT_lock_unset(&b->batIdxLock);
4098 : }
4099 4 : if (oidxh == NULL &&
4100 34 : VIEWtparent(b) &&
4101 4 : (pb = BATdescriptor(VIEWtparent(b))) != NULL) {
4102 4 : if (BATcheckorderidx(pb)) {
4103 : /* no lock on b needed since it's a view */
4104 0 : MT_lock_set(&pb->batIdxLock);
4105 0 : MT_lock_set(&pb->theaplock);
4106 0 : if (pb->tbaseoff == b->tbaseoff &&
4107 0 : BATcount(pb) == BATcount(b) &&
4108 0 : pb->hseqbase == b->hseqbase &&
4109 0 : (oidxh = pb->torderidx) != NULL) {
4110 0 : HEAPincref(oidxh);
4111 : }
4112 0 : MT_lock_unset(&pb->batIdxLock);
4113 0 : MT_lock_unset(&pb->theaplock);
4114 : }
4115 4 : BBPunfix(pb->batCacheid);
4116 : }
4117 34 : if (oidxh != NULL) {
4118 0 : MT_thread_setalgorithm(pb ? "using parent oidx" : "using oids");
4119 0 : ords = (const oid *) oidxh->base + ORDERIDXOFF;
4120 : } else {
4121 34 : if (BATsort(NULL, &t1, NULL, b, NULL, g, false, false, false) != GDK_SUCCEED)
4122 0 : goto bailout;
4123 34 : if (BATtdense(t1))
4124 : ords = NULL;
4125 : else
4126 11 : ords = (const oid *) Tloc(t1, 0);
4127 : }
4128 :
4129 34 : bi = bat_iterator(b);
4130 :
4131 34 : if (skip_nils && !bi.nonil)
4132 7 : r = binsearch(ords, 0, tp, bi.base,
4133 7 : bi.vh ? bi.vh->base : NULL,
4134 7 : bi.width, 0, p,
4135 : nil, 1, 1);
4136 : else
4137 : r = 0;
4138 :
4139 34 : if (r == p) {
4140 : /* no non-nil values, so quantile is nil */
4141 : v = dnil;
4142 : nils++;
4143 30 : } else if (average) {
4144 4 : double f = (p - r - 1) * quantile;
4145 4 : double lo = floor(f);
4146 4 : double hi = ceil(f);
4147 8 : switch (ATOMbasetype(tp)) {
4148 4 : case TYPE_bte:
4149 4 : DO_QUANTILE_AVG(bte);
4150 : break;
4151 0 : case TYPE_sht:
4152 0 : DO_QUANTILE_AVG(sht);
4153 : break;
4154 0 : case TYPE_int:
4155 0 : DO_QUANTILE_AVG(int);
4156 : break;
4157 0 : case TYPE_lng:
4158 0 : DO_QUANTILE_AVG(lng);
4159 : break;
4160 : #ifdef HAVE_HGE
4161 0 : case TYPE_hge:
4162 0 : DO_QUANTILE_AVG(hge);
4163 : break;
4164 : #endif
4165 0 : case TYPE_flt:
4166 0 : DO_QUANTILE_AVG(flt);
4167 : break;
4168 0 : case TYPE_dbl:
4169 0 : DO_QUANTILE_AVG(dbl);
4170 : break;
4171 : }
4172 : v = &val;
4173 : } else {
4174 26 : double f;
4175 : /* round (p-r-1)*quantile *down* to nearest
4176 : * integer (i.e., 1.49 and 1.5 are rounded to
4177 : * 1, 1.51 is rounded to 2) */
4178 26 : f = (p - r - 1) * quantile;
4179 26 : index = r + p - (BUN) (p + 0.5 - f);
4180 26 : if (ords)
4181 10 : index = ords[index] - b->hseqbase;
4182 : else
4183 16 : index = index + t1->tseqbase;
4184 26 : v = BUNtail(bi, index);
4185 26 : nils += (*atomcmp)(v, dnil) == 0;
4186 : }
4187 34 : if (oidxh != NULL)
4188 0 : HEAPdecref(oidxh, false);
4189 34 : BBPreclaim(t1);
4190 34 : gdk_return rc = BUNappend(bn, v, false);
4191 34 : bat_iterator_end(&bi);
4192 34 : if (rc != GDK_SUCCEED)
4193 0 : goto bailout;
4194 : }
4195 :
4196 51 : if (b != origb)
4197 17 : BBPunfix(b->batCacheid);
4198 :
4199 51 : bn->tkey = BATcount(bn) <= 1;
4200 51 : bn->tsorted = BATcount(bn) <= 1;
4201 51 : bn->trevsorted = BATcount(bn) <= 1;
4202 51 : bn->tnil = nils != 0;
4203 51 : bn->tnonil = nils == 0;
4204 51 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
4205 : "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT
4206 : ",quantile=%g,average=%s -> " ALGOOPTBATFMT
4207 : "; start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
4208 : ALGOBATPAR(origb), ALGOOPTBATPAR(origg), ALGOOPTBATPAR(e),
4209 : ALGOOPTBATPAR(s), quantile, average ? "true" : "false",
4210 : ALGOOPTBATPAR(bn), ci.seq, ci.ncand, GDKusec() - t0);
4211 : return bn;
4212 :
4213 0 : bunins_failed:
4214 0 : bat_iterator_end(&bi);
4215 0 : bailout:
4216 0 : if (b && b != origb)
4217 0 : BBPunfix(b->batCacheid);
4218 0 : if (g && g != origg)
4219 0 : BBPunfix(g->batCacheid);
4220 0 : BBPreclaim(bn);
4221 : return NULL;
4222 : }
4223 :
4224 : BAT *
4225 29 : BATgroupmedian(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
4226 : bool skip_nils)
4227 : {
4228 29 : return doBATgroupquantile(b, g, e, s, tp, 0.5,
4229 : skip_nils, false);
4230 : }
4231 :
4232 : BAT *
4233 31 : BATgroupquantile(BAT *b, BAT *g, BAT *e, BAT *s, int tp, double quantile,
4234 : bool skip_nils)
4235 : {
4236 31 : return doBATgroupquantile(b, g, e, s, tp, quantile,
4237 : skip_nils, false);
4238 : }
4239 :
4240 : BAT *
4241 3 : BATgroupmedian_avg(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
4242 : {
4243 3 : return doBATgroupquantile(b, g, e, s, tp, 0.5, skip_nils, true);
4244 : }
4245 :
4246 : BAT *
4247 3 : BATgroupquantile_avg(BAT *b, BAT *g, BAT *e, BAT *s, int tp, double quantile,
4248 : bool skip_nils)
4249 : {
4250 3 : return doBATgroupquantile(b, g, e, s, tp, quantile,
4251 : skip_nils, true);
4252 : }
4253 :
4254 : /* ---------------------------------------------------------------------- */
4255 : /* standard deviation (both biased and non-biased) */
4256 :
4257 : #define AGGR_STDEV_SINGLE(TYPE) \
4258 : do { \
4259 : TYPE x; \
4260 : TIMEOUT_LOOP_IDX(i, cnt, qry_ctx) { \
4261 : x = ((const TYPE *) values)[i]; \
4262 : if (is_##TYPE##_nil(x)) \
4263 : continue; \
4264 : n++; \
4265 : delta = (dbl) x - mean; \
4266 : mean += delta / n; \
4267 : m2 += delta * ((dbl) x - mean); \
4268 : if (isinf(m2)) \
4269 : goto overflow; \
4270 : } \
4271 : TIMEOUT_CHECK(qry_ctx, \
4272 : TIMEOUT_HANDLER(dbl_nil, qry_ctx)); \
4273 : } while (0)
4274 :
4275 : static dbl
4276 27 : calcvariance(dbl *restrict avgp, const void *restrict values, BUN cnt, int tp, bool issample)
4277 : {
4278 27 : BUN n = 0, i;
4279 27 : dbl mean = 0;
4280 27 : dbl m2 = 0;
4281 27 : dbl delta;
4282 :
4283 27 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
4284 :
4285 27 : switch (tp) {
4286 : case TYPE_bte:
4287 18 : AGGR_STDEV_SINGLE(bte);
4288 : break;
4289 : case TYPE_sht:
4290 18 : AGGR_STDEV_SINGLE(sht);
4291 : break;
4292 : case TYPE_int:
4293 48 : AGGR_STDEV_SINGLE(int);
4294 : break;
4295 : case TYPE_lng:
4296 36 : AGGR_STDEV_SINGLE(lng);
4297 : break;
4298 : #ifdef HAVE_HGE
4299 : case TYPE_hge:
4300 0 : AGGR_STDEV_SINGLE(hge);
4301 : break;
4302 : #endif
4303 : case TYPE_flt:
4304 0 : AGGR_STDEV_SINGLE(flt);
4305 : break;
4306 : case TYPE_dbl:
4307 39 : AGGR_STDEV_SINGLE(dbl);
4308 : break;
4309 0 : default:
4310 0 : GDKerror("type (%s) not supported.\n", ATOMname(tp));
4311 0 : return dbl_nil;
4312 : }
4313 25 : if (n <= (BUN) issample) {
4314 8 : if (avgp)
4315 0 : *avgp = dbl_nil;
4316 8 : return dbl_nil;
4317 : }
4318 17 : if (avgp)
4319 0 : *avgp = mean;
4320 17 : return m2 / (n - issample);
4321 2 : overflow:
4322 2 : GDKerror("22003!overflow in calculation.\n");
4323 2 : return dbl_nil;
4324 : }
4325 :
4326 : dbl
4327 13 : BATcalcstdev_population(dbl *avgp, BAT *b)
4328 : {
4329 13 : lng t0 = 0;
4330 :
4331 13 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4332 13 : BATiter bi = bat_iterator(b);
4333 13 : dbl v = calcvariance(avgp, bi.base, bi.count, bi.type, false);
4334 13 : bat_iterator_end(&bi);
4335 13 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
4336 : ALGOBATPAR(b), GDKusec() - t0);
4337 13 : return is_dbl_nil(v) ? dbl_nil : sqrt(v);
4338 : }
4339 :
4340 : dbl
4341 7 : BATcalcstdev_sample(dbl *avgp, BAT *b)
4342 : {
4343 7 : lng t0 = 0;
4344 :
4345 7 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4346 7 : BATiter bi = bat_iterator(b);
4347 7 : dbl v = calcvariance(avgp, bi.base, bi.count, bi.type, true);
4348 7 : bat_iterator_end(&bi);
4349 7 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
4350 : ALGOBATPAR(b), GDKusec() - t0);
4351 7 : return is_dbl_nil(v) ? dbl_nil : sqrt(v);
4352 : }
4353 :
4354 : dbl
4355 5 : BATcalcvariance_population(dbl *avgp, BAT *b)
4356 : {
4357 5 : lng t0 = 0;
4358 :
4359 5 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4360 5 : BATiter bi = bat_iterator(b);
4361 5 : dbl v = calcvariance(avgp, bi.base, bi.count, bi.type, false);
4362 5 : bat_iterator_end(&bi);
4363 5 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
4364 : ALGOBATPAR(b), GDKusec() - t0);
4365 5 : return v;
4366 : }
4367 :
4368 : dbl
4369 2 : BATcalcvariance_sample(dbl *avgp, BAT *b)
4370 : {
4371 2 : lng t0 = 0;
4372 :
4373 2 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4374 2 : BATiter bi = bat_iterator(b);
4375 2 : dbl v = calcvariance(avgp, bi.base, bi.count, bi.type, true);
4376 2 : bat_iterator_end(&bi);
4377 2 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
4378 : ALGOBATPAR(b), GDKusec() - t0);
4379 2 : return v;
4380 : }
4381 :
4382 : #define AGGR_COVARIANCE_SINGLE(TYPE) \
4383 : do { \
4384 : TYPE x, y; \
4385 : TIMEOUT_LOOP_IDX(i, cnt, qry_ctx) { \
4386 : x = ((const TYPE *) v1)[i]; \
4387 : y = ((const TYPE *) v2)[i]; \
4388 : if (is_##TYPE##_nil(x) || is_##TYPE##_nil(y)) \
4389 : continue; \
4390 : n++; \
4391 : delta1 = (dbl) x - mean1; \
4392 : mean1 += delta1 / n; \
4393 : delta2 = (dbl) y - mean2; \
4394 : mean2 += delta2 / n; \
4395 : m2 += delta1 * ((dbl) y - mean2); \
4396 : if (isinf(m2)) \
4397 : goto overflow; \
4398 : } \
4399 : TIMEOUT_CHECK(qry_ctx, \
4400 : TIMEOUT_HANDLER(dbl_nil, qry_ctx)); \
4401 : } while (0)
4402 :
4403 : static dbl
4404 14 : calccovariance(const void *v1, const void *v2, BUN cnt, int tp, bool issample)
4405 : {
4406 14 : BUN n = 0, i;
4407 14 : dbl mean1 = 0, mean2 = 0, m2 = 0, delta1, delta2;
4408 :
4409 14 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
4410 :
4411 :
4412 14 : switch (tp) {
4413 : case TYPE_bte:
4414 0 : AGGR_COVARIANCE_SINGLE(bte);
4415 : break;
4416 : case TYPE_sht:
4417 0 : AGGR_COVARIANCE_SINGLE(sht);
4418 : break;
4419 : case TYPE_int:
4420 111 : AGGR_COVARIANCE_SINGLE(int);
4421 : break;
4422 : case TYPE_lng:
4423 0 : AGGR_COVARIANCE_SINGLE(lng);
4424 : break;
4425 : #ifdef HAVE_HGE
4426 : case TYPE_hge:
4427 0 : AGGR_COVARIANCE_SINGLE(hge);
4428 : break;
4429 : #endif
4430 : case TYPE_flt:
4431 0 : AGGR_COVARIANCE_SINGLE(flt);
4432 : break;
4433 : case TYPE_dbl:
4434 18 : AGGR_COVARIANCE_SINGLE(dbl);
4435 : break;
4436 0 : default:
4437 0 : GDKerror("type (%s) not supported.\n", ATOMname(tp));
4438 0 : return dbl_nil;
4439 : }
4440 13 : if (n <= (BUN) issample)
4441 1 : return dbl_nil;
4442 12 : return m2 / (n - issample);
4443 1 : overflow:
4444 1 : GDKerror("22003!overflow in calculation.\n");
4445 1 : return dbl_nil;
4446 : }
4447 :
4448 : dbl
4449 9 : BATcalccovariance_population(BAT *b1, BAT *b2)
4450 : {
4451 9 : lng t0 = 0;
4452 :
4453 9 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4454 9 : BATiter b1i = bat_iterator(b1), b2i = bat_iterator(b2);
4455 9 : dbl v = calccovariance(b1i.base, b2i.base, b1i.count, b1i.type, false);
4456 9 : bat_iterator_end(&b1i);
4457 9 : bat_iterator_end(&b2i);
4458 9 : TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT " (" LLFMT " usec)\n",
4459 : ALGOBATPAR(b1), ALGOBATPAR(b2), GDKusec() - t0);
4460 9 : return v;
4461 : }
4462 :
4463 : dbl
4464 5 : BATcalccovariance_sample(BAT *b1, BAT *b2)
4465 : {
4466 5 : lng t0 = 0;
4467 :
4468 5 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4469 5 : BATiter b1i = bat_iterator(b1), b2i = bat_iterator(b2);
4470 5 : dbl v = calccovariance(b1i.base, b2i.base, b1i.count, b1i.type, true);
4471 5 : bat_iterator_end(&b1i);
4472 5 : bat_iterator_end(&b2i);
4473 5 : TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT " (" LLFMT " usec)\n",
4474 : ALGOBATPAR(b1), ALGOBATPAR(b2), GDKusec() - t0);
4475 5 : return v;
4476 : }
4477 :
4478 : #define AGGR_CORRELATION_SINGLE(TYPE) \
4479 : do { \
4480 : TYPE x, y; \
4481 : TIMEOUT_LOOP_IDX(i, cnt, qry_ctx) { \
4482 : x = ((const TYPE *) v1)[i]; \
4483 : y = ((const TYPE *) v2)[i]; \
4484 : if (is_##TYPE##_nil(x) || is_##TYPE##_nil(y)) \
4485 : continue; \
4486 : n++; \
4487 : delta1 = (dbl) x - mean1; \
4488 : mean1 += delta1 / n; \
4489 : delta2 = (dbl) y - mean2; \
4490 : mean2 += delta2 / n; \
4491 : aux = (dbl) y - mean2; \
4492 : up += delta1 * aux; \
4493 : down1 += delta1 * ((dbl) x - mean1); \
4494 : down2 += delta2 * aux; \
4495 : if (isinf(up) || isinf(down1) || isinf(down2)) \
4496 : goto overflow; \
4497 : } \
4498 : TIMEOUT_CHECK(qry_ctx, \
4499 : GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx)); \
4500 : } while (0)
4501 :
4502 : dbl
4503 19 : BATcalccorrelation(BAT *b1, BAT *b2)
4504 : {
4505 19 : BUN n = 0, i, cnt = BATcount(b1);
4506 19 : dbl mean1 = 0, mean2 = 0, up = 0, down1 = 0, down2 = 0, delta1, delta2, aux;
4507 19 : BATiter b1i = bat_iterator(b1), b2i = bat_iterator(b2);
4508 19 : const void *v1 = b1i.base, *v2 = b2i.base;
4509 19 : int tp = b1i.type;
4510 19 : lng t0 = 0;
4511 :
4512 19 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
4513 :
4514 19 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4515 :
4516 19 : switch (tp) {
4517 : case TYPE_bte:
4518 11 : AGGR_CORRELATION_SINGLE(bte);
4519 : break;
4520 : case TYPE_sht:
4521 0 : AGGR_CORRELATION_SINGLE(sht);
4522 : break;
4523 : case TYPE_int:
4524 95 : AGGR_CORRELATION_SINGLE(int);
4525 : break;
4526 : case TYPE_lng:
4527 0 : AGGR_CORRELATION_SINGLE(lng);
4528 : break;
4529 : #ifdef HAVE_HGE
4530 : case TYPE_hge:
4531 0 : AGGR_CORRELATION_SINGLE(hge);
4532 : break;
4533 : #endif
4534 : case TYPE_flt:
4535 0 : AGGR_CORRELATION_SINGLE(flt);
4536 : break;
4537 : case TYPE_dbl:
4538 17 : AGGR_CORRELATION_SINGLE(dbl);
4539 : break;
4540 0 : default:
4541 0 : GDKerror("type (%s) not supported.\n", ATOMname(tp));
4542 0 : goto bailout;
4543 : }
4544 18 : bat_iterator_end(&b1i);
4545 18 : bat_iterator_end(&b2i);
4546 18 : if (n != 0 && down1 != 0 && down2 != 0)
4547 9 : aux = (up / n) / (sqrt(down1 / n) * sqrt(down2 / n));
4548 : else
4549 9 : aux = dbl_nil;
4550 18 : TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT " (" LLFMT " usec)\n",
4551 : ALGOBATPAR(b1), ALGOBATPAR(b2), GDKusec() - t0);
4552 : return aux;
4553 1 : overflow:
4554 1 : GDKerror("22003!overflow in calculation.\n");
4555 1 : bailout:
4556 1 : bat_iterator_end(&b1i);
4557 1 : bat_iterator_end(&b2i);
4558 1 : return dbl_nil;
4559 : }
4560 :
4561 : #define AGGR_STDEV(TYPE) \
4562 : do { \
4563 : const TYPE *restrict vals = (const TYPE *) bi.base; \
4564 : TIMEOUT_LOOP(ci.ncand, qry_ctx) { \
4565 : i = canditer_next(&ci) - b->hseqbase; \
4566 : if (gids == NULL || \
4567 : (gids[i] >= min && gids[i] <= max)) { \
4568 : if (gids) \
4569 : gid = gids[i] - min; \
4570 : else \
4571 : gid = (oid) i; \
4572 : if (is_##TYPE##_nil(vals[i])) { \
4573 : if (!skip_nils) \
4574 : cnts[gid] = BUN_NONE; \
4575 : } else if (cnts[gid] != BUN_NONE) { \
4576 : cnts[gid]++; \
4577 : delta[gid] = (dbl) vals[i] - mean[gid]; \
4578 : mean[gid] += delta[gid] / cnts[gid]; \
4579 : m2[gid] += delta[gid] * ((dbl) vals[i] - mean[gid]); \
4580 : } \
4581 : } \
4582 : } \
4583 : TIMEOUT_CHECK(qry_ctx, \
4584 : GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx)); \
4585 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) { \
4586 : if (cnts[i] == 0 || cnts[i] == BUN_NONE) { \
4587 : dbls[i] = dbl_nil; \
4588 : mean[i] = dbl_nil; \
4589 : nils++; \
4590 : } else if (cnts[i] == 1) { \
4591 : dbls[i] = issample ? dbl_nil : 0; \
4592 : nils2++; \
4593 : } else if (isinf(m2[i])) { \
4594 : goto overflow; \
4595 : } else if (variance) { \
4596 : dbls[i] = m2[i] / (cnts[i] - issample); \
4597 : } else { \
4598 : dbls[i] = sqrt(m2[i] / (cnts[i] - issample)); \
4599 : } \
4600 : } \
4601 : } while (0)
4602 :
4603 : /* Calculate group standard deviation (population (i.e. biased) or
4604 : * sample (i.e. non-biased)) with optional candidates list.
4605 : *
4606 : * Note that this helper function is prepared to return two BATs: one
4607 : * (as return value) with the standard deviation per group, and one
4608 : * (as return argument) with the average per group. This isn't
4609 : * currently used since it doesn't fit into the mold of grouped
4610 : * aggregates. */
4611 : static BAT *
4612 17 : dogroupstdev(BAT **avgb, BAT *b, BAT *g, BAT *e, BAT *s, int tp,
4613 : bool skip_nils, bool issample, bool variance, const char *func)
4614 : {
4615 17 : const oid *restrict gids;
4616 17 : oid gid;
4617 17 : oid min, max;
4618 17 : BUN i, ngrp;
4619 17 : BUN nils = 0, nils2 = 0;
4620 17 : BUN *restrict cnts = NULL;
4621 17 : dbl *restrict dbls, *restrict mean, *restrict delta, *restrict m2;
4622 17 : BAT *bn = NULL, *an = NULL;
4623 17 : struct canditer ci;
4624 17 : const char *err;
4625 17 : lng t0 = 0;
4626 17 : BATiter bi = {0};
4627 :
4628 17 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
4629 :
4630 16 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4631 :
4632 16 : assert(tp == TYPE_dbl);
4633 16 : (void) tp; /* compatibility (with other BATgroup*
4634 : * functions) argument */
4635 :
4636 16 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
4637 0 : GDKerror("%s: %s\n", func, err);
4638 0 : return NULL;
4639 : }
4640 16 : if (g == NULL) {
4641 0 : GDKerror("%s: b and g must be aligned\n", func);
4642 0 : return NULL;
4643 : }
4644 :
4645 16 : if (BATcount(b) == 0 || ngrp == 0) {
4646 : /* trivial: no products, so return bat aligned with g
4647 : * with nil in the tail */
4648 2 : bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
4649 2 : goto doreturn;
4650 : }
4651 :
4652 14 : if ((e == NULL ||
4653 15 : (BATcount(e) == ci.ncand && e->hseqbase == ci.hseq)) &&
4654 6 : (BATtdense(g) || (g->tkey && g->tnonil)) &&
4655 4 : (issample || b->tnonil)) {
4656 : /* trivial: singleton groups, so all results are equal
4657 : * to zero (population) or nil (sample) */
4658 2 : dbl v = issample ? dbl_nil : 0;
4659 5 : bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &v, ngrp, TRANSIENT);
4660 6 : goto doreturn;
4661 : }
4662 :
4663 9 : delta = GDKmalloc(ngrp * sizeof(dbl));
4664 9 : m2 = GDKmalloc(ngrp * sizeof(dbl));
4665 9 : cnts = GDKzalloc(ngrp * sizeof(BUN));
4666 9 : if (avgb) {
4667 0 : an = COLnew(0, TYPE_dbl, ngrp, TRANSIENT);
4668 0 : *avgb = an;
4669 0 : if (an == NULL) {
4670 0 : mean = NULL;
4671 0 : goto alloc_fail;
4672 : }
4673 0 : mean = (dbl *) Tloc(an, 0);
4674 : } else {
4675 9 : mean = GDKmalloc(ngrp * sizeof(dbl));
4676 : }
4677 9 : if (mean == NULL || delta == NULL || m2 == NULL || cnts == NULL)
4678 0 : goto alloc_fail;
4679 :
4680 9 : bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
4681 9 : if (bn == NULL)
4682 0 : goto alloc_fail;
4683 9 : dbls = (dbl *) Tloc(bn, 0);
4684 :
4685 1080116 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) {
4686 1080034 : mean[i] = 0;
4687 1080034 : delta[i] = 0;
4688 1080034 : m2[i] = 0;
4689 : }
4690 :
4691 9 : bi = bat_iterator(b);
4692 9 : if (BATtdense(g))
4693 : gids = NULL;
4694 : else
4695 8 : gids = (const oid *) Tloc(g, 0);
4696 :
4697 9 : switch (b->ttype) {
4698 0 : case TYPE_bte:
4699 0 : AGGR_STDEV(bte);
4700 : break;
4701 0 : case TYPE_sht:
4702 0 : AGGR_STDEV(sht);
4703 : break;
4704 4 : case TYPE_int:
4705 5760376 : AGGR_STDEV(int);
4706 : break;
4707 0 : case TYPE_lng:
4708 0 : AGGR_STDEV(lng);
4709 : break;
4710 : #ifdef HAVE_HGE
4711 0 : case TYPE_hge:
4712 0 : AGGR_STDEV(hge);
4713 : break;
4714 : #endif
4715 0 : case TYPE_flt:
4716 0 : AGGR_STDEV(flt);
4717 : break;
4718 5 : case TYPE_dbl:
4719 88 : AGGR_STDEV(dbl);
4720 : break;
4721 0 : default:
4722 0 : GDKerror("%s: type (%s) not supported.\n",
4723 : func, ATOMname(b->ttype));
4724 0 : goto bailout;
4725 : }
4726 7 : bat_iterator_end(&bi);
4727 7 : if (an) {
4728 0 : BATsetcount(an, ngrp);
4729 0 : an->tkey = ngrp <= 1;
4730 0 : an->tsorted = ngrp <= 1;
4731 0 : an->trevsorted = ngrp <= 1;
4732 0 : an->tnil = nils != 0;
4733 0 : an->tnonil = nils == 0;
4734 : } else {
4735 7 : GDKfree(mean);
4736 : }
4737 7 : if (issample)
4738 3 : nils += nils2;
4739 7 : GDKfree(delta);
4740 7 : GDKfree(m2);
4741 7 : GDKfree(cnts);
4742 7 : BATsetcount(bn, ngrp);
4743 7 : bn->tkey = ngrp <= 1;
4744 7 : bn->tsorted = ngrp <= 1;
4745 7 : bn->trevsorted = ngrp <= 1;
4746 7 : bn->tnil = nils != 0;
4747 7 : bn->tnonil = nils == 0;
4748 15 : doreturn:
4749 15 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOBATFMT ",e=" ALGOOPTBATFMT
4750 : ",s=" ALGOOPTBATFMT
4751 : ",skip_nils=%s,issample=%s,variance=%s -> " ALGOOPTBATFMT
4752 : ",avgb=" ALGOOPTBATFMT " (%s -- " LLFMT " usec)\n",
4753 : ALGOBATPAR(b), ALGOBATPAR(g), ALGOOPTBATPAR(e),
4754 : ALGOOPTBATPAR(s),
4755 : skip_nils ? "true" : "false",
4756 : issample ? "true" : "false",
4757 : variance ? "true" : "false",
4758 : ALGOOPTBATPAR(bn), ALGOOPTBATPAR(an),
4759 : func, GDKusec() - t0);
4760 : return bn;
4761 2 : overflow:
4762 2 : GDKerror("22003!overflow in calculation.\n");
4763 2 : bailout:
4764 2 : bat_iterator_end(&bi);
4765 2 : alloc_fail:
4766 2 : if (an)
4767 0 : BBPreclaim(an);
4768 : else
4769 2 : GDKfree(mean);
4770 2 : BBPreclaim(bn);
4771 2 : GDKfree(delta);
4772 2 : GDKfree(m2);
4773 2 : GDKfree(cnts);
4774 2 : return NULL;
4775 : }
4776 :
4777 : BAT *
4778 7 : BATgroupstdev_sample(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
4779 : {
4780 7 : return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, true, false,
4781 : __func__);
4782 : }
4783 :
4784 : BAT *
4785 5 : BATgroupstdev_population(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
4786 : bool skip_nils)
4787 : {
4788 5 : return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, false, false,
4789 : __func__);
4790 : }
4791 :
4792 : BAT *
4793 1 : BATgroupvariance_sample(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
4794 : bool skip_nils)
4795 : {
4796 1 : return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, true, true,
4797 : __func__);
4798 : }
4799 :
4800 : BAT *
4801 4 : BATgroupvariance_population(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
4802 : bool skip_nils)
4803 : {
4804 4 : return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, false, true,
4805 : __func__);
4806 : }
4807 :
4808 : #define AGGR_COVARIANCE(TYPE) \
4809 : do { \
4810 : const TYPE *vals1 = (const TYPE *) b1i.base; \
4811 : const TYPE *vals2 = (const TYPE *) b2i.base; \
4812 : TIMEOUT_LOOP(ci.ncand, qry_ctx) { \
4813 : i = canditer_next(&ci) - b1->hseqbase; \
4814 : if (gids == NULL || \
4815 : (gids[i] >= min && gids[i] <= max)) { \
4816 : if (gids) \
4817 : gid = gids[i] - min; \
4818 : else \
4819 : gid = (oid) i; \
4820 : if (is_##TYPE##_nil(vals1[i]) || is_##TYPE##_nil(vals2[i])) { \
4821 : if (!skip_nils) \
4822 : cnts[gid] = BUN_NONE; \
4823 : } else if (cnts[gid] != BUN_NONE) { \
4824 : cnts[gid]++; \
4825 : delta1[gid] = (dbl) vals1[i] - mean1[gid]; \
4826 : mean1[gid] += delta1[gid] / cnts[gid]; \
4827 : delta2[gid] = (dbl) vals2[i] - mean2[gid]; \
4828 : mean2[gid] += delta2[gid] / cnts[gid]; \
4829 : m2[gid] += delta1[gid] * ((dbl) vals2[i] - mean2[gid]); \
4830 : } \
4831 : } \
4832 : } \
4833 : TIMEOUT_CHECK(qry_ctx, \
4834 : GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx)); \
4835 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) { \
4836 : if (cnts[i] == 0 || cnts[i] == BUN_NONE) { \
4837 : dbls[i] = dbl_nil; \
4838 : nils++; \
4839 : } else if (cnts[i] == 1) { \
4840 : dbls[i] = issample ? dbl_nil : 0; \
4841 : nils2++; \
4842 : } else if (isinf(m2[i])) { \
4843 : goto overflow; \
4844 : } else { \
4845 : dbls[i] = m2[i] / (cnts[i] - issample); \
4846 : } \
4847 : } \
4848 : } while (0)
4849 :
4850 : static BAT *
4851 60 : dogroupcovariance(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp,
4852 : bool skip_nils, bool issample, const char *func)
4853 : {
4854 60 : const oid *restrict gids;
4855 60 : oid gid, min, max;
4856 60 : BUN i, ngrp, nils = 0, nils2 = 0;
4857 60 : BUN *restrict cnts = NULL;
4858 60 : dbl *restrict dbls, *restrict mean1, *restrict mean2, *restrict delta1, *restrict delta2, *restrict m2;
4859 60 : BAT *bn = NULL;
4860 60 : struct canditer ci;
4861 60 : const char *err;
4862 60 : lng t0 = 0;
4863 60 : BATiter b1i = {0}, b2i = {0};
4864 :
4865 60 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
4866 :
4867 :
4868 60 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4869 :
4870 60 : assert(tp == TYPE_dbl && BATcount(b1) == BATcount(b2) && b1->ttype == b2->ttype && BATtdense(b1) == BATtdense(b2));
4871 60 : (void) tp;
4872 :
4873 60 : if ((err = BATgroupaggrinit(b1, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
4874 0 : GDKerror("%s: %s\n", func, err);
4875 0 : return NULL;
4876 : }
4877 60 : if (g == NULL) {
4878 0 : GDKerror("%s: b1, b2 and g must be aligned\n", func);
4879 0 : return NULL;
4880 : }
4881 :
4882 60 : if (BATcount(b1) == 0 || ngrp == 0) {
4883 0 : bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
4884 0 : goto doreturn;
4885 : }
4886 :
4887 60 : if ((e == NULL ||
4888 60 : (BATcount(e) == BATcount(b1) && (e->hseqbase == b1->hseqbase || e->hseqbase == b2->hseqbase))) &&
4889 20 : (BATtdense(g) || (g->tkey && g->tnonil)) &&
4890 10 : (issample || (b1->tnonil && b2->tnonil))) {
4891 : /* trivial: singleton groups, so all results are equal
4892 : * to zero (population) or nil (sample) */
4893 10 : dbl v = issample ? dbl_nil : 0;
4894 13 : bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &v, ngrp, TRANSIENT);
4895 12 : goto doreturn;
4896 : }
4897 :
4898 47 : delta1 = GDKmalloc(ngrp * sizeof(dbl));
4899 47 : delta2 = GDKmalloc(ngrp * sizeof(dbl));
4900 47 : m2 = GDKmalloc(ngrp * sizeof(dbl));
4901 47 : cnts = GDKzalloc(ngrp * sizeof(BUN));
4902 47 : mean1 = GDKmalloc(ngrp * sizeof(dbl));
4903 47 : mean2 = GDKmalloc(ngrp * sizeof(dbl));
4904 :
4905 47 : if (mean1 == NULL || mean2 == NULL || delta1 == NULL || delta2 == NULL || m2 == NULL || cnts == NULL)
4906 0 : goto alloc_fail;
4907 :
4908 404 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) {
4909 310 : m2[i] = 0;
4910 310 : mean1[i] = 0;
4911 310 : mean2[i] = 0;
4912 : }
4913 :
4914 47 : bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
4915 47 : if (bn == NULL)
4916 0 : goto alloc_fail;
4917 47 : dbls = (dbl *) Tloc(bn, 0);
4918 :
4919 47 : if (!g || BATtdense(g))
4920 : gids = NULL;
4921 : else
4922 47 : gids = (const oid *) Tloc(g, 0);
4923 :
4924 47 : b1i = bat_iterator(b1);
4925 45 : b2i = bat_iterator(b2);
4926 46 : switch (b1->ttype) {
4927 9 : case TYPE_bte:
4928 184 : AGGR_COVARIANCE(bte);
4929 : break;
4930 0 : case TYPE_sht:
4931 0 : AGGR_COVARIANCE(sht);
4932 : break;
4933 37 : case TYPE_int:
4934 745 : AGGR_COVARIANCE(int);
4935 : break;
4936 0 : case TYPE_lng:
4937 0 : AGGR_COVARIANCE(lng);
4938 : break;
4939 : #ifdef HAVE_HGE
4940 0 : case TYPE_hge:
4941 0 : AGGR_COVARIANCE(hge);
4942 : break;
4943 : #endif
4944 0 : case TYPE_flt:
4945 0 : AGGR_COVARIANCE(flt);
4946 : break;
4947 0 : case TYPE_dbl:
4948 0 : AGGR_COVARIANCE(dbl);
4949 : break;
4950 0 : default:
4951 0 : GDKerror("%s: type (%s) not supported.\n", func, ATOMname(b1->ttype));
4952 0 : goto bailout;
4953 : }
4954 47 : bat_iterator_end(&b1i);
4955 44 : bat_iterator_end(&b2i);
4956 47 : GDKfree(mean1);
4957 47 : GDKfree(mean2);
4958 :
4959 46 : if (issample)
4960 19 : nils += nils2;
4961 46 : GDKfree(delta1);
4962 47 : GDKfree(delta2);
4963 47 : GDKfree(m2);
4964 47 : GDKfree(cnts);
4965 47 : BATsetcount(bn, ngrp);
4966 46 : bn->tkey = ngrp <= 1;
4967 46 : bn->tsorted = ngrp <= 1;
4968 46 : bn->trevsorted = ngrp <= 1;
4969 46 : bn->tnil = nils != 0;
4970 46 : bn->tnonil = nils == 0;
4971 58 : doreturn:
4972 58 : TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT ",g=" ALGOBATFMT
4973 : ",e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT
4974 : ",skip_nils=%s,issample=%s -> " ALGOOPTBATFMT
4975 : " (%s -- " LLFMT " usec)\n",
4976 : ALGOBATPAR(b1), ALGOBATPAR(b2), ALGOBATPAR(g),
4977 : ALGOOPTBATPAR(e), ALGOOPTBATPAR(s),
4978 : skip_nils ? "true" : "false",
4979 : issample ? "true" : "false",
4980 : ALGOOPTBATPAR(bn),
4981 : func, GDKusec() - t0);
4982 : return bn;
4983 0 : overflow:
4984 0 : GDKerror("22003!overflow in calculation.\n");
4985 0 : bailout:
4986 0 : bat_iterator_end(&b1i);
4987 0 : bat_iterator_end(&b2i);
4988 0 : alloc_fail:
4989 0 : BBPreclaim(bn);
4990 0 : GDKfree(mean1);
4991 0 : GDKfree(mean2);
4992 0 : GDKfree(delta1);
4993 0 : GDKfree(delta2);
4994 0 : GDKfree(m2);
4995 0 : GDKfree(cnts);
4996 0 : return NULL;
4997 : }
4998 :
4999 : BAT *
5000 30 : BATgroupcovariance_sample(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
5001 : {
5002 30 : return dogroupcovariance(b1, b2, g, e, s, tp, skip_nils, true,
5003 : __func__);
5004 : }
5005 :
5006 : BAT *
5007 30 : BATgroupcovariance_population(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
5008 : {
5009 30 : return dogroupcovariance(b1, b2, g, e, s, tp, skip_nils, false,
5010 : __func__);
5011 : }
5012 :
5013 : #define AGGR_CORRELATION(TYPE) \
5014 : do { \
5015 : const TYPE *vals1 = (const TYPE *) b1i.base; \
5016 : const TYPE *vals2 = (const TYPE *) b2i.base; \
5017 : TIMEOUT_LOOP(ci.ncand, qry_ctx) { \
5018 : i = canditer_next(&ci) - b1->hseqbase; \
5019 : if (gids == NULL || \
5020 : (gids[i] >= min && gids[i] <= max)) { \
5021 : if (gids) \
5022 : gid = gids[i] - min; \
5023 : else \
5024 : gid = (oid) i; \
5025 : if (is_##TYPE##_nil(vals1[i]) || is_##TYPE##_nil(vals2[i])) { \
5026 : if (!skip_nils) \
5027 : cnts[gid] = BUN_NONE; \
5028 : } else if (cnts[gid] != BUN_NONE) { \
5029 : cnts[gid]++; \
5030 : delta1[gid] = (dbl) vals1[i] - mean1[gid]; \
5031 : mean1[gid] += delta1[gid] / cnts[gid]; \
5032 : delta2[gid] = (dbl) vals2[i] - mean2[gid]; \
5033 : mean2[gid] += delta2[gid] / cnts[gid]; \
5034 : aux = (dbl) vals2[i] - mean2[gid]; \
5035 : up[gid] += delta1[gid] * aux; \
5036 : down1[gid] += delta1[gid] * ((dbl) vals1[i] - mean1[gid]); \
5037 : down2[gid] += delta2[gid] * aux; \
5038 : } \
5039 : } \
5040 : } \
5041 : TIMEOUT_CHECK(qry_ctx, \
5042 : GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx)); \
5043 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) { \
5044 : if (cnts[i] <= 1 || cnts[i] == BUN_NONE || down1[i] == 0 || down2[i] == 0) { \
5045 : dbls[i] = dbl_nil; \
5046 : nils++; \
5047 : } else if (isinf(up[i]) || isinf(down1[i]) || isinf(down2[i])) { \
5048 : goto overflow; \
5049 : } else { \
5050 : dbls[i] = (up[i] / cnts[i]) / (sqrt(down1[i] / cnts[i]) * sqrt(down2[i] / cnts[i])); \
5051 : assert(!is_dbl_nil(dbls[i])); \
5052 : } \
5053 : } \
5054 : } while (0)
5055 :
5056 : BAT *
5057 49 : BATgroupcorrelation(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
5058 : {
5059 49 : const oid *restrict gids;
5060 49 : oid gid, min, max;
5061 49 : BUN i, ngrp, nils = 0;
5062 49 : BUN *restrict cnts = NULL;
5063 49 : dbl *restrict dbls, *restrict mean1, *restrict mean2, *restrict delta1, *restrict delta2, *restrict up, *restrict down1, *restrict down2, aux;
5064 49 : BAT *bn = NULL;
5065 49 : struct canditer ci;
5066 49 : const char *err;
5067 49 : lng t0 = 0;
5068 49 : BATiter b1i = {0}, b2i = {0};
5069 :
5070 49 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
5071 :
5072 49 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
5073 :
5074 49 : assert(tp == TYPE_dbl && BATcount(b1) == BATcount(b2) && b1->ttype == b2->ttype && BATtdense(b1) == BATtdense(b2));
5075 49 : (void) tp;
5076 :
5077 49 : if ((err = BATgroupaggrinit(b1, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
5078 0 : GDKerror("%s\n", err);
5079 0 : return NULL;
5080 : }
5081 48 : if (g == NULL) {
5082 0 : GDKerror("b1, b2 and g must be aligned\n");
5083 0 : return NULL;
5084 : }
5085 :
5086 48 : if (BATcount(b1) == 0 || ngrp == 0) {
5087 1 : bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
5088 1 : goto doreturn;
5089 : }
5090 :
5091 47 : if ((e == NULL ||
5092 48 : (BATcount(e) == BATcount(b1) && (e->hseqbase == b1->hseqbase || e->hseqbase == b2->hseqbase))) &&
5093 13 : (BATtdense(g) || (g->tkey && g->tnonil))) {
5094 13 : dbl v = dbl_nil;
5095 13 : bn = BATconstant(min, TYPE_dbl, &v, ngrp, TRANSIENT);
5096 13 : goto doreturn;
5097 : }
5098 :
5099 34 : delta1 = GDKmalloc(ngrp * sizeof(dbl));
5100 34 : delta2 = GDKmalloc(ngrp * sizeof(dbl));
5101 34 : up = GDKmalloc(ngrp * sizeof(dbl));
5102 34 : down1 = GDKmalloc(ngrp * sizeof(dbl));
5103 34 : down2 = GDKmalloc(ngrp * sizeof(dbl));
5104 34 : cnts = GDKzalloc(ngrp * sizeof(BUN));
5105 32 : mean1 = GDKmalloc(ngrp * sizeof(dbl));
5106 34 : mean2 = GDKmalloc(ngrp * sizeof(dbl));
5107 :
5108 34 : if (mean1 == NULL || mean2 == NULL || delta1 == NULL || delta2 == NULL || up == NULL || down1 == NULL || down2 == NULL || cnts == NULL)
5109 0 : goto alloc_fail;
5110 :
5111 273 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) {
5112 205 : up[i] = 0;
5113 205 : down1[i] = 0;
5114 205 : down2[i] = 0;
5115 205 : mean1[i] = 0;
5116 205 : mean2[i] = 0;
5117 : }
5118 :
5119 34 : bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
5120 32 : if (bn == NULL)
5121 0 : goto alloc_fail;
5122 32 : dbls = (dbl *) Tloc(bn, 0);
5123 :
5124 32 : if (!g || BATtdense(g))
5125 : gids = NULL;
5126 : else
5127 32 : gids = (const oid *) Tloc(g, 0);
5128 :
5129 32 : b1i = bat_iterator(b1);
5130 34 : b2i = bat_iterator(b2);
5131 31 : switch (b1->ttype) {
5132 4 : case TYPE_bte:
5133 80 : AGGR_CORRELATION(bte);
5134 : break;
5135 0 : case TYPE_sht:
5136 0 : AGGR_CORRELATION(sht);
5137 : break;
5138 26 : case TYPE_int:
5139 1218 : AGGR_CORRELATION(int);
5140 : break;
5141 0 : case TYPE_lng:
5142 0 : AGGR_CORRELATION(lng);
5143 : break;
5144 : #ifdef HAVE_HGE
5145 0 : case TYPE_hge:
5146 0 : AGGR_CORRELATION(hge);
5147 : break;
5148 : #endif
5149 0 : case TYPE_flt:
5150 0 : AGGR_CORRELATION(flt);
5151 : break;
5152 1 : case TYPE_dbl:
5153 4 : AGGR_CORRELATION(dbl);
5154 : break;
5155 0 : default:
5156 0 : GDKerror("type (%s) not supported.\n", ATOMname(b1->ttype));
5157 0 : goto bailout;
5158 : }
5159 33 : bat_iterator_end(&b1i);
5160 31 : bat_iterator_end(&b2i);
5161 33 : GDKfree(mean1);
5162 33 : GDKfree(mean2);
5163 33 : GDKfree(delta1);
5164 32 : GDKfree(delta2);
5165 33 : GDKfree(up);
5166 33 : GDKfree(down1);
5167 32 : GDKfree(down2);
5168 33 : GDKfree(cnts);
5169 33 : BATsetcount(bn, ngrp);
5170 32 : bn->tkey = ngrp <= 1;
5171 32 : bn->tsorted = ngrp <= 1;
5172 32 : bn->trevsorted = ngrp <= 1;
5173 32 : bn->tnil = nils != 0;
5174 32 : bn->tnonil = nils == 0;
5175 46 : doreturn:
5176 46 : TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT ",g=" ALGOBATFMT
5177 : ",e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT
5178 : ",skip_nils=%s -> " ALGOOPTBATFMT
5179 : " (" LLFMT " usec)\n",
5180 : ALGOBATPAR(b1), ALGOBATPAR(b2), ALGOBATPAR(g),
5181 : ALGOOPTBATPAR(e), ALGOOPTBATPAR(s),
5182 : skip_nils ? "true" : "false",
5183 : ALGOOPTBATPAR(bn),
5184 : GDKusec() - t0);
5185 : return bn;
5186 1 : overflow:
5187 1 : GDKerror("22003!overflow in calculation.\n");
5188 1 : bailout:
5189 1 : bat_iterator_end(&b1i);
5190 1 : bat_iterator_end(&b2i);
5191 1 : alloc_fail:
5192 1 : BBPreclaim(bn);
5193 1 : GDKfree(mean1);
5194 1 : GDKfree(mean2);
5195 1 : GDKfree(delta1);
5196 1 : GDKfree(delta2);
5197 1 : GDKfree(up);
5198 1 : GDKfree(down1);
5199 1 : GDKfree(down2);
5200 1 : GDKfree(cnts);
5201 1 : return NULL;
5202 : }
|
