1 #include <aggregate/reducer.h>
2 #include <math.h>
3
4 typedef struct {
5 const RLookupKey *srckey;
6 size_t n;
7 double oldM, newM, oldS, newS;
8 } devCtx;
9
10 #define BLOCK_SIZE 1024 * sizeof(devCtx)
11
stddevNewInstance(Reducer * rbase)12 static void *stddevNewInstance(Reducer *rbase) {
13 devCtx *dctx = BlkAlloc_Alloc(&rbase->alloc, sizeof(*dctx), BLOCK_SIZE);
14 memset(dctx, 0, sizeof(*dctx));
15 dctx->srckey = rbase->srckey;
16 return dctx;
17 }
18
stddevAddInternal(devCtx * dctx,double d)19 static void stddevAddInternal(devCtx *dctx, double d) {
20 // https://www.johndcook.com/blog/standard_deviation/
21 dctx->n++;
22 if (dctx->n == 1) {
23 dctx->oldM = dctx->newM = d;
24 dctx->oldS = 0.0;
25 } else {
26 dctx->newM = dctx->oldM + (d - dctx->oldM) / dctx->n;
27 dctx->newS = dctx->oldS + (d - dctx->oldM) * (d - dctx->newM);
28
29 // set up for next iteration
30 dctx->oldM = dctx->newM;
31 dctx->oldS = dctx->newS;
32 }
33 }
34
stddevAdd(Reducer * r,void * ctx,const RLookupRow * srcrow)35 static int stddevAdd(Reducer *r, void *ctx, const RLookupRow *srcrow) {
36 devCtx *dctx = ctx;
37 double d;
38 RSValue *v = RLookup_GetItem(dctx->srckey, srcrow);
39 if (v) {
40 if (v->t != RSValue_Array) {
41 if (RSValue_ToNumber(v, &d)) {
42 stddevAddInternal(dctx, d);
43 }
44 } else {
45 uint32_t sz = RSValue_ArrayLen(v);
46 for (uint32_t i = 0; i < sz; i++) {
47 if (RSValue_ToNumber(RSValue_ArrayItem(v, i), &d)) {
48 stddevAddInternal(dctx, d);
49 }
50 }
51 }
52 }
53 return 1;
54 }
55
stddevFinalize(Reducer * parent,void * instance)56 static RSValue *stddevFinalize(Reducer *parent, void *instance) {
57 devCtx *dctx = instance;
58 double variance = ((dctx->n > 1) ? dctx->newS / (dctx->n - 1) : 0.0);
59 double stddev = sqrt(variance);
60 return RS_NumVal(stddev);
61 }
62
RDCRStdDev_New(const ReducerOptions * options)63 Reducer *RDCRStdDev_New(const ReducerOptions *options) {
64 Reducer *r = rm_calloc(1, sizeof(*r));
65 if (!ReducerOptions_GetKey(options, &r->srckey)) {
66 rm_free(r);
67 return NULL;
68 }
69 r->Add = stddevAdd;
70 r->Finalize = stddevFinalize;
71 r->Free = Reducer_GenericFree;
72 r->NewInstance = stddevNewInstance;
73 r->reducerId = REDUCER_T_STDDEV;
74 return r;
75 }
76