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 43881 : BATgroupaggrinit(BAT *b, BAT *g, BAT *e, BAT *s,
66 : /* outputs: */
67 : oid *minp, oid *maxp, BUN *ngrpp,
68 : struct canditer *ci)
69 : {
70 43881 : oid min, max;
71 43881 : BUN i, ngrp;
72 43881 : const oid *restrict gids;
73 :
74 43881 : if (b == NULL)
75 : return "b must exist";
76 43881 : canditer_init(ci, b, s);
77 43907 : if (g) {
78 34875 : if (ci->ncand != BATcount(g) ||
79 19314 : (ci->ncand != 0 && ci->seq != g->hseqbase))
80 : return "b with s and g must be aligned";
81 34858 : assert(BATttype(g) == TYPE_oid);
82 : }
83 : if (g == NULL) {
84 : min = 0;
85 : max = 0;
86 : ngrp = 1;
87 34858 : } 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 34855 : ngrp = BATcount(e);
138 34855 : min = e->hseqbase;
139 34855 : max = e->hseqbase + ngrp - 1;
140 : }
141 43890 : *minp = min;
142 43890 : *maxp = max;
143 43890 : *ngrpp = ngrp;
144 :
145 43890 : 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 31674 : twosum(volatile double *hi, volatile double *lo, double x, double y)
163 : {
164 31674 : volatile double yr;
165 :
166 31674 : assert(fabs(x) >= fabs(y));
167 :
168 31674 : *hi = x + y;
169 31674 : yr = *hi - x;
170 31674 : *lo = y - yr;
171 31674 : }
172 :
173 : static inline void
174 13769 : exchange(double *x, double *y)
175 : {
176 13769 : double t = *x;
177 13769 : *x = *y;
178 13769 : *y = t;
179 13769 : }
180 :
181 : /* this function was adapted from https://bugs.python.org/file10357/msum4.py */
182 : BUN
183 1097 : 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 1097 : 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 1097 : BUN listi;
201 1097 : int parti;
202 1097 : int i;
203 1097 : BUN grp;
204 1097 : double x, y;
205 1097 : volatile double lo, hi;
206 1097 : double twopow = pow((double) FLT_RADIX, (double) (DBL_MAX_EXP - 1));
207 1097 : BUN nils = 0;
208 1097 : volatile flt f;
209 :
210 1097 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
211 :
212 : /* we only deal with the two floating point types */
213 1097 : assert(tp1 == TYPE_flt || tp1 == TYPE_dbl);
214 1097 : assert(tp2 == TYPE_flt || tp2 == TYPE_dbl);
215 : /* if input is dbl, then so it output */
216 1097 : assert(tp2 == TYPE_flt || tp1 == TYPE_dbl);
217 : /* if no gids, then we have a single group */
218 1097 : assert(ngrp == 1 || gids != NULL);
219 1097 : if (gids == NULL || ngrp == 1) {
220 1080 : min = max = 0;
221 1080 : ngrp = 1;
222 1080 : gids = NULL;
223 : }
224 1097 : pergroup = GDKmalloc(ngrp * sizeof(*pergroup));
225 1097 : if (pergroup == NULL)
226 : return BUN_NONE;
227 10725 : for (grp = 0; grp < ngrp; grp++) {
228 19256 : pergroup[grp] = (struct pergroup) {
229 : .maxpartials = 2,
230 9628 : .partials = GDKmalloc(2 * sizeof(double)),
231 : };
232 9628 : 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 28150 : TIMEOUT_LOOP(ci->ncand, qry_ctx) {
240 25956 : listi = canditer_next(ci) - seqb;
241 25956 : grp = gids ? gids[listi] : 0;
242 25956 : if (grp < min || grp > max)
243 0 : continue;
244 25956 : if (pergroup[grp].partials == NULL)
245 0 : continue;
246 25956 : if (tp1 == TYPE_flt && !is_flt_nil(((const flt *) values)[listi]))
247 16 : x = ((const flt *) values)[listi];
248 25940 : 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 25869 : pergroup[grp].valseen = true;
265 : #ifdef INFINITES_ALLOWED
266 : if (isinf(x)) {
267 : pergroup[grp].infs += x;
268 : continue;
269 : }
270 : #endif
271 25869 : i = 0;
272 56002 : for (parti = 0; parti < pergroup[grp].npartials; parti++) {
273 30133 : y = pergroup[grp].partials[parti];
274 30133 : if (fabs(x) < fabs(y))
275 13761 : exchange(&x, &y);
276 30133 : twosum(&hi, &lo, x, y);
277 30133 : 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 30133 : if (lo != 0)
287 16428 : pergroup[grp].partials[i++] = lo;
288 30133 : x = hi;
289 : }
290 25869 : if (x != 0) {
291 25798 : 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 25798 : pergroup[grp].partials[i++] = x;
300 : }
301 25869 : pergroup[grp].npartials = i;
302 : }
303 1097 : TIMEOUT_CHECK(qry_ctx, GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx));
304 10697 : for (grp = 0; grp < ngrp; grp++) {
305 9628 : if (pergroup[grp].partials == NULL)
306 0 : continue;
307 9628 : 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 9597 : 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 9591 : if (pergroup[grp].infs != 0)
373 28 : goto overflow;
374 :
375 9563 : 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 9549 : hi = pergroup[grp].partials[--pergroup[grp].npartials];
387 9549 : 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 9549 : 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 9549 : GDKfree(pergroup[grp].partials);
403 9549 : pergroup[grp].partials = NULL;
404 9549 : 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 9538 : } else if (is_dbl_nil(hi)) {
412 0 : goto overflow;
413 : } else {
414 9538 : ((dbl *) results)[grp] = hi;
415 : }
416 : }
417 1069 : GDKfree(pergroup);
418 1069 : 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 : ADD_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 : ADD_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 : ADD_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 : ADD_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 : ADD_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 11144 : 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 11144 : BUN nils = 0;
716 11144 : BUN i;
717 11144 : oid gid;
718 11144 : unsigned int *restrict seen = NULL; /* bitmask for groups that we've seen */
719 :
720 11144 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
721 :
722 11150 : 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 1080 : if (tp1 != TYPE_flt && tp1 != TYPE_dbl)
729 0 : goto unsupported;
730 1092 : *algo = "sum: floating point";
731 1092 : 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 10058 : seen = GDKzalloc(((ngrp + 31) / 32) * sizeof(int));
738 10054 : if (seen == NULL) {
739 : return BUN_NONE;
740 : }
741 :
742 10054 : switch (tp2) {
743 6 : case TYPE_bte: {
744 6 : bte *restrict sums = (bte *) results;
745 6 : switch (tp1) {
746 6 : case TYPE_bte:
747 27 : AGGR_SUM(bte, bte);
748 : break;
749 0 : default:
750 0 : goto unsupported;
751 : }
752 : break;
753 : }
754 21 : case TYPE_sht: {
755 21 : sht *restrict sums = (sht *) results;
756 21 : switch (tp1) {
757 13 : case TYPE_bte:
758 13 : if (ci->ncand < ((BUN) 1 << ((sizeof(sht) - sizeof(bte)) << 3)))
759 179 : AGGR_SUM_NOOVL(bte, sht);
760 : else
761 0 : AGGR_SUM(bte, sht);
762 : break;
763 8 : case TYPE_sht:
764 47 : AGGR_SUM(sht, sht);
765 : break;
766 0 : default:
767 0 : goto unsupported;
768 : }
769 : break;
770 : }
771 915 : case TYPE_int: {
772 915 : int *restrict sums = (int *) results;
773 915 : 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 915 : case TYPE_int:
787 5908 : AGGR_SUM(int, int);
788 : break;
789 0 : default:
790 0 : goto unsupported;
791 : }
792 : break;
793 : }
794 7500 : case TYPE_lng: {
795 7500 : lng *restrict sums = (lng *) results;
796 7500 : 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 113 : case TYPE_int:
821 113 : if (ci->ncand < ((BUN) 1 << ((sizeof(lng) - sizeof(int)) << 3)))
822 466338 : AGGR_SUM_NOOVL(int, lng);
823 : else
824 0 : AGGR_SUM(int, lng);
825 : break;
826 : #endif
827 7387 : case TYPE_lng:
828 8547890 : AGGR_SUM(lng, lng);
829 : break;
830 0 : default:
831 0 : goto unsupported;
832 : }
833 : break;
834 : }
835 : #ifdef HAVE_HGE
836 1612 : case TYPE_hge: {
837 1612 : hge *sums = (hge *) results;
838 1612 : switch (tp1) {
839 165 : case TYPE_bte:
840 4035576 : AGGR_SUM_NOOVL(bte, hge);
841 : break;
842 8 : case TYPE_sht:
843 178 : AGGR_SUM_NOOVL(sht, hge);
844 : break;
845 817 : case TYPE_int:
846 53182756 : AGGR_SUM_NOOVL(int, hge);
847 : break;
848 53 : case TYPE_lng:
849 4774834 : AGGR_SUM_NOOVL(lng, hge);
850 : break;
851 569 : case TYPE_hge:
852 35454436 : 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 10064 : if (nils == 0 && nil_if_empty) {
865 : /* figure out whether there were any empty groups
866 : * (that result in a nil value) */
867 10055 : if (ngrp & 0x1F) {
868 : /* fill last slot */
869 10044 : seen[ngrp >> 5] |= ~0U << (ngrp & 0x1F);
870 : }
871 231817 : for (i = 0, ngrp = (ngrp + 31) / 32; i < ngrp; i++) {
872 222244 : if (seen[i] != ~0U) {
873 : nils = 1;
874 : break;
875 : }
876 : }
877 : }
878 10064 : GDKfree(seen);
879 :
880 10064 : 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 6241 : BATgroupsum(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
901 : {
902 6241 : const oid *restrict gids;
903 6241 : oid min, max;
904 6241 : BUN ngrp;
905 6241 : BUN nils;
906 6241 : BAT *bn;
907 6241 : struct canditer ci;
908 6241 : const char *err;
909 6241 : const char *algo = NULL;
910 6241 : lng t0 = 0;
911 :
912 6241 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
913 :
914 6241 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
915 0 : GDKerror("%s\n", err);
916 0 : return NULL;
917 : }
918 6238 : if (g == NULL) {
919 0 : GDKerror("b and g must be aligned\n");
920 0 : return NULL;
921 : }
922 :
923 6238 : if (ci.ncand == 0 || ngrp == 0) {
924 : /* trivial: no sums, so return bat aligned with g with
925 : * nil in the tail */
926 1986 : return BATconstant(ngrp == 0 ? 0 : min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
927 : }
928 :
929 4252 : if ((e == NULL ||
930 4251 : (BATcount(e) == ci.ncand && e->hseqbase == ci.hseq)) &&
931 1390 : (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 1389 : return BATconvert(b, s, tp, 0, 0, 0);
935 : }
936 :
937 2863 : bn = BATconstant(min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
938 2856 : if (bn == NULL) {
939 : return NULL;
940 : }
941 :
942 2856 : if (BATtdense(g))
943 : gids = NULL;
944 : else
945 2536 : gids = (const oid *) Tloc(g, 0);
946 :
947 2856 : BATiter bi = bat_iterator(b);
948 5729 : nils = dosum(bi.base, bi.nonil, b->hseqbase, &ci,
949 2863 : Tloc(bn, 0), ngrp, bi.type, tp, gids, min, max,
950 : skip_nils, true, __func__, &algo);
951 2866 : bat_iterator_end(&bi);
952 :
953 2857 : if (nils < BUN_NONE) {
954 2857 : BATsetcount(bn, ngrp);
955 2865 : bn->tkey = BATcount(bn) <= 1;
956 2865 : bn->tsorted = BATcount(bn) <= 1;
957 2865 : bn->trevsorted = BATcount(bn) <= 1;
958 2865 : bn->tnil = nils != 0;
959 2865 : bn->tnonil = nils == 0;
960 : } else {
961 0 : BBPunfix(bn->batCacheid);
962 0 : bn = NULL;
963 : }
964 :
965 2866 : if (algo)
966 2866 : MT_thread_setalgorithm(algo);
967 2860 : 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 8718 : BATsum(void *res, int tp, BAT *b, BAT *s, bool skip_nils, bool nil_if_empty)
1019 : {
1020 8718 : oid min, max;
1021 8718 : BUN ngrp;
1022 8718 : struct canditer ci;
1023 8718 : const char *err;
1024 8718 : const char *algo = NULL;
1025 8718 : lng t0 = 0;
1026 :
1027 8718 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
1028 :
1029 8718 : 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 8710 : 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 8680 : switch (tp) {
1085 9 : case TYPE_bte:
1086 9 : * (bte *) res = nil_if_empty ? bte_nil : 0;
1087 9 : break;
1088 16 : case TYPE_sht:
1089 16 : * (sht *) res = nil_if_empty ? sht_nil : 0;
1090 16 : break;
1091 505 : case TYPE_int:
1092 505 : * (int *) res = nil_if_empty ? int_nil : 0;
1093 505 : break;
1094 6406 : case TYPE_lng:
1095 6406 : * (lng *) res = nil_if_empty ? lng_nil : 0;
1096 6406 : break;
1097 : #ifdef HAVE_HGE
1098 556 : case TYPE_hge:
1099 556 : * (hge *) res = nil_if_empty ? hge_nil : 0;
1100 556 : break;
1101 : #endif
1102 1188 : case TYPE_flt:
1103 : case TYPE_dbl:
1104 1188 : 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 1187 : break;
1158 : }
1159 1187 : if (b->ttype == TYPE_dbl)
1160 1174 : * (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 8679 : if (ci.ncand == 0)
1170 : return GDK_SUCCEED;
1171 8286 : BATiter bi = bat_iterator(b);
1172 16611 : BUN nils = dosum(bi.base, bi.nonil, b->hseqbase, &ci,
1173 8306 : res, true, bi.type, tp, &min, min, max,
1174 : skip_nils, nil_if_empty, __func__, &algo);
1175 8305 : bat_iterator_end(&bi);
1176 8293 : if (algo)
1177 8293 : MT_thread_setalgorithm(algo);
1178 8291 : 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 8291 : 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 : MUL4_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 42 : 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 42 : BUN nils = 0;
1355 42 : BUN i;
1356 42 : oid gid;
1357 42 : unsigned int *restrict seen; /* bitmask for groups that we've seen */
1358 :
1359 42 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
1360 :
1361 : /* allocate bitmap for seen group ids */
1362 42 : seen = GDKzalloc(((ngrp + 31) / 32) * sizeof(int));
1363 42 : if (seen == NULL) {
1364 : return BUN_NONE;
1365 : }
1366 :
1367 42 : switch (tp2) {
1368 0 : case TYPE_bte: {
1369 0 : bte *restrict prods = (bte *) results;
1370 0 : switch (tp1) {
1371 0 : case TYPE_bte:
1372 0 : AGGR_PROD(bte, bte, sht);
1373 : break;
1374 0 : default:
1375 0 : goto unsupported;
1376 : }
1377 : break;
1378 : }
1379 0 : case TYPE_sht: {
1380 0 : sht *restrict prods = (sht *) results;
1381 0 : switch (tp1) {
1382 0 : case TYPE_bte:
1383 0 : AGGR_PROD(bte, sht, int);
1384 : break;
1385 0 : case TYPE_sht:
1386 0 : AGGR_PROD(sht, sht, int);
1387 : break;
1388 0 : default:
1389 0 : goto unsupported;
1390 : }
1391 : break;
1392 : }
1393 0 : case TYPE_int: {
1394 0 : int *restrict prods = (int *) results;
1395 0 : switch (tp1) {
1396 0 : case TYPE_bte:
1397 0 : AGGR_PROD(bte, int, lng);
1398 : break;
1399 0 : case TYPE_sht:
1400 0 : AGGR_PROD(sht, int, lng);
1401 : break;
1402 0 : case TYPE_int:
1403 0 : AGGR_PROD(int, int, lng);
1404 : break;
1405 0 : default:
1406 0 : goto unsupported;
1407 : }
1408 : break;
1409 : }
1410 : #ifdef HAVE_HGE
1411 0 : case TYPE_lng: {
1412 0 : lng *prods = (lng *) results;
1413 0 : switch (tp1) {
1414 0 : case TYPE_bte:
1415 0 : AGGR_PROD(bte, lng, hge);
1416 : break;
1417 0 : case TYPE_sht:
1418 0 : AGGR_PROD(sht, lng, hge);
1419 : break;
1420 0 : case TYPE_int:
1421 0 : AGGR_PROD(int, lng, hge);
1422 : break;
1423 0 : case TYPE_lng:
1424 0 : AGGR_PROD(lng, lng, hge);
1425 : break;
1426 0 : default:
1427 0 : goto unsupported;
1428 : }
1429 : break;
1430 : }
1431 35 : case TYPE_hge: {
1432 35 : hge *prods = (hge *) results;
1433 35 : switch (tp1) {
1434 6 : case TYPE_bte:
1435 18 : AGGR_PROD_HGE(bte);
1436 : break;
1437 6 : case TYPE_sht:
1438 18 : AGGR_PROD_HGE(sht);
1439 : break;
1440 6 : case TYPE_int:
1441 18 : AGGR_PROD_HGE(int);
1442 : break;
1443 12 : case TYPE_lng:
1444 36 : AGGR_PROD_HGE(lng);
1445 : break;
1446 5 : case TYPE_hge:
1447 40 : AGGR_PROD_HGE(hge);
1448 : break;
1449 0 : default:
1450 0 : goto unsupported;
1451 : }
1452 : break;
1453 : }
1454 : #else
1455 : case TYPE_lng: {
1456 : lng *restrict prods = (lng *) results;
1457 : switch (tp1) {
1458 : case TYPE_bte:
1459 : AGGR_PROD_LNG(bte);
1460 : break;
1461 : case TYPE_sht:
1462 : AGGR_PROD_LNG(sht);
1463 : break;
1464 : case TYPE_int:
1465 : AGGR_PROD_LNG(int);
1466 : break;
1467 : case TYPE_lng:
1468 : AGGR_PROD_LNG(lng);
1469 : break;
1470 : default:
1471 : goto unsupported;
1472 : }
1473 : break;
1474 : }
1475 : #endif
1476 0 : case TYPE_flt: {
1477 0 : flt *restrict prods = (flt *) results;
1478 0 : switch (tp1) {
1479 0 : case TYPE_bte:
1480 0 : AGGR_PROD_FLOAT(bte, flt);
1481 : break;
1482 0 : case TYPE_sht:
1483 0 : AGGR_PROD_FLOAT(sht, flt);
1484 : break;
1485 0 : case TYPE_int:
1486 0 : AGGR_PROD_FLOAT(int, flt);
1487 : break;
1488 0 : case TYPE_lng:
1489 0 : AGGR_PROD_FLOAT(lng, flt);
1490 : break;
1491 : #ifdef HAVE_HGE
1492 0 : case TYPE_hge:
1493 0 : AGGR_PROD_FLOAT(hge, flt);
1494 : break;
1495 : #endif
1496 0 : case TYPE_flt:
1497 0 : AGGR_PROD_FLOAT(flt, flt);
1498 : break;
1499 0 : default:
1500 0 : goto unsupported;
1501 : }
1502 : break;
1503 : }
1504 7 : case TYPE_dbl: {
1505 7 : dbl *restrict prods = (dbl *) results;
1506 7 : 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 7 : case TYPE_dbl:
1528 35 : 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 42 : if (nils == 0 && nil_if_empty) {
1540 : /* figure out whether there were any empty groups
1541 : * (that result in a nil value) */
1542 42 : if (ngrp & 0x1F) {
1543 : /* fill last slot */
1544 42 : seen[ngrp >> 5] |= ~0U << (ngrp & 0x1F);
1545 : }
1546 84 : for (i = 0, ngrp = (ngrp + 31) / 32; i < ngrp; i++) {
1547 42 : if (seen[i] != ~0U) {
1548 : nils = 1;
1549 : break;
1550 : }
1551 : }
1552 : }
1553 42 : GDKfree(seen);
1554 :
1555 42 : 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 21 : BATgroupprod(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
1576 : {
1577 21 : const oid *restrict gids;
1578 21 : oid min, max;
1579 21 : BUN ngrp;
1580 21 : BUN nils;
1581 21 : BAT *bn;
1582 21 : struct canditer ci;
1583 21 : const char *err;
1584 21 : lng t0 = 0;
1585 :
1586 21 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
1587 :
1588 21 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
1589 0 : GDKerror("%s\n", err);
1590 0 : return NULL;
1591 : }
1592 21 : if (g == NULL) {
1593 0 : GDKerror("b and g must be aligned\n");
1594 0 : return NULL;
1595 : }
1596 :
1597 21 : if (ci.ncand == 0 || ngrp == 0) {
1598 : /* trivial: no products, so return bat aligned with g
1599 : * with nil in the tail */
1600 10 : return BATconstant(ngrp == 0 ? 0 : min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
1601 : }
1602 :
1603 11 : if ((e == NULL ||
1604 11 : (BATcount(e) == ci.ncand && e->hseqbase == ci.hseq)) &&
1605 9 : (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 9 : return BATconvert(b, s, tp, 0, 0, 0);
1609 : }
1610 :
1611 2 : bn = BATconstant(min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
1612 2 : if (bn == NULL) {
1613 : return NULL;
1614 : }
1615 :
1616 2 : if (BATtdense(g))
1617 : gids = NULL;
1618 : else
1619 2 : gids = (const oid *) Tloc(g, 0);
1620 :
1621 2 : BATiter bi = bat_iterator(b);
1622 4 : nils = doprod(bi.base, b->hseqbase, &ci, Tloc(bn, 0), ngrp,
1623 2 : bi.type, tp, gids, true, min, max, skip_nils,
1624 : true, __func__);
1625 2 : bat_iterator_end(&bi);
1626 :
1627 2 : if (nils < BUN_NONE) {
1628 2 : BATsetcount(bn, ngrp);
1629 2 : bn->tkey = BATcount(bn) <= 1;
1630 2 : bn->tsorted = BATcount(bn) <= 1;
1631 2 : bn->trevsorted = BATcount(bn) <= 1;
1632 2 : bn->tnil = nils != 0;
1633 2 : bn->tnonil = nils == 0;
1634 : } else {
1635 0 : BBPunfix(bn->batCacheid);
1636 0 : bn = NULL;
1637 : }
1638 :
1639 2 : 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 40 : BATprod(void *res, int tp, BAT *b, BAT *s, bool skip_nils, bool nil_if_empty)
1651 : {
1652 40 : oid min, max;
1653 40 : BUN ngrp;
1654 40 : BUN nils;
1655 40 : struct canditer ci;
1656 40 : const char *err;
1657 40 : lng t0 = 0;
1658 :
1659 40 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
1660 :
1661 40 : 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 40 : 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 5 : case TYPE_dbl:
1687 5 : * (dbl *) res = nil_if_empty ? dbl_nil : (dbl) 1;
1688 5 : 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 40 : if (ci.ncand == 0)
1695 : return GDK_SUCCEED;
1696 40 : BATiter bi = bat_iterator(b);
1697 80 : nils = doprod(bi.base, b->hseqbase, &ci, res, true,
1698 40 : bi.type, tp, &min, false, min, max,
1699 : skip_nils, nil_if_empty, __func__);
1700 40 : bat_iterator_end(&bi);
1701 40 : 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 40 : return nils < BUN_NONE ? GDK_SUCCEED : GDK_FAIL;
1706 : }
1707 :
1708 : /* ---------------------------------------------------------------------- */
1709 : /* average */
1710 :
1711 : #define GOTO_BAILOUT() \
1712 : do { \
1713 : GDKfree(avgs); \
1714 : GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx); \
1715 : } while (0)
1716 :
1717 : #define AGGR_AVG(TYPE) \
1718 : do { \
1719 : const TYPE *restrict vals = (const TYPE *) bi.base; \
1720 : TYPE *restrict avgs = GDKzalloc(ngrp * sizeof(TYPE)); \
1721 : if (avgs == NULL) \
1722 : goto bailout; \
1723 : TIMEOUT_LOOP(ci.ncand, qry_ctx) { \
1724 : i = canditer_next(&ci) - b->hseqbase; \
1725 : if (gids == NULL || \
1726 : (gids[i] >= min && gids[i] <= max)) { \
1727 : if (gids) \
1728 : gid = gids[i] - min; \
1729 : else \
1730 : gid = (oid) i; \
1731 : if (is_##TYPE##_nil(vals[i])) { \
1732 : if (!skip_nils) \
1733 : cnts[gid] = lng_nil; \
1734 : } else if (!is_lng_nil(cnts[gid])) { \
1735 : AVERAGE_ITER(TYPE, vals[i], \
1736 : avgs[gid], \
1737 : rems[gid], \
1738 : cnts[gid]); \
1739 : } \
1740 : } \
1741 : } \
1742 : TIMEOUT_CHECK(qry_ctx, GOTO_BAILOUT()); \
1743 : for (i = 0; i < ngrp; i++) { \
1744 : if (cnts[i] == 0 || is_lng_nil(cnts[i])) { \
1745 : dbls[i] = dbl_nil; \
1746 : cnts[i] = 0; \
1747 : nils++; \
1748 : } else { \
1749 : dbls[i] = avgs[i] + (dbl) rems[i] / cnts[i]; \
1750 : } \
1751 : } \
1752 : GDKfree(avgs); \
1753 : } while (0)
1754 :
1755 : #define AGGR_AVG_FLOAT(TYPE) \
1756 : do { \
1757 : const TYPE *restrict vals = (const TYPE *) bi.base; \
1758 : for (i = 0; i < ngrp; i++) \
1759 : dbls[i] = 0; \
1760 : TIMEOUT_LOOP(ci.ncand, qry_ctx) { \
1761 : i = canditer_next(&ci) - b->hseqbase; \
1762 : if (gids == NULL || \
1763 : (gids[i] >= min && gids[i] <= max)) { \
1764 : if (gids) \
1765 : gid = gids[i] - min; \
1766 : else \
1767 : gid = (oid) i; \
1768 : if (is_##TYPE##_nil(vals[i])) { \
1769 : if (!skip_nils) \
1770 : cnts[gid] = lng_nil; \
1771 : } else if (!is_lng_nil(cnts[gid])) { \
1772 : AVERAGE_ITER_FLOAT(TYPE, vals[i], \
1773 : dbls[gid], \
1774 : cnts[gid]); \
1775 : } \
1776 : } \
1777 : } \
1778 : TIMEOUT_CHECK(qry_ctx, \
1779 : GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx)); \
1780 : for (i = 0; i < ngrp; i++) { \
1781 : if (cnts[i] == 0 || is_lng_nil(cnts[i])) { \
1782 : dbls[i] = dbl_nil; \
1783 : cnts[i] = 0; \
1784 : nils++; \
1785 : } \
1786 : } \
1787 : } while (0)
1788 :
1789 : /* There are three functions that are used for calculating averages.
1790 : * The first one (BATgroupavg) returns averages as a floating point
1791 : * value, the other two (BATgroupavg3 and BATgroupavg3combine) work
1792 : * together to return averages in the domain type (which should be an
1793 : * integer type). */
1794 :
1795 : /* Calculate group averages with optional candidates list. The average
1796 : * that is calculated is returned in a dbl, independent of the type of
1797 : * the input. The average is calculated exactly, so not in a floating
1798 : * point which could potentially losse bits during processing
1799 : * (e.g. average of 2**62 and a billion 1's). */
1800 : gdk_return
1801 319 : BATgroupavg(BAT **bnp, BAT **cntsp, BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, int scale)
1802 : {
1803 319 : const oid *restrict gids;
1804 319 : oid gid;
1805 319 : oid min, max;
1806 319 : BUN i, ngrp;
1807 319 : BUN nils = 0;
1808 319 : lng *restrict rems = NULL;
1809 319 : lng *restrict cnts = NULL;
1810 319 : dbl *restrict dbls;
1811 319 : BAT *bn = NULL, *cn = NULL;
1812 319 : struct canditer ci;
1813 319 : const char *err;
1814 319 : lng t0 = 0;
1815 319 : BATiter bi = {0};
1816 :
1817 319 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
1818 :
1819 319 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
1820 :
1821 319 : assert(tp == TYPE_dbl);
1822 319 : (void) tp; /* compatibility (with other BATgroup*
1823 : * functions) argument */
1824 :
1825 319 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
1826 0 : GDKerror("%s\n", err);
1827 0 : return GDK_FAIL;
1828 : }
1829 317 : if (g == NULL) {
1830 0 : GDKerror("b and g must be aligned\n");
1831 0 : return GDK_FAIL;
1832 : }
1833 :
1834 317 : 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 308 : if ((!skip_nils || cntsp == NULL || b->tnonil) &&
1854 195 : (e == NULL ||
1855 195 : (BATcount(e) == ci.ncand && e->hseqbase == b->hseqbase)) &&
1856 115 : (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 115 : if ((bn = BATconvert(b, s, TYPE_dbl, 0, 0, 0)) == NULL)
1860 : return GDK_FAIL;
1861 116 : if (cntsp) {
1862 15 : lng one = 1;
1863 15 : if ((cn = BATconstant(ngrp == 0 ? 0 : min, TYPE_lng, &one, ngrp, TRANSIENT)) == NULL) {
1864 0 : BBPreclaim(bn);
1865 0 : return GDK_FAIL;
1866 : }
1867 15 : *cntsp = cn;
1868 : }
1869 116 : *bnp = bn;
1870 116 : return GDK_SUCCEED;
1871 : }
1872 :
1873 : /* allocate temporary space to do per group calculations */
1874 193 : switch (b->ttype) {
1875 167 : 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 167 : rems = GDKzalloc(ngrp * sizeof(lng));
1883 168 : if (rems == NULL)
1884 0 : goto bailout1;
1885 : break;
1886 : default:
1887 : break;
1888 : }
1889 194 : if (cntsp) {
1890 172 : if ((cn = COLnew(min, TYPE_lng, ngrp, TRANSIENT)) == NULL)
1891 0 : goto bailout1;
1892 172 : cnts = (lng *) Tloc(cn, 0);
1893 172 : memset(cnts, 0, ngrp * sizeof(lng));
1894 : } else {
1895 22 : cnts = GDKzalloc(ngrp * sizeof(lng));
1896 22 : if (cnts == NULL)
1897 0 : goto bailout1;
1898 : }
1899 :
1900 194 : bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
1901 194 : if (bn == NULL)
1902 0 : goto bailout1;
1903 194 : dbls = (dbl *) Tloc(bn, 0);
1904 :
1905 194 : if (BATtdense(g))
1906 : gids = NULL;
1907 : else
1908 109 : gids = (const oid *) Tloc(g, 0);
1909 :
1910 194 : bi = bat_iterator(b);
1911 194 : switch (b->ttype) {
1912 0 : case TYPE_bte:
1913 0 : AGGR_AVG(bte);
1914 0 : break;
1915 5 : case TYPE_sht:
1916 25 : AGGR_AVG(sht);
1917 5 : break;
1918 142 : case TYPE_int:
1919 21487123 : AGGR_AVG(int);
1920 142 : break;
1921 21 : case TYPE_lng:
1922 111 : AGGR_AVG(lng);
1923 21 : 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 328 : 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 193 : bat_iterator_end(&bi);
1940 192 : GDKfree(rems);
1941 192 : if (cn == NULL)
1942 22 : GDKfree(cnts);
1943 : else {
1944 170 : BATsetcount(cn, ngrp);
1945 170 : cn->tkey = BATcount(cn) <= 1;
1946 170 : cn->tsorted = BATcount(cn) <= 1;
1947 170 : cn->trevsorted = BATcount(cn) <= 1;
1948 170 : cn->tnil = false;
1949 170 : cn->tnonil = true;
1950 170 : *cntsp = cn;
1951 : }
1952 192 : 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 192 : BATsetcount(bn, ngrp);
1960 192 : bn->tkey = BATcount(bn) <= 1;
1961 192 : bn->tsorted = BATcount(bn) <= 1;
1962 192 : bn->trevsorted = BATcount(bn) <= 1;
1963 192 : bn->tnil = nils != 0;
1964 192 : bn->tnonil = nils == 0;
1965 192 : *bnp = bn;
1966 192 : 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 411 : BATgroupavg3(BAT **avgp, BAT **remp, BAT **cntp, BAT *b, BAT *g, BAT *e, BAT *s, bool skip_nils)
1997 : {
1998 411 : const char *err;
1999 411 : oid min, max;
2000 411 : BUN ngrp;
2001 411 : struct canditer ci;
2002 411 : BAT *bn, *rn, *cn;
2003 411 : BUN i;
2004 411 : oid o;
2005 :
2006 411 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
2007 :
2008 411 : 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 411 : 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 370 : ValRecord zero;
2031 370 : (void) VALinit(&zero, TYPE_bte, &(bte){0});
2032 370 : bn = BATconstant(min, b->ttype, VALconvert(b->ttype, &zero),
2033 : ngrp, TRANSIENT);
2034 371 : rn = BATconstant(min, TYPE_lng, &(lng){0}, ngrp, TRANSIENT);
2035 374 : cn = BATconstant(min, TYPE_lng, &(lng){0}, ngrp, TRANSIENT);
2036 372 : 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 372 : lng *rems = Tloc(rn, 0);
2043 372 : lng *cnts = Tloc(cn, 0);
2044 372 : const oid *gids = g && !BATtdense(g) ? Tloc(g, 0) : NULL;
2045 372 : oid gid = ngrp == 1 && gids ? gids[0] - min : 0;
2046 :
2047 372 : BATiter bi = bat_iterator(b);
2048 :
2049 732 : 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 284 : TIMEOUT_LOOP(ci.ncand, qry_ctx) {
2054 252 : o = canditer_next(&ci) - b->hseqbase;
2055 252 : if (ngrp > 1)
2056 0 : gid = gids ? gids[o] - min : o;
2057 252 : 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 252 : } else if (!is_lng_nil(cnts[gid])) {
2067 252 : AVERAGE_ITER(bte, vals[o], avgs[gid], rems[gid], cnts[gid]);
2068 : }
2069 : }
2070 48 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) {
2071 16 : if (cnts[i] == 0) {
2072 0 : avgs[i] = bte_nil;
2073 0 : bn->tnil = true;
2074 : } else
2075 : #ifdef TRUNCATE_NUMBERS
2076 : if (!is_bte_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0) {
2077 : avgs[i]++;
2078 : rems[i] -= cnts[i];
2079 : }
2080 : #else
2081 16 : if (!is_bte_nil(avgs[i]) && rems[i] > 0) {
2082 15 : if (avgs[i] < 0) {
2083 0 : if (2*rems[i] > cnts[i]) {
2084 0 : avgs[i]++;
2085 0 : rems[i] -= cnts[i];
2086 : }
2087 : } else {
2088 15 : if (2*rems[i] >= cnts[i]) {
2089 8 : avgs[i]++;
2090 8 : rems[i] -= cnts[i];
2091 : }
2092 : }
2093 : }
2094 : #endif
2095 : }
2096 : break;
2097 : }
2098 1 : case TYPE_sht: {
2099 1 : const sht *vals = (const sht *) bi.base;
2100 1 : sht *avgs = Tloc(bn, 0);
2101 21 : TIMEOUT_LOOP(ci.ncand, qry_ctx) {
2102 19 : o = canditer_next(&ci) - b->hseqbase;
2103 19 : if (ngrp > 1)
2104 0 : gid = gids ? gids[o] - min : o;
2105 19 : if (is_sht_nil(vals[o])) {
2106 0 : if (!skip_nils) {
2107 0 : avgs[gid] = sht_nil;
2108 0 : rems[gid] = lng_nil;
2109 0 : cnts[gid] = lng_nil;
2110 0 : bn->tnil = true;
2111 0 : rn->tnil = true;
2112 0 : cn->tnil = true;
2113 : }
2114 19 : } else if (!is_lng_nil(cnts[gid])) {
2115 19 : AVERAGE_ITER(sht, vals[o], avgs[gid], rems[gid], cnts[gid]);
2116 : }
2117 : }
2118 3 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) {
2119 1 : if (cnts[i] == 0) {
2120 0 : avgs[i] = sht_nil;
2121 0 : bn->tnil = true;
2122 : } else
2123 : #ifdef TRUNCATE_NUMBERS
2124 : if (!is_sht_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0) {
2125 : avgs[i]++;
2126 : rems[i] -= cnts[i];
2127 : }
2128 : #else
2129 1 : if (!is_sht_nil(avgs[i]) && rems[i] > 0) {
2130 1 : if (avgs[i] < 0) {
2131 0 : if (2*rems[i] > cnts[i]) {
2132 0 : avgs[i]++;
2133 0 : rems[i] -= cnts[i];
2134 : }
2135 : } else {
2136 1 : if (2*rems[i] >= cnts[i]) {
2137 0 : avgs[i]++;
2138 0 : rems[i] -= cnts[i];
2139 : }
2140 : }
2141 : }
2142 : #endif
2143 : }
2144 : break;
2145 : }
2146 239 : case TYPE_int: {
2147 239 : const int *vals = (const int *) bi.base;
2148 239 : int *avgs = Tloc(bn, 0);
2149 1826195 : TIMEOUT_LOOP(ci.ncand, qry_ctx) {
2150 1825369 : o = canditer_next(&ci) - b->hseqbase;
2151 1825372 : if (ngrp > 1)
2152 197459 : gid = gids ? gids[o] - min : o;
2153 1825372 : if (is_int_nil(vals[o])) {
2154 122721 : 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 1702651 : } else if (!is_lng_nil(cnts[gid])) {
2163 1844437 : AVERAGE_ITER(int, vals[o], avgs[gid], rems[gid], cnts[gid]);
2164 : }
2165 : }
2166 211071 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) {
2167 210579 : if (cnts[i] == 0) {
2168 1032 : avgs[i] = int_nil;
2169 1032 : 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 209547 : if (!is_int_nil(avgs[i]) && rems[i] > 0) {
2178 77645 : if (avgs[i] < 0) {
2179 59061 : if (2*rems[i] > cnts[i]) {
2180 22039 : avgs[i]++;
2181 22039 : rems[i] -= cnts[i];
2182 : }
2183 : } else {
2184 18584 : if (2*rems[i] >= cnts[i]) {
2185 14860 : avgs[i]++;
2186 14860 : rems[i] -= cnts[i];
2187 : }
2188 : }
2189 : }
2190 : #endif
2191 : }
2192 : break;
2193 : }
2194 100 : case TYPE_lng: {
2195 100 : const lng *vals = (const lng *) bi.base;
2196 100 : lng *avgs = Tloc(bn, 0);
2197 143653 : TIMEOUT_LOOP(ci.ncand, qry_ctx) {
2198 143454 : o = canditer_next(&ci) - b->hseqbase;
2199 143453 : if (ngrp > 1)
2200 115941 : gid = gids ? gids[o] - min : o;
2201 143453 : 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 143244 : } else if (!is_lng_nil(cnts[gid])) {
2211 143983 : AVERAGE_ITER(lng, vals[o], avgs[gid], rems[gid], cnts[gid]);
2212 : }
2213 : }
2214 59301 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) {
2215 59098 : if (cnts[i] == 0) {
2216 158 : avgs[i] = lng_nil;
2217 158 : 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 58940 : if (!is_lng_nil(avgs[i]) && rems[i] > 0) {
2226 844 : 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 715 : if (2*rems[i] >= cnts[i]) {
2233 682 : avgs[i]++;
2234 682 : 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 371 : bat_iterator_end(&bi);
2294 370 : TIMEOUT_CHECK(qry_ctx, GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx));
2295 371 : BATsetcount(bn, ngrp);
2296 373 : BATsetcount(rn, ngrp);
2297 374 : BATsetcount(cn, ngrp);
2298 374 : bn->tnonil = !bn->tnil;
2299 374 : rn->tnonil = !rn->tnil;
2300 374 : cn->tnonil = !cn->tnil;
2301 374 : bn->tkey = rn->tkey = cn->tkey = ngrp == 1;
2302 374 : bn->tsorted = rn->tsorted = cn->tsorted = ngrp == 1;
2303 374 : bn->trevsorted = rn->trevsorted = cn->trevsorted = ngrp == 1;
2304 374 : *avgp = bn;
2305 374 : *remp = rn;
2306 374 : *cntp = cn;
2307 374 : 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 141753 : 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 141753 : lng cnt = cnt1 + cnt2;
2454 :
2455 141753 : if (rem2 < 0) {
2456 28945 : avg2--;
2457 28945 : rem2 += cnt2;
2458 : }
2459 141753 : *cntp = cnt;
2460 : #ifdef HAVE_HGE
2461 141753 : hge v = avg1 * cnt1 + rem1 + avg2 * cnt2 + rem2;
2462 141753 : int a = (int) (v / cnt);
2463 141753 : v %= cnt;
2464 141753 : if (v < 0) {
2465 104947 : a--;
2466 104947 : v += cnt;
2467 : }
2468 141753 : *avgp = a;
2469 141753 : *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 141753 : }
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 39 : BATgroupavg3combine(BAT *avg, BAT *rem, BAT *cnt, BAT *g, BAT *e, bool skip_nils)
2635 : {
2636 39 : const char *err;
2637 39 : oid min, max;
2638 39 : BUN ngrp;
2639 39 : struct canditer ci;
2640 39 : BUN i;
2641 39 : BAT *bn, *rn, *cn;
2642 :
2643 39 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
2644 :
2645 39 : if ((err = BATgroupaggrinit(avg, g, e, NULL, &min, &max, &ngrp, &ci)) != NULL) {
2646 0 : GDKerror("%s\n", err);
2647 0 : return NULL;
2648 : }
2649 39 : assert(ci.tpe == cand_dense);
2650 39 : if (ci.ncand != BATcount(rem) || ci.ncand != BATcount(cnt)) {
2651 0 : GDKerror("input bats not aligned");
2652 0 : return NULL;
2653 : }
2654 39 : 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 38 : ValRecord zero;
2659 38 : (void) VALinit(&zero, TYPE_bte, &(bte){0});
2660 38 : bn = BATconstant(min, avg->ttype, VALconvert(avg->ttype, &zero),
2661 : ngrp, TRANSIENT);
2662 : /* rn and cn are temporary storage of intermediates */
2663 39 : rn = BATconstant(min, TYPE_lng, &(lng){0}, ngrp, TRANSIENT);
2664 38 : cn = BATconstant(min, TYPE_lng, &(lng){0}, ngrp, TRANSIENT);
2665 39 : 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 39 : lng *rems = Tloc(rn, 0);
2672 39 : lng *cnts = Tloc(cn, 0);
2673 39 : const lng *orems = Tloc(rem, 0);
2674 39 : const lng *ocnts = Tloc(cnt, 0);
2675 39 : const oid *gids = g && !BATtdense(g) ? Tloc(g, 0) : NULL;
2676 39 : oid gid = ngrp == 1 && gids ? gids[0] - min : 0;
2677 :
2678 39 : BATiter bi = bat_iterator(avg);
2679 :
2680 76 : 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 31 : case TYPE_int: {
2764 31 : const int *vals = (const int *) bi.base;
2765 31 : int *avgs = Tloc(bn, 0);
2766 142560 : TIMEOUT_LOOP_IDX(i, ci.ncand, qry_ctx) {
2767 142458 : if (ngrp > 1)
2768 142299 : gid = gids ? gids[i] - min : i;
2769 142458 : if (is_int_nil(vals[i])) {
2770 705 : if (!skip_nils) {
2771 0 : avgs[gid] = int_nil;
2772 0 : bn->tnil = true;
2773 : }
2774 141753 : } else if (!is_int_nil(avgs[gid])) {
2775 141753 : combine_averages_int(&avgs[gid], &rems[gid],
2776 : &cnts[gid], avgs[gid],
2777 141753 : rems[gid], cnts[gid],
2778 141753 : vals[i], orems[i],
2779 141753 : ocnts[i]);
2780 : }
2781 : }
2782 38024 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) {
2783 37959 : if (cnts[i] == 0) {
2784 42 : avgs[i] = int_nil;
2785 42 : 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 37917 : if (!is_int_nil(avgs[i]) && rems[i] > 0) {
2792 20009 : if (avgs[i] < 0) {
2793 16717 : if (2*rems[i] > cnts[i])
2794 8094 : avgs[i]++;
2795 : } else {
2796 3292 : if (2*rems[i] >= cnts[i])
2797 2947 : 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 39 : bat_iterator_end(&bi);
2890 38 : BBPreclaim(rn);
2891 39 : BBPreclaim(cn);
2892 39 : TIMEOUT_CHECK(qry_ctx, GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx));
2893 39 : BATsetcount(bn, ngrp);
2894 39 : bn->tnonil = !bn->tnil;
2895 39 : bn->tkey = ngrp == 1;
2896 39 : bn->tsorted = ngrp == 1;
2897 39 : bn->trevsorted = ngrp == 1;
2898 39 : 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 : ADD_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 7850 : BATcalcavg(BAT *b, BAT *s, dbl *avg, BUN *vals, int scale)
2988 : {
2989 7850 : lng n = 0, r = 0;
2990 7850 : BUN i = 0;
2991 : #ifdef HAVE_HGE
2992 7850 : hge sum = 0;
2993 : #else
2994 : lng sum = 0;
2995 : #endif
2996 7850 : struct canditer ci;
2997 7850 : const void *restrict src;
2998 :
2999 7850 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
3000 :
3001 7850 : canditer_init(&ci, b, s);
3002 :
3003 7853 : BATiter bi = bat_iterator(b);
3004 7854 : src = bi.base;
3005 :
3006 7854 : 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 7790 : case TYPE_int:
3014 4749760 : 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 161 : 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 7849 : bat_iterator_end(&bi);
3036 7851 : if (scale != 0 && !is_dbl_nil(*avg))
3037 0 : *avg /= pow(10.0, (double) scale);
3038 7851 : if (vals)
3039 7851 : *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 26490 : BATgroupcount(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
3070 : {
3071 26490 : const oid *restrict gids;
3072 26490 : oid gid;
3073 26490 : oid min, max;
3074 26490 : BUN i, ngrp;
3075 26490 : lng *restrict cnts;
3076 26490 : BAT *bn = NULL;
3077 26490 : int t;
3078 26490 : const void *nil;
3079 26490 : int (*atomcmp)(const void *, const void *);
3080 26490 : struct canditer ci;
3081 26490 : const char *err;
3082 26490 : lng t0 = 0;
3083 26490 : BATiter bi = {0};
3084 :
3085 26490 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
3086 :
3087 26494 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
3088 :
3089 26494 : assert(tp == TYPE_lng);
3090 26494 : (void) tp; /* compatibility (with other BATgroup* */
3091 :
3092 26494 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
3093 0 : GDKerror("%s\n", err);
3094 0 : return NULL;
3095 : }
3096 26489 : if (g == NULL) {
3097 0 : GDKerror("b and g must be aligned\n");
3098 0 : return NULL;
3099 : }
3100 :
3101 26489 : if (ci.ncand == 0 || ngrp == 0) {
3102 : /* trivial: no counts, so return bat aligned with g
3103 : * with zero in the tail */
3104 13334 : lng zero = 0;
3105 13334 : return BATconstant(ngrp == 0 ? 0 : min, TYPE_lng, &zero, ngrp, TRANSIENT);
3106 : }
3107 :
3108 13155 : bn = COLnew(min, TYPE_lng, ngrp, TRANSIENT);
3109 13148 : if (bn == NULL)
3110 : return NULL;
3111 13148 : cnts = (lng *) Tloc(bn, 0);
3112 13148 : memset(cnts, 0, ngrp * sizeof(lng));
3113 :
3114 13148 : if (BATtdense(g))
3115 : gids = NULL;
3116 : else
3117 10476 : gids = (const oid *) Tloc(g, 0);
3118 :
3119 13148 : 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 12879 : if (gids) {
3123 17694755 : TIMEOUT_LOOP(ci.ncand, qry_ctx) {
3124 17672836 : i = canditer_next(&ci) - b->hseqbase;
3125 17672860 : if (gids[i] >= min && gids[i] <= max)
3126 17672601 : cnts[gids[i] - min]++;
3127 : }
3128 : } else {
3129 72441 : TIMEOUT_LOOP(ci.ncand, qry_ctx) {
3130 67543 : i = canditer_next(&ci) - b->hseqbase;
3131 67544 : cnts[i] = 1;
3132 : }
3133 : }
3134 : } else {
3135 269 : t = b->ttype;
3136 269 : nil = ATOMnilptr(t);
3137 269 : atomcmp = ATOMcompare(t);
3138 269 : t = ATOMbasetype(t);
3139 :
3140 269 : bi = bat_iterator(b);
3141 :
3142 269 : switch (t) {
3143 4 : case TYPE_bte:
3144 14 : AGGR_COUNT(bte);
3145 : break;
3146 0 : case TYPE_sht:
3147 0 : AGGR_COUNT(sht);
3148 : break;
3149 250 : case TYPE_int:
3150 16245 : AGGR_COUNT(int);
3151 : break;
3152 4 : case TYPE_lng:
3153 12 : 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 6 : default:
3167 154 : TIMEOUT_LOOP(ci.ncand, qry_ctx) {
3168 141 : i = canditer_next(&ci) - b->hseqbase;
3169 144 : if (gids == NULL ||
3170 144 : (gids[i] >= min && gids[i] <= max)) {
3171 145 : if (gids)
3172 145 : gid = gids[i] - min;
3173 : else
3174 : gid = (oid) i;
3175 145 : if ((*atomcmp)(BUNtail(bi, i), nil) != 0) {
3176 61 : cnts[gid]++;
3177 : }
3178 : }
3179 : }
3180 : break;
3181 : }
3182 272 : bat_iterator_end(&bi);
3183 : }
3184 13157 : TIMEOUT_CHECK(qry_ctx, GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx));
3185 13147 : BATsetcount(bn, ngrp);
3186 13150 : bn->tkey = BATcount(bn) <= 1;
3187 13150 : bn->tsorted = BATcount(bn) <= 1;
3188 13150 : bn->trevsorted = BATcount(bn) <= 1;
3189 13150 : bn->tnil = false;
3190 13150 : bn->tnonil = true;
3191 13150 : 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 641 : 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 641 : oid gid = 0;
3252 641 : BUN i, nils;
3253 641 : int t;
3254 641 : const void *nil;
3255 641 : int (*atomcmp)(const void *, const void *);
3256 :
3257 641 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
3258 :
3259 643 : nils = ngrp;
3260 61373 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx)
3261 60089 : oids[i] = oid_nil;
3262 642 : if (ci->ncand == 0)
3263 : return nils;
3264 :
3265 642 : t = bi->b->ttype;
3266 642 : nil = ATOMnilptr(t);
3267 642 : atomcmp = ATOMcompare(t);
3268 642 : t = ATOMbasetype(t);
3269 642 : oid hseq = bi->b->hseqbase;
3270 :
3271 642 : switch (t) {
3272 17 : case TYPE_bte:
3273 91 : AGGR_CMP(bte, LT);
3274 : break;
3275 1 : case TYPE_sht:
3276 1756 : AGGR_CMP(sht, LT);
3277 : break;
3278 484 : case TYPE_int:
3279 61224 : AGGR_CMP(int, LT);
3280 : break;
3281 56 : case TYPE_lng:
3282 95889 : 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 43 : case TYPE_dbl:
3293 354 : 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 41 : default:
3320 41 : assert(bi->b->ttype != TYPE_oid);
3321 :
3322 41 : 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 20568 : TIMEOUT_LOOP(ci->ncand, qry_ctx) {
3335 20502 : i = canditer_next(ci) - hseq;
3336 20502 : if (gids == NULL ||
3337 45 : (gids[i] >= min && gids[i] <= max)) {
3338 20502 : const void *v = BUNtail(*bi, i);
3339 20502 : if (gids)
3340 45 : gid = gids[i] - min;
3341 41004 : if (!skip_nils ||
3342 20502 : (*atomcmp)(v, nil) != 0) {
3343 18118 : if (is_oid_nil(oids[gid])) {
3344 50 : oids[gid] = i + hseq;
3345 50 : nils--;
3346 18068 : } else if (t != TYPE_void) {
3347 18068 : const void *g = BUNtail(*bi, (BUN) (oids[gid] - hseq));
3348 36136 : if ((*atomcmp)(g, nil) != 0 &&
3349 36136 : ((*atomcmp)(v, nil) == 0 ||
3350 18068 : LT((*atomcmp)(v, g), 0)))
3351 77 : oids[gid] = i + hseq;
3352 : }
3353 : }
3354 : }
3355 : }
3356 : }
3357 : break;
3358 : }
3359 640 : TIMEOUT_CHECK(qry_ctx, TIMEOUT_HANDLER(BUN_NONE, qry_ctx));
3360 :
3361 : return nils;
3362 : }
3363 :
3364 : /* calculate group maximums with optional candidates list
3365 : *
3366 : * note that this functions returns *positions* of where the maximum
3367 : * values occur */
3368 : static BUN
3369 661 : do_groupmax(oid *restrict oids, BATiter *bi, const oid *restrict gids, BUN ngrp,
3370 : oid min, oid max, struct canditer *restrict ci,
3371 : bool skip_nils, bool gdense)
3372 : {
3373 661 : oid gid = 0;
3374 661 : BUN i, nils;
3375 661 : int t;
3376 661 : const void *nil;
3377 661 : int (*atomcmp)(const void *, const void *);
3378 :
3379 661 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
3380 :
3381 661 : nils = ngrp;
3382 700959 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx)
3383 699562 : oids[i] = oid_nil;
3384 661 : if (ci->ncand == 0)
3385 : return nils;
3386 :
3387 661 : t = bi->b->ttype;
3388 661 : nil = ATOMnilptr(t);
3389 661 : atomcmp = ATOMcompare(t);
3390 661 : t = ATOMbasetype(t);
3391 661 : oid hseq = bi->b->hseqbase;
3392 :
3393 661 : switch (t) {
3394 24 : case TYPE_bte:
3395 128 : AGGR_CMP(bte, GT);
3396 : break;
3397 1 : case TYPE_sht:
3398 168 : AGGR_CMP(sht, GT);
3399 : break;
3400 488 : case TYPE_int:
3401 126183 : AGGR_CMP(int, GT);
3402 : break;
3403 73 : case TYPE_lng:
3404 94563 : AGGR_CMP(lng, GT);
3405 : break;
3406 : #ifdef HAVE_HGE
3407 5 : case TYPE_hge:
3408 35489608 : AGGR_CMP(hge, GT);
3409 : break;
3410 : #endif
3411 1 : case TYPE_flt:
3412 8 : AGGR_CMP(flt, GT);
3413 : break;
3414 34 : case TYPE_dbl:
3415 277 : 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 35 : default:
3441 35 : assert(bi->b->ttype != TYPE_oid);
3442 :
3443 35 : 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 21197 : TIMEOUT_LOOP(ci->ncand, qry_ctx) {
3456 21147 : i = canditer_next(ci) - hseq;
3457 21147 : if (gids == NULL ||
3458 20 : (gids[i] >= min && gids[i] <= max)) {
3459 21147 : const void *v = BUNtail(*bi, i);
3460 21147 : if (gids)
3461 20 : gid = gids[i] - min;
3462 42294 : if (!skip_nils ||
3463 21147 : (*atomcmp)(v, nil) != 0) {
3464 18765 : if (is_oid_nil(oids[gid])) {
3465 30 : oids[gid] = i + hseq;
3466 30 : nils--;
3467 : } else {
3468 18735 : const void *g = BUNtail(*bi, (BUN) (oids[gid] - hseq));
3469 18735 : if (t == TYPE_void ||
3470 37470 : ((*atomcmp)(g, nil) != 0 &&
3471 37470 : ((*atomcmp)(v, nil) == 0 ||
3472 18735 : GT((*atomcmp)(v, g), 0))))
3473 206 : oids[gid] = i + hseq;
3474 : }
3475 : }
3476 : }
3477 : }
3478 : }
3479 : break;
3480 : }
3481 660 : TIMEOUT_CHECK(qry_ctx, TIMEOUT_HANDLER(BUN_NONE, qry_ctx));
3482 :
3483 : return nils;
3484 : }
3485 :
3486 : static BAT *
3487 1240 : 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 1240 : const oid *restrict gids;
3494 1240 : oid min, max;
3495 1240 : BUN ngrp;
3496 1240 : oid *restrict oids;
3497 1240 : BAT *bn = NULL;
3498 1240 : BUN nils;
3499 1240 : struct canditer ci;
3500 1240 : const char *err;
3501 1240 : lng t0 = 0;
3502 :
3503 1240 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
3504 :
3505 1240 : assert(tp == TYPE_oid);
3506 1240 : (void) tp; /* compatibility (with other BATgroup* */
3507 :
3508 1240 : 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 1240 : 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 1241 : if (ci.ncand == 0 || ngrp == 0) {
3520 : /* trivial: no minimums, so return bat aligned with g
3521 : * with nil in the tail */
3522 124 : return BATconstant(ngrp == 0 ? 0 : min, TYPE_oid, &oid_nil, ngrp, TRANSIENT);
3523 : }
3524 :
3525 1117 : bn = COLnew(min, TYPE_oid, ngrp, TRANSIENT);
3526 1117 : if (bn == NULL)
3527 : return NULL;
3528 1117 : oids = (oid *) Tloc(bn, 0);
3529 :
3530 1117 : if (g == NULL || BATtdense(g))
3531 : gids = NULL;
3532 : else
3533 527 : gids = (const oid *) Tloc(g, 0);
3534 :
3535 1117 : BATiter bi = bat_iterator(b);
3536 1118 : nils = (*minmax)(oids, &bi, gids, ngrp, min, max, &ci,
3537 1118 : skip_nils, g && BATtdense(g));
3538 1116 : bat_iterator_end(&bi);
3539 1119 : if (nils == BUN_NONE) {
3540 0 : BBPreclaim(bn);
3541 0 : return NULL;
3542 : }
3543 :
3544 1119 : BATsetcount(bn, ngrp);
3545 :
3546 1119 : bn->tkey = BATcount(bn) <= 1;
3547 1119 : bn->tsorted = BATcount(bn) <= 1;
3548 1119 : bn->trevsorted = BATcount(bn) <= 1;
3549 1119 : bn->tnil = nils != 0;
3550 1119 : bn->tnonil = nils == 0;
3551 1119 : 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 599 : BATgroupmin(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
3562 : {
3563 599 : 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 1090 : BATmin_skipnil(BAT *b, void *aggr, bit skipnil)
3571 : {
3572 1090 : const void *res = NULL;
3573 1090 : size_t s;
3574 1090 : lng t0 = 0;
3575 :
3576 1090 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
3577 :
3578 1090 : 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 1090 : BATiter bi = bat_iterator(b);
3585 1090 : if (bi.count == 0) {
3586 225 : res = ATOMnilptr(bi.type);
3587 865 : } else if (bi.minpos != BUN_NONE) {
3588 381 : res = BUNtail(bi, bi.minpos);
3589 : } else {
3590 484 : oid pos;
3591 617 : BAT *pb = BATdescriptor(VIEWtparent(b));
3592 484 : Heap *oidxh = NULL;
3593 484 : bool usepoidx = false;
3594 :
3595 484 : if (BATordered(b)) {
3596 392 : if (skipnil && !bi.nonil) {
3597 304 : pos = binsearch(NULL, 0, bi.type, bi.base,
3598 304 : bi.vh ? bi.vh->base : NULL,
3599 304 : bi.width, 0, bi.count,
3600 304 : ATOMnilptr(bi.type), 1, 1);
3601 304 : if (pos == bi.count)
3602 44 : pos = oid_nil;
3603 : else
3604 260 : pos += b->hseqbase;
3605 : } else {
3606 88 : pos = b->hseqbase;
3607 : }
3608 92 : } else if (BATordered_rev(b)) {
3609 8 : if (skipnil && !bi.nonil) {
3610 4 : pos = binsearch(NULL, 0, bi.type, bi.base,
3611 4 : bi.vh ? bi.vh->base : NULL,
3612 4 : bi.width, 0, bi.count,
3613 4 : ATOMnilptr(bi.type), -1, 1);
3614 4 : if (pos == 0)
3615 0 : pos = oid_nil;
3616 : else
3617 4 : pos = pos + b->hseqbase - 1;
3618 : } else {
3619 4 : pos = bi.count + b->hseqbase - 1;
3620 : }
3621 : } else {
3622 84 : 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 84 : if (oidxh == NULL &&
3630 124 : pb != NULL &&
3631 40 : 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 84 : 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 84 : Imprints *imprints = NULL;
3675 144 : if ((pb == NULL || bi.count == BATcount(pb)) &&
3676 60 : BATcheckimprints(b)) {
3677 0 : if (pb != NULL) {
3678 0 : MT_lock_set(&pb->batIdxLock);
3679 0 : imprints = pb->timprints;
3680 0 : if (imprints != NULL)
3681 0 : IMPSincref(imprints);
3682 : else
3683 : imprints = NULL;
3684 0 : MT_lock_unset(&pb->batIdxLock);
3685 : } else {
3686 0 : MT_lock_set(&b->batIdxLock);
3687 0 : imprints = b->timprints;
3688 0 : if (imprints != NULL)
3689 0 : IMPSincref(imprints);
3690 : else
3691 : imprints = NULL;
3692 0 : MT_lock_unset(&b->batIdxLock);
3693 : }
3694 : }
3695 0 : if (imprints) {
3696 0 : int i;
3697 :
3698 0 : MT_thread_setalgorithm(VIEWtparent(b) ? "using parent imprints" : "using imprints");
3699 0 : pos = oid_nil;
3700 : /* find first non-empty bin */
3701 0 : for (i = 0; i < imprints->bits; i++) {
3702 0 : if (imprints->stats[i + 128]) {
3703 0 : pos = imprints->stats[i] + b->hseqbase;
3704 0 : break;
3705 : }
3706 : }
3707 0 : IMPSdecref(imprints, false);
3708 : } else {
3709 84 : struct canditer ci;
3710 84 : canditer_init(&ci, b, NULL);
3711 84 : (void) do_groupmin(&pos, &bi, NULL, 1, 0, 0, &ci,
3712 : skipnil, false);
3713 : }
3714 : }
3715 : }
3716 484 : if (is_oid_nil(pos)) {
3717 44 : res = ATOMnilptr(bi.type);
3718 : } else {
3719 440 : bi.minpos = pos - b->hseqbase;
3720 440 : res = BUNtail(bi, bi.minpos);
3721 440 : MT_lock_set(&b->theaplock);
3722 440 : if (bi.count == BATcount(b) && bi.h == b->theap)
3723 440 : b->tminpos = bi.minpos;
3724 440 : MT_lock_unset(&b->theaplock);
3725 440 : if (pb) {
3726 337 : MT_lock_set(&pb->theaplock);
3727 337 : if (bi.count == BATcount(pb) &&
3728 86 : bi.h == pb->theap)
3729 86 : pb->tminpos = bi.minpos;
3730 337 : MT_lock_unset(&pb->theaplock);
3731 : }
3732 : }
3733 835 : BBPreclaim(pb);
3734 : }
3735 1090 : if (aggr == NULL) {
3736 490 : s = ATOMlen(bi.type, res);
3737 490 : aggr = GDKmalloc(s);
3738 : } else {
3739 1200 : s = ATOMsize(ATOMtype(bi.type));
3740 : }
3741 1090 : if (aggr != NULL) /* else: malloc error */
3742 1091 : memcpy(aggr, res, s);
3743 1090 : bat_iterator_end(&bi);
3744 1091 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",skipnil=%d; (" LLFMT " usec)\n",
3745 : ALGOBATPAR(b), skipnil, GDKusec() - t0);
3746 : return aggr;
3747 : }
3748 :
3749 : void *
3750 464 : BATmin(BAT *b, void *aggr)
3751 : {
3752 464 : return BATmin_skipnil(b, aggr, 1);
3753 : }
3754 :
3755 : BAT *
3756 641 : BATgroupmax(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
3757 : {
3758 641 : return BATgroupminmax(b, g, e, s, tp, skip_nils,
3759 : do_groupmax, __func__);
3760 : }
3761 :
3762 : void *
3763 1409 : BATmax_skipnil(BAT *b, void *aggr, bit skipnil)
3764 : {
3765 1409 : const void *res = NULL;
3766 1409 : size_t s;
3767 1409 : lng t0 = 0;
3768 :
3769 1409 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
3770 :
3771 1409 : if (!ATOMlinear(b->ttype)) {
3772 : /* there is no such thing as a largest value if you
3773 : * can't compare values */
3774 0 : GDKerror("non-linear type");
3775 0 : return NULL;
3776 : }
3777 1409 : BATiter bi = bat_iterator(b);
3778 1411 : if (bi.count == 0) {
3779 201 : res = ATOMnilptr(bi.type);
3780 1210 : } else if (bi.maxpos != BUN_NONE) {
3781 727 : res = BUNtail(bi, bi.maxpos);
3782 : } else {
3783 483 : oid pos;
3784 622 : BAT *pb = BATdescriptor(VIEWtparent(b));
3785 483 : Heap *oidxh = NULL;
3786 483 : bool usepoidx = false;
3787 :
3788 483 : if (BATordered(b)) {
3789 373 : pos = bi.count - 1 + b->hseqbase;
3790 660 : if (skipnil && !bi.nonil &&
3791 287 : ATOMcmp(bi.type, BUNtail(bi, bi.count - 1),
3792 : ATOMnilptr(bi.type)) == 0)
3793 38 : pos = oid_nil; /* no non-nil values */
3794 110 : } else if (BATordered_rev(b)) {
3795 9 : pos = b->hseqbase;
3796 14 : if (skipnil && !bi.nonil &&
3797 5 : ATOMcmp(bi.type, BUNtail(bi, 0),
3798 : ATOMnilptr(bi.type)) == 0)
3799 0 : pos = oid_nil; /* no non-nil values */
3800 : } else {
3801 101 : if (BATcheckorderidx(b)) {
3802 0 : MT_lock_set(&b->batIdxLock);
3803 0 : oidxh = b->torderidx;
3804 0 : if (oidxh != NULL)
3805 0 : HEAPincref(oidxh);
3806 0 : MT_lock_unset(&b->batIdxLock);
3807 : }
3808 101 : if (oidxh == NULL &&
3809 147 : pb != NULL &&
3810 46 : BATcheckorderidx(pb)) {
3811 : /* no lock on b needed since it's a view */
3812 0 : MT_lock_set(&pb->batIdxLock);
3813 0 : MT_lock_set(&pb->theaplock);
3814 0 : if (pb->tbaseoff == bi.baseoff &&
3815 0 : BATcount(pb) == bi.count &&
3816 0 : pb->hseqbase == b->hseqbase &&
3817 0 : (oidxh = pb->torderidx) != NULL) {
3818 0 : HEAPincref(oidxh);
3819 0 : usepoidx = true;
3820 : }
3821 0 : MT_lock_unset(&pb->batIdxLock);
3822 0 : MT_lock_unset(&pb->theaplock);
3823 : }
3824 101 : if (oidxh != NULL) {
3825 0 : const oid *ords = (const oid *) oidxh->base + ORDERIDXOFF;
3826 :
3827 0 : MT_thread_setalgorithm(usepoidx ? "using parent oidx" : "using oids");
3828 0 : pos = ords[bi.count - 1];
3829 : /* nils are first, ie !skipnil, check for nils */
3830 0 : if (!skipnil) {
3831 0 : BUN z = ords[0];
3832 :
3833 0 : res = BUNtail(bi, z - b->hseqbase);
3834 :
3835 0 : if (ATOMcmp(bi.type, res, ATOMnilptr(bi.type)) == 0)
3836 0 : pos = z;
3837 : }
3838 0 : HEAPdecref(oidxh, false);
3839 : } else {
3840 101 : Imprints *imprints = NULL;
3841 170 : if ((pb == NULL || BATcount(b) == BATcount(pb)) &&
3842 69 : BATcheckimprints(b)) {
3843 0 : if (pb != NULL) {
3844 0 : MT_lock_set(&pb->batIdxLock);
3845 0 : imprints = pb->timprints;
3846 0 : if (imprints != NULL)
3847 0 : IMPSincref(imprints);
3848 : else
3849 : imprints = NULL;
3850 0 : MT_lock_unset(&pb->batIdxLock);
3851 : } else {
3852 0 : MT_lock_set(&b->batIdxLock);
3853 0 : imprints = b->timprints;
3854 0 : if (imprints != NULL)
3855 0 : IMPSincref(imprints);
3856 : else
3857 : imprints = NULL;
3858 0 : MT_lock_unset(&b->batIdxLock);
3859 : }
3860 : }
3861 0 : if (imprints) {
3862 0 : int i;
3863 :
3864 0 : MT_thread_setalgorithm(VIEWtparent(b) ? "using parent imprints" : "using imprints");
3865 0 : pos = oid_nil;
3866 : /* find last non-empty bin */
3867 0 : for (i = imprints->bits - 1; i >= 0; i--) {
3868 0 : if (imprints->stats[i + 128]) {
3869 0 : pos = imprints->stats[i + 64] + b->hseqbase;
3870 0 : break;
3871 : }
3872 : }
3873 0 : IMPSdecref(imprints, false);
3874 : } else {
3875 101 : struct canditer ci;
3876 101 : canditer_init(&ci, b, NULL);
3877 101 : (void) do_groupmax(&pos, &bi, NULL, 1, 0, 0, &ci,
3878 : skipnil, false);
3879 : }
3880 : }
3881 : }
3882 482 : if (is_oid_nil(pos)) {
3883 38 : res = ATOMnilptr(bi.type);
3884 : } else {
3885 444 : bi.maxpos = pos - b->hseqbase;
3886 444 : res = BUNtail(bi, bi.maxpos);
3887 443 : MT_lock_set(&b->theaplock);
3888 445 : if (bi.count == BATcount(b) && bi.h == b->theap)
3889 445 : b->tmaxpos = bi.maxpos;
3890 445 : MT_lock_unset(&b->theaplock);
3891 445 : if (pb) {
3892 336 : MT_lock_set(&pb->theaplock);
3893 336 : if (bi.count == BATcount(pb) &&
3894 111 : bi.h == pb->theap)
3895 111 : pb->tmaxpos = bi.maxpos;
3896 336 : MT_lock_unset(&pb->theaplock);
3897 : }
3898 : }
3899 827 : BBPreclaim(pb);
3900 : }
3901 1411 : if (aggr == NULL) {
3902 514 : s = ATOMlen(bi.type, res);
3903 514 : aggr = GDKmalloc(s);
3904 : } else {
3905 1793 : s = ATOMsize(ATOMtype(bi.type));
3906 : }
3907 1412 : if (aggr != NULL) /* else: malloc error */
3908 1411 : memcpy(aggr, res, s);
3909 1412 : bat_iterator_end(&bi);
3910 1411 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",skipnil=%d; (" LLFMT " usec)\n",
3911 : ALGOBATPAR(b), skipnil, GDKusec() - t0);
3912 : return aggr;
3913 : }
3914 :
3915 : void *
3916 778 : BATmax(BAT *b, void *aggr)
3917 : {
3918 778 : return BATmax_skipnil(b, aggr, 1);
3919 : }
3920 :
3921 :
3922 : /* ---------------------------------------------------------------------- */
3923 : /* quantiles/median */
3924 :
3925 : #if SIZEOF_OID == SIZEOF_INT
3926 : #define binsearch_oid(indir, offset, vals, lo, hi, v, ordering, last) binsearch_int(indir, offset, (const int *) vals, lo, hi, (int) (v), ordering, last)
3927 : #endif
3928 : #if SIZEOF_OID == SIZEOF_LNG
3929 : #define binsearch_oid(indir, offset, vals, lo, hi, v, ordering, last) binsearch_lng(indir, offset, (const lng *) vals, lo, hi, (lng) (v), ordering, last)
3930 : #endif
3931 :
3932 : #define DO_QUANTILE_AVG(TPE) \
3933 : do { \
3934 : BUN idxlo, idxhi; \
3935 : if (ords) { \
3936 : idxlo = ords[r + (BUN) lo] - b->hseqbase; \
3937 : idxhi = ords[r + (BUN) hi] - b->hseqbase; \
3938 : } else { \
3939 : idxlo = r + (BUN) lo; \
3940 : idxhi = r + (BUN) hi; \
3941 : } \
3942 : TPE low = *(TPE*) BUNtloc(bi, idxhi); \
3943 : TPE high = *(TPE*) BUNtloc(bi, idxlo); \
3944 : if (is_##TPE##_nil(low) || is_##TPE##_nil(high)) { \
3945 : val = dbl_nil; \
3946 : nils++; \
3947 : } else { \
3948 : val = (f - lo) * low + (lo + 1 - f) * high; \
3949 : } \
3950 : } while (0)
3951 :
3952 : static BAT *
3953 65 : doBATgroupquantile(BAT *b, BAT *g, BAT *e, BAT *s, int tp, double quantile,
3954 : bool skip_nils, bool average)
3955 : {
3956 65 : BAT *origb = b;
3957 65 : BAT *origg = g;
3958 65 : oid min, max;
3959 65 : BUN ngrp;
3960 65 : BUN nils = 0;
3961 65 : BAT *bn = NULL;
3962 65 : struct canditer ci;
3963 65 : BAT *t1, *t2;
3964 65 : BATiter bi = {0};
3965 65 : const void *v;
3966 65 : const void *nil = ATOMnilptr(tp);
3967 65 : const void *dnil = nil;
3968 65 : dbl val; /* only used for average */
3969 65 : int (*atomcmp)(const void *, const void *) = ATOMcompare(tp);
3970 65 : const char *err;
3971 65 : lng t0 = 0;
3972 :
3973 65 : size_t counter = 0;
3974 65 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
3975 :
3976 65 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
3977 :
3978 65 : if (average) {
3979 12 : switch (ATOMbasetype(b->ttype)) {
3980 : case TYPE_bte:
3981 : case TYPE_sht:
3982 : case TYPE_int:
3983 : case TYPE_lng:
3984 : #ifdef HAVE_HGE
3985 : case TYPE_hge:
3986 : #endif
3987 : case TYPE_flt:
3988 : case TYPE_dbl:
3989 : break;
3990 0 : default:
3991 0 : GDKerror("incompatible type\n");
3992 0 : return NULL;
3993 : }
3994 : dnil = &dbl_nil;
3995 : }
3996 65 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
3997 0 : GDKerror("%s\n", err);
3998 0 : return NULL;
3999 : }
4000 65 : assert(tp == b->ttype);
4001 65 : if (!ATOMlinear(tp)) {
4002 0 : GDKerror("cannot determine quantile on "
4003 : "non-linear type %s\n", ATOMname(tp));
4004 0 : return NULL;
4005 : }
4006 65 : if (quantile < 0 || quantile > 1) {
4007 2 : GDKerror("cannot determine quantile for "
4008 : "p=%f (p has to be in [0,1])\n", quantile);
4009 2 : return NULL;
4010 : }
4011 :
4012 63 : if (BATcount(b) == 0 || ngrp == 0 || is_dbl_nil(quantile)) {
4013 : /* trivial: no values, thus also no quantiles,
4014 : * so return bat aligned with e with nil in the tail
4015 : * The same happens for a NULL quantile */
4016 24 : return BATconstant(ngrp == 0 ? 0 : min, average ? TYPE_dbl : tp, dnil, ngrp, TRANSIENT);
4017 : }
4018 :
4019 51 : if (s) {
4020 : /* there is a candidate list, replace b (and g, if
4021 : * given) with just the values we're interested in */
4022 0 : b = BATproject(s, b);
4023 0 : if (b == NULL)
4024 : return NULL;
4025 0 : if (g) {
4026 0 : g = BATproject(s, g);
4027 0 : if (g == NULL)
4028 0 : goto bailout;
4029 : }
4030 : }
4031 :
4032 : /* we want to sort b so that we can figure out the quantile, but
4033 : * if g is given, sort g and subsort b so that we can get the
4034 : * quantile for each group */
4035 51 : if (g) {
4036 19 : const oid *restrict grps;
4037 19 : oid prev;
4038 19 : BUN p, q, r;
4039 :
4040 19 : if (BATtdense(g)) {
4041 : /* singleton groups, so calculating quantile is
4042 : * easy */
4043 1 : if (average)
4044 0 : bn = BATconvert(b, NULL, TYPE_dbl, 0, 0, 0);
4045 : else
4046 1 : bn = COLcopy(b, tp, false, TRANSIENT);
4047 1 : BAThseqbase(bn, g->tseqbase); /* deals with NULL */
4048 1 : if (b != origb)
4049 0 : BBPunfix(b->batCacheid);
4050 1 : if (g != origg)
4051 0 : BBPunfix(g->batCacheid);
4052 1 : return bn;
4053 : }
4054 18 : if (BATsort(&t1, &t2, NULL, g, NULL, NULL, false, false, false) != GDK_SUCCEED)
4055 0 : goto bailout;
4056 18 : if (g != origg)
4057 0 : BBPunfix(g->batCacheid);
4058 18 : g = t1;
4059 :
4060 18 : if (BATsort(&t1, NULL, NULL, b, t2, g, false, false, false) != GDK_SUCCEED) {
4061 0 : BBPunfix(t2->batCacheid);
4062 0 : goto bailout;
4063 : }
4064 18 : if (b != origb)
4065 0 : BBPunfix(b->batCacheid);
4066 18 : b = t1;
4067 18 : BBPunfix(t2->batCacheid);
4068 :
4069 18 : if (average)
4070 2 : bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
4071 : else
4072 16 : bn = COLnew(min, tp, ngrp, TRANSIENT);
4073 18 : if (bn == NULL)
4074 0 : goto bailout;
4075 :
4076 18 : bi = bat_iterator(b);
4077 :
4078 18 : grps = (const oid *) Tloc(g, 0);
4079 : /* for each group (r and p are the beginning and end
4080 : * of the current group, respectively) */
4081 68 : for (r = 0, q = BATcount(g); r < q; r = p) {
4082 50 : GDK_CHECK_TIMEOUT(qry_ctx, counter,
4083 : GOTO_LABEL_TIMEOUT_HANDLER(bunins_failed, qry_ctx));
4084 50 : BUN qindex;
4085 50 : prev = grps[r];
4086 : /* search for end of current group (grps is
4087 : * sorted so we can use binary search) */
4088 50 : p = binsearch_oid(NULL, 0, grps, r, q - 1, prev, 1, 1);
4089 50 : if (skip_nils && !bi.nonil) {
4090 : /* within group, locate start of non-nils */
4091 23 : r = binsearch(NULL, 0, tp, bi.base,
4092 23 : bi.vh ? bi.vh->base : NULL,
4093 23 : bi.width, r, p, nil,
4094 : 1, 1);
4095 : }
4096 50 : if (r == p) {
4097 : /* no non-nils found */
4098 8 : v = dnil;
4099 8 : nils++;
4100 42 : } else if (average) {
4101 4 : double f = (p - r - 1) * quantile;
4102 4 : double lo = floor(f);
4103 4 : double hi = ceil(f);
4104 4 : const oid *const ords = NULL;
4105 8 : switch (ATOMbasetype(tp)) {
4106 : case TYPE_bte:
4107 4 : DO_QUANTILE_AVG(bte);
4108 : break;
4109 : case TYPE_sht:
4110 0 : DO_QUANTILE_AVG(sht);
4111 : break;
4112 : case TYPE_int:
4113 0 : DO_QUANTILE_AVG(int);
4114 : break;
4115 : case TYPE_lng:
4116 0 : DO_QUANTILE_AVG(lng);
4117 : break;
4118 : #ifdef HAVE_HGE
4119 : case TYPE_hge:
4120 0 : DO_QUANTILE_AVG(hge);
4121 : break;
4122 : #endif
4123 : case TYPE_flt:
4124 0 : DO_QUANTILE_AVG(flt);
4125 : break;
4126 : case TYPE_dbl:
4127 0 : DO_QUANTILE_AVG(dbl);
4128 : break;
4129 : }
4130 : v = &val;
4131 : } else {
4132 : /* round *down* to nearest integer */
4133 38 : double f = (p - r - 1) * quantile;
4134 38 : qindex = r + p - (BUN) (p + 0.5 - f);
4135 : /* be a little paranoid about the index */
4136 38 : assert(qindex >= r && qindex < p);
4137 38 : v = BUNtail(bi, qindex);
4138 38 : if (!skip_nils && !bi.nonil)
4139 0 : nils += (*atomcmp)(v, dnil) == 0;
4140 : }
4141 50 : if (bunfastapp_nocheck(bn, v) != GDK_SUCCEED)
4142 0 : goto bunins_failed;
4143 : }
4144 18 : bat_iterator_end(&bi);
4145 18 : nils += ngrp - BATcount(bn);
4146 18 : while (BATcount(bn) < ngrp) {
4147 0 : if (bunfastapp_nocheck(bn, dnil) != GDK_SUCCEED)
4148 0 : goto bailout;
4149 : }
4150 18 : bn->theap->dirty = true;
4151 18 : BBPunfix(g->batCacheid);
4152 : } else {
4153 32 : BUN index, r, p = BATcount(b);
4154 32 : BAT *pb = NULL;
4155 32 : const oid *ords;
4156 32 : Heap *oidxh = NULL;
4157 :
4158 60 : bn = COLnew(0, average ? TYPE_dbl : tp, 1, TRANSIENT);
4159 32 : if (bn == NULL)
4160 0 : goto bailout;
4161 :
4162 32 : t1 = NULL;
4163 :
4164 32 : if (BATcheckorderidx(b)) {
4165 0 : MT_lock_set(&b->batIdxLock);
4166 0 : oidxh = b->torderidx;
4167 0 : if (oidxh != NULL)
4168 0 : HEAPincref(oidxh);
4169 0 : MT_lock_unset(&b->batIdxLock);
4170 : }
4171 4 : if (oidxh == NULL &&
4172 32 : VIEWtparent(b) &&
4173 4 : (pb = BATdescriptor(VIEWtparent(b))) != NULL) {
4174 4 : if (BATcheckorderidx(pb)) {
4175 : /* no lock on b needed since it's a view */
4176 0 : MT_lock_set(&pb->batIdxLock);
4177 0 : MT_lock_set(&pb->theaplock);
4178 0 : if (pb->tbaseoff == b->tbaseoff &&
4179 0 : BATcount(pb) == BATcount(b) &&
4180 0 : pb->hseqbase == b->hseqbase &&
4181 0 : (oidxh = pb->torderidx) != NULL) {
4182 0 : HEAPincref(oidxh);
4183 : }
4184 0 : MT_lock_unset(&pb->batIdxLock);
4185 0 : MT_lock_unset(&pb->theaplock);
4186 : }
4187 4 : BBPunfix(pb->batCacheid);
4188 : }
4189 32 : if (oidxh != NULL) {
4190 0 : MT_thread_setalgorithm(pb ? "using parent oidx" : "using oids");
4191 0 : ords = (const oid *) oidxh->base + ORDERIDXOFF;
4192 : } else {
4193 32 : if (BATsort(NULL, &t1, NULL, b, NULL, g, false, false, false) != GDK_SUCCEED)
4194 0 : goto bailout;
4195 32 : if (BATtdense(t1))
4196 : ords = NULL;
4197 : else
4198 11 : ords = (const oid *) Tloc(t1, 0);
4199 : }
4200 :
4201 32 : bi = bat_iterator(b);
4202 :
4203 32 : if (skip_nils && !bi.nonil)
4204 14 : r = binsearch(ords, 0, tp, bi.base,
4205 14 : bi.vh ? bi.vh->base : NULL,
4206 14 : bi.width, 0, p,
4207 : nil, 1, 1);
4208 : else
4209 : r = 0;
4210 :
4211 32 : if (r == p) {
4212 : /* no non-nil values, so quantile is nil */
4213 : v = dnil;
4214 : nils++;
4215 29 : } else if (average) {
4216 4 : double f = (p - r - 1) * quantile;
4217 4 : double lo = floor(f);
4218 4 : double hi = ceil(f);
4219 8 : switch (ATOMbasetype(tp)) {
4220 4 : case TYPE_bte:
4221 4 : DO_QUANTILE_AVG(bte);
4222 : break;
4223 0 : case TYPE_sht:
4224 0 : DO_QUANTILE_AVG(sht);
4225 : break;
4226 0 : case TYPE_int:
4227 0 : DO_QUANTILE_AVG(int);
4228 : break;
4229 0 : case TYPE_lng:
4230 0 : DO_QUANTILE_AVG(lng);
4231 : break;
4232 : #ifdef HAVE_HGE
4233 0 : case TYPE_hge:
4234 0 : DO_QUANTILE_AVG(hge);
4235 : break;
4236 : #endif
4237 0 : case TYPE_flt:
4238 0 : DO_QUANTILE_AVG(flt);
4239 : break;
4240 0 : case TYPE_dbl:
4241 0 : DO_QUANTILE_AVG(dbl);
4242 : break;
4243 : }
4244 : v = &val;
4245 : } else {
4246 25 : double f;
4247 : /* round (p-r-1)*quantile *down* to nearest
4248 : * integer (i.e., 1.49 and 1.5 are rounded to
4249 : * 1, 1.51 is rounded to 2) */
4250 25 : f = (p - r - 1) * quantile;
4251 25 : index = r + p - (BUN) (p + 0.5 - f);
4252 25 : if (ords)
4253 10 : index = ords[index] - b->hseqbase;
4254 : else
4255 15 : index = index + t1->tseqbase;
4256 25 : v = BUNtail(bi, index);
4257 25 : nils += (*atomcmp)(v, dnil) == 0;
4258 : }
4259 32 : if (oidxh != NULL)
4260 0 : HEAPdecref(oidxh, false);
4261 32 : BBPreclaim(t1);
4262 32 : gdk_return rc = BUNappend(bn, v, false);
4263 32 : bat_iterator_end(&bi);
4264 32 : if (rc != GDK_SUCCEED)
4265 0 : goto bailout;
4266 : }
4267 :
4268 50 : if (b != origb)
4269 18 : BBPunfix(b->batCacheid);
4270 :
4271 50 : bn->tkey = BATcount(bn) <= 1;
4272 50 : bn->tsorted = BATcount(bn) <= 1;
4273 50 : bn->trevsorted = BATcount(bn) <= 1;
4274 50 : bn->tnil = nils != 0;
4275 50 : bn->tnonil = nils == 0;
4276 50 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
4277 : "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT
4278 : ",quantile=%g,average=%s -> " ALGOOPTBATFMT
4279 : "; start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
4280 : ALGOBATPAR(origb), ALGOOPTBATPAR(origg), ALGOOPTBATPAR(e),
4281 : ALGOOPTBATPAR(s), quantile, average ? "true" : "false",
4282 : ALGOOPTBATPAR(bn), ci.seq, ci.ncand, GDKusec() - t0);
4283 : return bn;
4284 :
4285 0 : bunins_failed:
4286 0 : bat_iterator_end(&bi);
4287 0 : bailout:
4288 0 : if (b && b != origb)
4289 0 : BBPunfix(b->batCacheid);
4290 0 : if (g && g != origg)
4291 0 : BBPunfix(g->batCacheid);
4292 0 : BBPreclaim(bn);
4293 : return NULL;
4294 : }
4295 :
4296 : BAT *
4297 28 : BATgroupmedian(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
4298 : bool skip_nils)
4299 : {
4300 28 : return doBATgroupquantile(b, g, e, s, tp, 0.5,
4301 : skip_nils, false);
4302 : }
4303 :
4304 : BAT *
4305 31 : BATgroupquantile(BAT *b, BAT *g, BAT *e, BAT *s, int tp, double quantile,
4306 : bool skip_nils)
4307 : {
4308 31 : return doBATgroupquantile(b, g, e, s, tp, quantile,
4309 : skip_nils, false);
4310 : }
4311 :
4312 : BAT *
4313 3 : BATgroupmedian_avg(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
4314 : {
4315 3 : return doBATgroupquantile(b, g, e, s, tp, 0.5, skip_nils, true);
4316 : }
4317 :
4318 : BAT *
4319 3 : BATgroupquantile_avg(BAT *b, BAT *g, BAT *e, BAT *s, int tp, double quantile,
4320 : bool skip_nils)
4321 : {
4322 3 : return doBATgroupquantile(b, g, e, s, tp, quantile,
4323 : skip_nils, true);
4324 : }
4325 :
4326 : /* ---------------------------------------------------------------------- */
4327 : /* standard deviation (both biased and non-biased) */
4328 :
4329 : #define AGGR_STDEV_SINGLE(TYPE) \
4330 : do { \
4331 : TYPE x; \
4332 : TIMEOUT_LOOP_IDX(i, cnt, qry_ctx) { \
4333 : x = ((const TYPE *) values)[i]; \
4334 : if (is_##TYPE##_nil(x)) \
4335 : continue; \
4336 : n++; \
4337 : delta = (dbl) x - mean; \
4338 : mean += delta / n; \
4339 : m2 += delta * ((dbl) x - mean); \
4340 : if (isinf(m2)) \
4341 : goto overflow; \
4342 : } \
4343 : TIMEOUT_CHECK(qry_ctx, \
4344 : TIMEOUT_HANDLER(dbl_nil, qry_ctx)); \
4345 : } while (0)
4346 :
4347 : static dbl
4348 25 : calcvariance(dbl *restrict avgp, const void *restrict values, BUN cnt, int tp, bool issample)
4349 : {
4350 25 : BUN n = 0, i;
4351 25 : dbl mean = 0;
4352 25 : dbl m2 = 0;
4353 25 : dbl delta;
4354 :
4355 25 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
4356 :
4357 25 : switch (tp) {
4358 : case TYPE_bte:
4359 18 : AGGR_STDEV_SINGLE(bte);
4360 : break;
4361 : case TYPE_sht:
4362 18 : AGGR_STDEV_SINGLE(sht);
4363 : break;
4364 : case TYPE_int:
4365 47 : AGGR_STDEV_SINGLE(int);
4366 : break;
4367 : case TYPE_lng:
4368 36 : AGGR_STDEV_SINGLE(lng);
4369 : break;
4370 : #ifdef HAVE_HGE
4371 : case TYPE_hge:
4372 0 : AGGR_STDEV_SINGLE(hge);
4373 : break;
4374 : #endif
4375 : case TYPE_flt:
4376 0 : AGGR_STDEV_SINGLE(flt);
4377 : break;
4378 : case TYPE_dbl:
4379 29 : AGGR_STDEV_SINGLE(dbl);
4380 : break;
4381 0 : default:
4382 0 : GDKerror("type (%s) not supported.\n", ATOMname(tp));
4383 0 : return dbl_nil;
4384 : }
4385 23 : if (n <= (BUN) issample) {
4386 7 : if (avgp)
4387 0 : *avgp = dbl_nil;
4388 7 : return dbl_nil;
4389 : }
4390 16 : if (avgp)
4391 0 : *avgp = mean;
4392 16 : return m2 / (n - issample);
4393 2 : overflow:
4394 2 : GDKerror("22003!overflow in calculation.\n");
4395 2 : return dbl_nil;
4396 : }
4397 :
4398 : dbl
4399 13 : BATcalcstdev_population(dbl *avgp, BAT *b)
4400 : {
4401 13 : lng t0 = 0;
4402 :
4403 13 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4404 13 : BATiter bi = bat_iterator(b);
4405 13 : dbl v = calcvariance(avgp, bi.base, bi.count, bi.type, false);
4406 13 : bat_iterator_end(&bi);
4407 13 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
4408 : ALGOBATPAR(b), GDKusec() - t0);
4409 13 : return is_dbl_nil(v) ? dbl_nil : sqrt(v);
4410 : }
4411 :
4412 : dbl
4413 7 : BATcalcstdev_sample(dbl *avgp, BAT *b)
4414 : {
4415 7 : lng t0 = 0;
4416 :
4417 7 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4418 7 : BATiter bi = bat_iterator(b);
4419 7 : dbl v = calcvariance(avgp, bi.base, bi.count, bi.type, true);
4420 7 : bat_iterator_end(&bi);
4421 7 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
4422 : ALGOBATPAR(b), GDKusec() - t0);
4423 7 : return is_dbl_nil(v) ? dbl_nil : sqrt(v);
4424 : }
4425 :
4426 : dbl
4427 3 : BATcalcvariance_population(dbl *avgp, BAT *b)
4428 : {
4429 3 : lng t0 = 0;
4430 :
4431 3 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4432 3 : BATiter bi = bat_iterator(b);
4433 3 : dbl v = calcvariance(avgp, bi.base, bi.count, bi.type, false);
4434 3 : bat_iterator_end(&bi);
4435 3 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
4436 : ALGOBATPAR(b), GDKusec() - t0);
4437 3 : return v;
4438 : }
4439 :
4440 : dbl
4441 2 : BATcalcvariance_sample(dbl *avgp, BAT *b)
4442 : {
4443 2 : lng t0 = 0;
4444 :
4445 2 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4446 2 : BATiter bi = bat_iterator(b);
4447 2 : dbl v = calcvariance(avgp, bi.base, bi.count, bi.type, true);
4448 2 : bat_iterator_end(&bi);
4449 2 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
4450 : ALGOBATPAR(b), GDKusec() - t0);
4451 2 : return v;
4452 : }
4453 :
4454 : #define AGGR_COVARIANCE_SINGLE(TYPE) \
4455 : do { \
4456 : TYPE x, y; \
4457 : TIMEOUT_LOOP_IDX(i, cnt, qry_ctx) { \
4458 : x = ((const TYPE *) v1)[i]; \
4459 : y = ((const TYPE *) v2)[i]; \
4460 : if (is_##TYPE##_nil(x) || is_##TYPE##_nil(y)) \
4461 : continue; \
4462 : n++; \
4463 : delta1 = (dbl) x - mean1; \
4464 : mean1 += delta1 / n; \
4465 : delta2 = (dbl) y - mean2; \
4466 : mean2 += delta2 / n; \
4467 : m2 += delta1 * ((dbl) y - mean2); \
4468 : if (isinf(m2)) \
4469 : goto overflow; \
4470 : } \
4471 : TIMEOUT_CHECK(qry_ctx, \
4472 : TIMEOUT_HANDLER(dbl_nil, qry_ctx)); \
4473 : } while (0)
4474 :
4475 : static dbl
4476 13 : calccovariance(const void *v1, const void *v2, BUN cnt, int tp, bool issample)
4477 : {
4478 13 : BUN n = 0, i;
4479 13 : dbl mean1 = 0, mean2 = 0, m2 = 0, delta1, delta2;
4480 :
4481 13 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
4482 :
4483 :
4484 13 : switch (tp) {
4485 : case TYPE_bte:
4486 0 : AGGR_COVARIANCE_SINGLE(bte);
4487 : break;
4488 : case TYPE_sht:
4489 0 : AGGR_COVARIANCE_SINGLE(sht);
4490 : break;
4491 : case TYPE_int:
4492 111 : AGGR_COVARIANCE_SINGLE(int);
4493 : break;
4494 : case TYPE_lng:
4495 0 : AGGR_COVARIANCE_SINGLE(lng);
4496 : break;
4497 : #ifdef HAVE_HGE
4498 : case TYPE_hge:
4499 0 : AGGR_COVARIANCE_SINGLE(hge);
4500 : break;
4501 : #endif
4502 : case TYPE_flt:
4503 0 : AGGR_COVARIANCE_SINGLE(flt);
4504 : break;
4505 : case TYPE_dbl:
4506 14 : AGGR_COVARIANCE_SINGLE(dbl);
4507 : break;
4508 0 : default:
4509 0 : GDKerror("type (%s) not supported.\n", ATOMname(tp));
4510 0 : return dbl_nil;
4511 : }
4512 12 : if (n <= (BUN) issample)
4513 1 : return dbl_nil;
4514 11 : return m2 / (n - issample);
4515 1 : overflow:
4516 1 : GDKerror("22003!overflow in calculation.\n");
4517 1 : return dbl_nil;
4518 : }
4519 :
4520 : dbl
4521 8 : BATcalccovariance_population(BAT *b1, BAT *b2)
4522 : {
4523 8 : lng t0 = 0;
4524 :
4525 8 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4526 8 : BATiter b1i = bat_iterator(b1), b2i = bat_iterator(b2);
4527 8 : dbl v = calccovariance(b1i.base, b2i.base, b1i.count, b1i.type, false);
4528 8 : bat_iterator_end(&b1i);
4529 8 : bat_iterator_end(&b2i);
4530 8 : TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT " (" LLFMT " usec)\n",
4531 : ALGOBATPAR(b1), ALGOBATPAR(b2), GDKusec() - t0);
4532 8 : return v;
4533 : }
4534 :
4535 : dbl
4536 5 : BATcalccovariance_sample(BAT *b1, BAT *b2)
4537 : {
4538 5 : lng t0 = 0;
4539 :
4540 5 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4541 5 : BATiter b1i = bat_iterator(b1), b2i = bat_iterator(b2);
4542 5 : dbl v = calccovariance(b1i.base, b2i.base, b1i.count, b1i.type, true);
4543 5 : bat_iterator_end(&b1i);
4544 5 : bat_iterator_end(&b2i);
4545 5 : TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT " (" LLFMT " usec)\n",
4546 : ALGOBATPAR(b1), ALGOBATPAR(b2), GDKusec() - t0);
4547 5 : return v;
4548 : }
4549 :
4550 : #define AGGR_CORRELATION_SINGLE(TYPE) \
4551 : do { \
4552 : TYPE x, y; \
4553 : TIMEOUT_LOOP_IDX(i, cnt, qry_ctx) { \
4554 : x = ((const TYPE *) v1)[i]; \
4555 : y = ((const TYPE *) v2)[i]; \
4556 : if (is_##TYPE##_nil(x) || is_##TYPE##_nil(y)) \
4557 : continue; \
4558 : n++; \
4559 : delta1 = (dbl) x - mean1; \
4560 : mean1 += delta1 / n; \
4561 : delta2 = (dbl) y - mean2; \
4562 : mean2 += delta2 / n; \
4563 : aux = (dbl) y - mean2; \
4564 : up += delta1 * aux; \
4565 : down1 += delta1 * ((dbl) x - mean1); \
4566 : down2 += delta2 * aux; \
4567 : if (isinf(up) || isinf(down1) || isinf(down2)) \
4568 : goto overflow; \
4569 : } \
4570 : TIMEOUT_CHECK(qry_ctx, \
4571 : GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx)); \
4572 : } while (0)
4573 :
4574 : dbl
4575 18 : BATcalccorrelation(BAT *b1, BAT *b2)
4576 : {
4577 18 : BUN n = 0, i, cnt = BATcount(b1);
4578 18 : dbl mean1 = 0, mean2 = 0, up = 0, down1 = 0, down2 = 0, delta1, delta2, aux;
4579 18 : BATiter b1i = bat_iterator(b1), b2i = bat_iterator(b2);
4580 19 : const void *v1 = b1i.base, *v2 = b2i.base;
4581 19 : int tp = b1i.type;
4582 19 : lng t0 = 0;
4583 :
4584 19 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
4585 :
4586 19 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4587 :
4588 19 : switch (tp) {
4589 : case TYPE_bte:
4590 11 : AGGR_CORRELATION_SINGLE(bte);
4591 : break;
4592 : case TYPE_sht:
4593 0 : AGGR_CORRELATION_SINGLE(sht);
4594 : break;
4595 : case TYPE_int:
4596 95 : AGGR_CORRELATION_SINGLE(int);
4597 : break;
4598 : case TYPE_lng:
4599 0 : AGGR_CORRELATION_SINGLE(lng);
4600 : break;
4601 : #ifdef HAVE_HGE
4602 : case TYPE_hge:
4603 0 : AGGR_CORRELATION_SINGLE(hge);
4604 : break;
4605 : #endif
4606 : case TYPE_flt:
4607 0 : AGGR_CORRELATION_SINGLE(flt);
4608 : break;
4609 : case TYPE_dbl:
4610 17 : AGGR_CORRELATION_SINGLE(dbl);
4611 : break;
4612 0 : default:
4613 0 : GDKerror("type (%s) not supported.\n", ATOMname(tp));
4614 0 : goto bailout;
4615 : }
4616 18 : bat_iterator_end(&b1i);
4617 18 : bat_iterator_end(&b2i);
4618 18 : if (n != 0 && down1 != 0 && down2 != 0)
4619 9 : aux = (up / n) / (sqrt(down1 / n) * sqrt(down2 / n));
4620 : else
4621 9 : aux = dbl_nil;
4622 18 : TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT " (" LLFMT " usec)\n",
4623 : ALGOBATPAR(b1), ALGOBATPAR(b2), GDKusec() - t0);
4624 : return aux;
4625 1 : overflow:
4626 1 : GDKerror("22003!overflow in calculation.\n");
4627 1 : bailout:
4628 1 : bat_iterator_end(&b1i);
4629 1 : bat_iterator_end(&b2i);
4630 1 : return dbl_nil;
4631 : }
4632 :
4633 : #define AGGR_STDEV(TYPE) \
4634 : do { \
4635 : const TYPE *restrict vals = (const TYPE *) bi.base; \
4636 : TIMEOUT_LOOP(ci.ncand, qry_ctx) { \
4637 : i = canditer_next(&ci) - b->hseqbase; \
4638 : if (gids == NULL || \
4639 : (gids[i] >= min && gids[i] <= max)) { \
4640 : if (gids) \
4641 : gid = gids[i] - min; \
4642 : else \
4643 : gid = (oid) i; \
4644 : if (is_##TYPE##_nil(vals[i])) { \
4645 : if (!skip_nils) \
4646 : cnts[gid] = BUN_NONE; \
4647 : } else if (cnts[gid] != BUN_NONE) { \
4648 : cnts[gid]++; \
4649 : delta[gid] = (dbl) vals[i] - mean[gid]; \
4650 : mean[gid] += delta[gid] / cnts[gid]; \
4651 : m2[gid] += delta[gid] * ((dbl) vals[i] - mean[gid]); \
4652 : } \
4653 : } \
4654 : } \
4655 : TIMEOUT_CHECK(qry_ctx, \
4656 : GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx)); \
4657 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) { \
4658 : if (cnts[i] == 0 || cnts[i] == BUN_NONE) { \
4659 : dbls[i] = dbl_nil; \
4660 : mean[i] = dbl_nil; \
4661 : nils++; \
4662 : } else if (cnts[i] == 1) { \
4663 : dbls[i] = issample ? dbl_nil : 0; \
4664 : nils2++; \
4665 : } else if (isinf(m2[i])) { \
4666 : goto overflow; \
4667 : } else if (variance) { \
4668 : dbls[i] = m2[i] / (cnts[i] - issample); \
4669 : } else { \
4670 : dbls[i] = sqrt(m2[i] / (cnts[i] - issample)); \
4671 : } \
4672 : } \
4673 : } while (0)
4674 :
4675 : /* Calculate group standard deviation (population (i.e. biased) or
4676 : * sample (i.e. non-biased)) with optional candidates list.
4677 : *
4678 : * Note that this helper function is prepared to return two BATs: one
4679 : * (as return value) with the standard deviation per group, and one
4680 : * (as return argument) with the average per group. This isn't
4681 : * currently used since it doesn't fit into the mold of grouped
4682 : * aggregates. */
4683 : static BAT *
4684 18 : dogroupstdev(BAT **avgb, BAT *b, BAT *g, BAT *e, BAT *s, int tp,
4685 : bool skip_nils, bool issample, bool variance, const char *func)
4686 : {
4687 18 : const oid *restrict gids;
4688 18 : oid gid;
4689 18 : oid min, max;
4690 18 : BUN i, ngrp;
4691 18 : BUN nils = 0, nils2 = 0;
4692 18 : BUN *restrict cnts = NULL;
4693 18 : dbl *restrict dbls, *restrict mean, *restrict delta, *restrict m2;
4694 18 : BAT *bn = NULL, *an = NULL;
4695 18 : struct canditer ci;
4696 18 : const char *err;
4697 18 : lng t0 = 0;
4698 18 : BATiter bi = {0};
4699 :
4700 18 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
4701 :
4702 18 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4703 :
4704 18 : assert(tp == TYPE_dbl);
4705 18 : (void) tp; /* compatibility (with other BATgroup*
4706 : * functions) argument */
4707 :
4708 18 : if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
4709 0 : GDKerror("%s: %s\n", func, err);
4710 0 : return NULL;
4711 : }
4712 18 : if (g == NULL) {
4713 0 : GDKerror("%s: b and g must be aligned\n", func);
4714 0 : return NULL;
4715 : }
4716 :
4717 18 : if (BATcount(b) == 0 || ngrp == 0) {
4718 : /* trivial: no products, so return bat aligned with g
4719 : * with nil in the tail */
4720 3 : bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
4721 3 : goto doreturn;
4722 : }
4723 :
4724 15 : if ((e == NULL ||
4725 15 : (BATcount(e) == ci.ncand && e->hseqbase == ci.hseq)) &&
4726 7 : (BATtdense(g) || (g->tkey && g->tnonil)) &&
4727 4 : (issample || b->tnonil)) {
4728 : /* trivial: singleton groups, so all results are equal
4729 : * to zero (population) or nil (sample) */
4730 3 : dbl v = issample ? dbl_nil : 0;
4731 4 : bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &v, ngrp, TRANSIENT);
4732 3 : goto doreturn;
4733 : }
4734 :
4735 11 : delta = GDKmalloc(ngrp * sizeof(dbl));
4736 11 : m2 = GDKmalloc(ngrp * sizeof(dbl));
4737 11 : cnts = GDKzalloc(ngrp * sizeof(BUN));
4738 11 : if (avgb) {
4739 0 : an = COLnew(0, TYPE_dbl, ngrp, TRANSIENT);
4740 0 : *avgb = an;
4741 0 : if (an == NULL) {
4742 0 : mean = NULL;
4743 0 : goto alloc_fail;
4744 : }
4745 0 : mean = (dbl *) Tloc(an, 0);
4746 : } else {
4747 11 : mean = GDKmalloc(ngrp * sizeof(dbl));
4748 : }
4749 11 : if (mean == NULL || delta == NULL || m2 == NULL || cnts == NULL)
4750 0 : goto alloc_fail;
4751 :
4752 11 : bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
4753 11 : if (bn == NULL)
4754 0 : goto alloc_fail;
4755 11 : dbls = (dbl *) Tloc(bn, 0);
4756 :
4757 1080124 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) {
4758 1080038 : mean[i] = 0;
4759 1080038 : delta[i] = 0;
4760 1080038 : m2[i] = 0;
4761 : }
4762 :
4763 11 : bi = bat_iterator(b);
4764 11 : if (BATtdense(g))
4765 : gids = NULL;
4766 : else
4767 8 : gids = (const oid *) Tloc(g, 0);
4768 :
4769 11 : switch (b->ttype) {
4770 0 : case TYPE_bte:
4771 0 : AGGR_STDEV(bte);
4772 : break;
4773 0 : case TYPE_sht:
4774 0 : AGGR_STDEV(sht);
4775 : break;
4776 4 : case TYPE_int:
4777 5760376 : AGGR_STDEV(int);
4778 : break;
4779 0 : case TYPE_lng:
4780 0 : AGGR_STDEV(lng);
4781 : break;
4782 : #ifdef HAVE_HGE
4783 0 : case TYPE_hge:
4784 0 : AGGR_STDEV(hge);
4785 : break;
4786 : #endif
4787 0 : case TYPE_flt:
4788 0 : AGGR_STDEV(flt);
4789 : break;
4790 7 : case TYPE_dbl:
4791 104 : AGGR_STDEV(dbl);
4792 : break;
4793 0 : default:
4794 0 : GDKerror("%s: type (%s) not supported.\n",
4795 : func, ATOMname(b->ttype));
4796 0 : goto bailout;
4797 : }
4798 9 : bat_iterator_end(&bi);
4799 9 : if (an) {
4800 0 : BATsetcount(an, ngrp);
4801 0 : an->tkey = ngrp <= 1;
4802 0 : an->tsorted = ngrp <= 1;
4803 0 : an->trevsorted = ngrp <= 1;
4804 0 : an->tnil = nils != 0;
4805 0 : an->tnonil = nils == 0;
4806 : } else {
4807 9 : GDKfree(mean);
4808 : }
4809 9 : if (issample)
4810 3 : nils += nils2;
4811 9 : GDKfree(delta);
4812 9 : GDKfree(m2);
4813 9 : GDKfree(cnts);
4814 9 : BATsetcount(bn, ngrp);
4815 9 : bn->tkey = ngrp <= 1;
4816 9 : bn->tsorted = ngrp <= 1;
4817 9 : bn->trevsorted = ngrp <= 1;
4818 9 : bn->tnil = nils != 0;
4819 9 : bn->tnonil = nils == 0;
4820 15 : doreturn:
4821 15 : TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOBATFMT ",e=" ALGOOPTBATFMT
4822 : ",s=" ALGOOPTBATFMT
4823 : ",skip_nils=%s,issample=%s,variance=%s -> " ALGOOPTBATFMT
4824 : ",avgb=" ALGOOPTBATFMT " (%s -- " LLFMT " usec)\n",
4825 : ALGOBATPAR(b), ALGOBATPAR(g), ALGOOPTBATPAR(e),
4826 : ALGOOPTBATPAR(s),
4827 : skip_nils ? "true" : "false",
4828 : issample ? "true" : "false",
4829 : variance ? "true" : "false",
4830 : ALGOOPTBATPAR(bn), ALGOOPTBATPAR(an),
4831 : func, GDKusec() - t0);
4832 : return bn;
4833 2 : overflow:
4834 2 : GDKerror("22003!overflow in calculation.\n");
4835 2 : bailout:
4836 2 : bat_iterator_end(&bi);
4837 2 : alloc_fail:
4838 2 : if (an)
4839 0 : BBPreclaim(an);
4840 : else
4841 2 : GDKfree(mean);
4842 2 : BBPreclaim(bn);
4843 2 : GDKfree(delta);
4844 2 : GDKfree(m2);
4845 2 : GDKfree(cnts);
4846 2 : return NULL;
4847 : }
4848 :
4849 : BAT *
4850 7 : BATgroupstdev_sample(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
4851 : {
4852 7 : return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, true, false,
4853 : __func__);
4854 : }
4855 :
4856 : BAT *
4857 5 : BATgroupstdev_population(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
4858 : bool skip_nils)
4859 : {
4860 5 : return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, false, false,
4861 : __func__);
4862 : }
4863 :
4864 : BAT *
4865 1 : BATgroupvariance_sample(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
4866 : bool skip_nils)
4867 : {
4868 1 : return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, true, true,
4869 : __func__);
4870 : }
4871 :
4872 : BAT *
4873 5 : BATgroupvariance_population(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
4874 : bool skip_nils)
4875 : {
4876 5 : return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, false, true,
4877 : __func__);
4878 : }
4879 :
4880 : #define AGGR_COVARIANCE(TYPE) \
4881 : do { \
4882 : const TYPE *vals1 = (const TYPE *) b1i.base; \
4883 : const TYPE *vals2 = (const TYPE *) b2i.base; \
4884 : TIMEOUT_LOOP(ci.ncand, qry_ctx) { \
4885 : i = canditer_next(&ci) - b1->hseqbase; \
4886 : if (gids == NULL || \
4887 : (gids[i] >= min && gids[i] <= max)) { \
4888 : if (gids) \
4889 : gid = gids[i] - min; \
4890 : else \
4891 : gid = (oid) i; \
4892 : if (is_##TYPE##_nil(vals1[i]) || is_##TYPE##_nil(vals2[i])) { \
4893 : if (!skip_nils) \
4894 : cnts[gid] = BUN_NONE; \
4895 : } else if (cnts[gid] != BUN_NONE) { \
4896 : cnts[gid]++; \
4897 : delta1[gid] = (dbl) vals1[i] - mean1[gid]; \
4898 : mean1[gid] += delta1[gid] / cnts[gid]; \
4899 : delta2[gid] = (dbl) vals2[i] - mean2[gid]; \
4900 : mean2[gid] += delta2[gid] / cnts[gid]; \
4901 : m2[gid] += delta1[gid] * ((dbl) vals2[i] - mean2[gid]); \
4902 : } \
4903 : } \
4904 : } \
4905 : TIMEOUT_CHECK(qry_ctx, \
4906 : GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx)); \
4907 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) { \
4908 : if (cnts[i] == 0 || cnts[i] == BUN_NONE) { \
4909 : dbls[i] = dbl_nil; \
4910 : nils++; \
4911 : } else if (cnts[i] == 1) { \
4912 : dbls[i] = issample ? dbl_nil : 0; \
4913 : nils2++; \
4914 : } else if (isinf(m2[i])) { \
4915 : goto overflow; \
4916 : } else { \
4917 : dbls[i] = m2[i] / (cnts[i] - issample); \
4918 : } \
4919 : } \
4920 : } while (0)
4921 :
4922 : static BAT *
4923 60 : dogroupcovariance(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp,
4924 : bool skip_nils, bool issample, const char *func)
4925 : {
4926 60 : const oid *restrict gids;
4927 60 : oid gid, min, max;
4928 60 : BUN i, ngrp, nils = 0, nils2 = 0;
4929 60 : BUN *restrict cnts = NULL;
4930 60 : dbl *restrict dbls, *restrict mean1, *restrict mean2, *restrict delta1, *restrict delta2, *restrict m2;
4931 60 : BAT *bn = NULL;
4932 60 : struct canditer ci;
4933 60 : const char *err;
4934 60 : lng t0 = 0;
4935 60 : BATiter b1i = {0}, b2i = {0};
4936 :
4937 60 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
4938 :
4939 :
4940 60 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
4941 :
4942 60 : assert(tp == TYPE_dbl && BATcount(b1) == BATcount(b2) && b1->ttype == b2->ttype && BATtdense(b1) == BATtdense(b2));
4943 60 : (void) tp;
4944 :
4945 60 : if ((err = BATgroupaggrinit(b1, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
4946 0 : GDKerror("%s: %s\n", func, err);
4947 0 : return NULL;
4948 : }
4949 60 : if (g == NULL) {
4950 0 : GDKerror("%s: b1, b2 and g must be aligned\n", func);
4951 0 : return NULL;
4952 : }
4953 :
4954 60 : if (BATcount(b1) == 0 || ngrp == 0) {
4955 0 : bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
4956 0 : goto doreturn;
4957 : }
4958 :
4959 60 : if ((e == NULL ||
4960 60 : (BATcount(e) == BATcount(b1) && (e->hseqbase == b1->hseqbase || e->hseqbase == b2->hseqbase))) &&
4961 20 : (BATtdense(g) || (g->tkey && g->tnonil)) &&
4962 10 : (issample || (b1->tnonil && b2->tnonil))) {
4963 : /* trivial: singleton groups, so all results are equal
4964 : * to zero (population) or nil (sample) */
4965 10 : dbl v = issample ? dbl_nil : 0;
4966 11 : bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &v, ngrp, TRANSIENT);
4967 11 : goto doreturn;
4968 : }
4969 :
4970 49 : delta1 = GDKmalloc(ngrp * sizeof(dbl));
4971 49 : delta2 = GDKmalloc(ngrp * sizeof(dbl));
4972 49 : m2 = GDKmalloc(ngrp * sizeof(dbl));
4973 49 : cnts = GDKzalloc(ngrp * sizeof(BUN));
4974 48 : mean1 = GDKmalloc(ngrp * sizeof(dbl));
4975 49 : mean2 = GDKmalloc(ngrp * sizeof(dbl));
4976 :
4977 49 : if (mean1 == NULL || mean2 == NULL || delta1 == NULL || delta2 == NULL || m2 == NULL || cnts == NULL)
4978 0 : goto alloc_fail;
4979 :
4980 415 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) {
4981 318 : m2[i] = 0;
4982 318 : mean1[i] = 0;
4983 318 : mean2[i] = 0;
4984 : }
4985 :
4986 49 : bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
4987 48 : if (bn == NULL)
4988 0 : goto alloc_fail;
4989 48 : dbls = (dbl *) Tloc(bn, 0);
4990 :
4991 48 : if (!g || BATtdense(g))
4992 : gids = NULL;
4993 : else
4994 48 : gids = (const oid *) Tloc(g, 0);
4995 :
4996 48 : b1i = bat_iterator(b1);
4997 48 : b2i = bat_iterator(b2);
4998 46 : switch (b1->ttype) {
4999 9 : case TYPE_bte:
5000 184 : AGGR_COVARIANCE(bte);
5001 : break;
5002 0 : case TYPE_sht:
5003 0 : AGGR_COVARIANCE(sht);
5004 : break;
5005 37 : case TYPE_int:
5006 794 : AGGR_COVARIANCE(int);
5007 : break;
5008 0 : case TYPE_lng:
5009 0 : AGGR_COVARIANCE(lng);
5010 : break;
5011 : #ifdef HAVE_HGE
5012 0 : case TYPE_hge:
5013 0 : AGGR_COVARIANCE(hge);
5014 : break;
5015 : #endif
5016 0 : case TYPE_flt:
5017 0 : AGGR_COVARIANCE(flt);
5018 : break;
5019 0 : case TYPE_dbl:
5020 0 : AGGR_COVARIANCE(dbl);
5021 : break;
5022 0 : default:
5023 0 : GDKerror("%s: type (%s) not supported.\n", func, ATOMname(b1->ttype));
5024 0 : goto bailout;
5025 : }
5026 48 : bat_iterator_end(&b1i);
5027 46 : bat_iterator_end(&b2i);
5028 48 : GDKfree(mean1);
5029 48 : GDKfree(mean2);
5030 :
5031 48 : if (issample)
5032 20 : nils += nils2;
5033 48 : GDKfree(delta1);
5034 48 : GDKfree(delta2);
5035 48 : GDKfree(m2);
5036 48 : GDKfree(cnts);
5037 49 : BATsetcount(bn, ngrp);
5038 49 : bn->tkey = ngrp <= 1;
5039 49 : bn->tsorted = ngrp <= 1;
5040 49 : bn->trevsorted = ngrp <= 1;
5041 49 : bn->tnil = nils != 0;
5042 49 : bn->tnonil = nils == 0;
5043 60 : doreturn:
5044 60 : TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT ",g=" ALGOBATFMT
5045 : ",e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT
5046 : ",skip_nils=%s,issample=%s -> " ALGOOPTBATFMT
5047 : " (%s -- " LLFMT " usec)\n",
5048 : ALGOBATPAR(b1), ALGOBATPAR(b2), ALGOBATPAR(g),
5049 : ALGOOPTBATPAR(e), ALGOOPTBATPAR(s),
5050 : skip_nils ? "true" : "false",
5051 : issample ? "true" : "false",
5052 : ALGOOPTBATPAR(bn),
5053 : func, GDKusec() - t0);
5054 : return bn;
5055 0 : overflow:
5056 0 : GDKerror("22003!overflow in calculation.\n");
5057 0 : bailout:
5058 0 : bat_iterator_end(&b1i);
5059 0 : bat_iterator_end(&b2i);
5060 0 : alloc_fail:
5061 0 : BBPreclaim(bn);
5062 0 : GDKfree(mean1);
5063 0 : GDKfree(mean2);
5064 0 : GDKfree(delta1);
5065 0 : GDKfree(delta2);
5066 0 : GDKfree(m2);
5067 0 : GDKfree(cnts);
5068 0 : return NULL;
5069 : }
5070 :
5071 : BAT *
5072 30 : BATgroupcovariance_sample(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
5073 : {
5074 30 : return dogroupcovariance(b1, b2, g, e, s, tp, skip_nils, true,
5075 : __func__);
5076 : }
5077 :
5078 : BAT *
5079 30 : BATgroupcovariance_population(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
5080 : {
5081 30 : return dogroupcovariance(b1, b2, g, e, s, tp, skip_nils, false,
5082 : __func__);
5083 : }
5084 :
5085 : #define AGGR_CORRELATION(TYPE) \
5086 : do { \
5087 : const TYPE *vals1 = (const TYPE *) b1i.base; \
5088 : const TYPE *vals2 = (const TYPE *) b2i.base; \
5089 : TIMEOUT_LOOP(ci.ncand, qry_ctx) { \
5090 : i = canditer_next(&ci) - b1->hseqbase; \
5091 : if (gids == NULL || \
5092 : (gids[i] >= min && gids[i] <= max)) { \
5093 : if (gids) \
5094 : gid = gids[i] - min; \
5095 : else \
5096 : gid = (oid) i; \
5097 : if (is_##TYPE##_nil(vals1[i]) || is_##TYPE##_nil(vals2[i])) { \
5098 : if (!skip_nils) \
5099 : cnts[gid] = BUN_NONE; \
5100 : } else if (cnts[gid] != BUN_NONE) { \
5101 : cnts[gid]++; \
5102 : delta1[gid] = (dbl) vals1[i] - mean1[gid]; \
5103 : mean1[gid] += delta1[gid] / cnts[gid]; \
5104 : delta2[gid] = (dbl) vals2[i] - mean2[gid]; \
5105 : mean2[gid] += delta2[gid] / cnts[gid]; \
5106 : aux = (dbl) vals2[i] - mean2[gid]; \
5107 : up[gid] += delta1[gid] * aux; \
5108 : down1[gid] += delta1[gid] * ((dbl) vals1[i] - mean1[gid]); \
5109 : down2[gid] += delta2[gid] * aux; \
5110 : } \
5111 : } \
5112 : } \
5113 : TIMEOUT_CHECK(qry_ctx, \
5114 : GOTO_LABEL_TIMEOUT_HANDLER(bailout, qry_ctx)); \
5115 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) { \
5116 : if (cnts[i] <= 1 || cnts[i] == BUN_NONE || down1[i] == 0 || down2[i] == 0) { \
5117 : dbls[i] = dbl_nil; \
5118 : nils++; \
5119 : } else if (isinf(up[i]) || isinf(down1[i]) || isinf(down2[i])) { \
5120 : goto overflow; \
5121 : } else { \
5122 : dbls[i] = (up[i] / cnts[i]) / (sqrt(down1[i] / cnts[i]) * sqrt(down2[i] / cnts[i])); \
5123 : assert(!is_dbl_nil(dbls[i])); \
5124 : } \
5125 : } \
5126 : } while (0)
5127 :
5128 : BAT *
5129 49 : BATgroupcorrelation(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils)
5130 : {
5131 49 : const oid *restrict gids;
5132 49 : oid gid, min, max;
5133 49 : BUN i, ngrp, nils = 0;
5134 49 : BUN *restrict cnts = NULL;
5135 49 : dbl *restrict dbls, *restrict mean1, *restrict mean2, *restrict delta1, *restrict delta2, *restrict up, *restrict down1, *restrict down2, aux;
5136 49 : BAT *bn = NULL;
5137 49 : struct canditer ci;
5138 49 : const char *err;
5139 49 : lng t0 = 0;
5140 49 : BATiter b1i = {0}, b2i = {0};
5141 :
5142 49 : QryCtx *qry_ctx = MT_thread_get_qry_ctx();
5143 :
5144 49 : TRC_DEBUG_IF(ALGO) t0 = GDKusec();
5145 :
5146 49 : assert(tp == TYPE_dbl && BATcount(b1) == BATcount(b2) && b1->ttype == b2->ttype && BATtdense(b1) == BATtdense(b2));
5147 49 : (void) tp;
5148 :
5149 49 : if ((err = BATgroupaggrinit(b1, g, e, s, &min, &max, &ngrp, &ci)) != NULL) {
5150 0 : GDKerror("%s\n", err);
5151 0 : return NULL;
5152 : }
5153 49 : if (g == NULL) {
5154 0 : GDKerror("b1, b2 and g must be aligned\n");
5155 0 : return NULL;
5156 : }
5157 :
5158 49 : if (BATcount(b1) == 0 || ngrp == 0) {
5159 1 : bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
5160 1 : goto doreturn;
5161 : }
5162 :
5163 48 : if ((e == NULL ||
5164 48 : (BATcount(e) == BATcount(b1) && (e->hseqbase == b1->hseqbase || e->hseqbase == b2->hseqbase))) &&
5165 14 : (BATtdense(g) || (g->tkey && g->tnonil))) {
5166 14 : dbl v = dbl_nil;
5167 14 : bn = BATconstant(min, TYPE_dbl, &v, ngrp, TRANSIENT);
5168 14 : goto doreturn;
5169 : }
5170 :
5171 34 : delta1 = GDKmalloc(ngrp * sizeof(dbl));
5172 34 : delta2 = GDKmalloc(ngrp * sizeof(dbl));
5173 34 : up = GDKmalloc(ngrp * sizeof(dbl));
5174 34 : down1 = GDKmalloc(ngrp * sizeof(dbl));
5175 34 : down2 = GDKmalloc(ngrp * sizeof(dbl));
5176 34 : cnts = GDKzalloc(ngrp * sizeof(BUN));
5177 34 : mean1 = GDKmalloc(ngrp * sizeof(dbl));
5178 34 : mean2 = GDKmalloc(ngrp * sizeof(dbl));
5179 :
5180 34 : if (mean1 == NULL || mean2 == NULL || delta1 == NULL || delta2 == NULL || up == NULL || down1 == NULL || down2 == NULL || cnts == NULL)
5181 0 : goto alloc_fail;
5182 :
5183 265 : TIMEOUT_LOOP_IDX(i, ngrp, qry_ctx) {
5184 197 : up[i] = 0;
5185 197 : down1[i] = 0;
5186 197 : down2[i] = 0;
5187 197 : mean1[i] = 0;
5188 197 : mean2[i] = 0;
5189 : }
5190 :
5191 34 : bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
5192 34 : if (bn == NULL)
5193 0 : goto alloc_fail;
5194 34 : dbls = (dbl *) Tloc(bn, 0);
5195 :
5196 34 : if (!g || BATtdense(g))
5197 : gids = NULL;
5198 : else
5199 34 : gids = (const oid *) Tloc(g, 0);
5200 :
5201 34 : b1i = bat_iterator(b1);
5202 34 : b2i = bat_iterator(b2);
5203 34 : switch (b1->ttype) {
5204 4 : case TYPE_bte:
5205 80 : AGGR_CORRELATION(bte);
5206 : break;
5207 0 : case TYPE_sht:
5208 0 : AGGR_CORRELATION(sht);
5209 : break;
5210 29 : case TYPE_int:
5211 1235 : AGGR_CORRELATION(int);
5212 : break;
5213 0 : case TYPE_lng:
5214 0 : AGGR_CORRELATION(lng);
5215 : break;
5216 : #ifdef HAVE_HGE
5217 0 : case TYPE_hge:
5218 0 : AGGR_CORRELATION(hge);
5219 : break;
5220 : #endif
5221 0 : case TYPE_flt:
5222 0 : AGGR_CORRELATION(flt);
5223 : break;
5224 1 : case TYPE_dbl:
5225 4 : AGGR_CORRELATION(dbl);
5226 : break;
5227 0 : default:
5228 0 : GDKerror("type (%s) not supported.\n", ATOMname(b1->ttype));
5229 0 : goto bailout;
5230 : }
5231 33 : bat_iterator_end(&b1i);
5232 33 : bat_iterator_end(&b2i);
5233 33 : GDKfree(mean1);
5234 33 : GDKfree(mean2);
5235 33 : GDKfree(delta1);
5236 33 : GDKfree(delta2);
5237 33 : GDKfree(up);
5238 33 : GDKfree(down1);
5239 33 : GDKfree(down2);
5240 32 : GDKfree(cnts);
5241 33 : BATsetcount(bn, ngrp);
5242 33 : bn->tkey = ngrp <= 1;
5243 33 : bn->tsorted = ngrp <= 1;
5244 33 : bn->trevsorted = ngrp <= 1;
5245 33 : bn->tnil = nils != 0;
5246 33 : bn->tnonil = nils == 0;
5247 48 : doreturn:
5248 48 : TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT ",g=" ALGOBATFMT
5249 : ",e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT
5250 : ",skip_nils=%s -> " ALGOOPTBATFMT
5251 : " (" LLFMT " usec)\n",
5252 : ALGOBATPAR(b1), ALGOBATPAR(b2), ALGOBATPAR(g),
5253 : ALGOOPTBATPAR(e), ALGOOPTBATPAR(s),
5254 : skip_nils ? "true" : "false",
5255 : ALGOOPTBATPAR(bn),
5256 : GDKusec() - t0);
5257 : return bn;
5258 1 : overflow:
5259 1 : GDKerror("22003!overflow in calculation.\n");
5260 1 : bailout:
5261 1 : bat_iterator_end(&b1i);
5262 1 : bat_iterator_end(&b2i);
5263 1 : alloc_fail:
5264 1 : BBPreclaim(bn);
5265 1 : GDKfree(mean1);
5266 1 : GDKfree(mean2);
5267 1 : GDKfree(delta1);
5268 1 : GDKfree(delta2);
5269 1 : GDKfree(up);
5270 1 : GDKfree(down1);
5271 1 : GDKfree(down2);
5272 1 : GDKfree(cnts);
5273 1 : return NULL;
5274 : }
|
