/** Histogram functions. * * @author Steffen Vogel * @copyright 2017, Institute for Automation of Complex Power Systems, EONERC * @license GNU General Public License (version 3) * * VILLASnode * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program. If not, see . *********************************************************************************/ #include #include #include #include #include #include #include #include #include #include #define VAL(h, i) ((h)->low + (i) * (h)->resolution) #define INDEX(h, v) round((v - (h)->low) / (h)->resolution) int hist_init(struct hist *h, int buckets, hist_cnt_t warmup) { h->length = buckets; h->warmup = warmup; h->data = buckets ? alloc(h->length * sizeof(hist_cnt_t)) : NULL; hist_reset(h); return 0; } int hist_destroy(struct hist *h) { if (h->data) { free(h->data); h->data = NULL; } return 0; } void hist_put(struct hist *h, double value) { h->last = value; /* Update min/max */ if (value > h->highest) h->highest = value; if (value < h->lowest) h->lowest = value; if (h->total < h->warmup) { } else if (h->total == h->warmup) { h->low = hist_mean(h) - 3 * hist_stddev(h); h->high = hist_mean(h) + 3 * hist_stddev(h); h->resolution = (h->high - h->low) / h->length; } else { int idx = INDEX(h, value); /* Check bounds and increment */ if (idx >= h->length) h->higher++; else if (idx < 0) h->lower++; else if (h->data != NULL) h->data[idx]++; } h->total++; /* Online / running calculation of variance and mean * by Donald Knuth’s Art of Computer Programming, Vol 2, page 232, 3rd edition */ if (h->total == 1) { h->_m[1] = h->_m[0] = value; h->_s[1] = 0.0; } else { h->_m[0] = h->_m[1] + (value - h->_m[1]) / h->total; h->_s[0] = h->_s[1] + (value - h->_m[1]) * (value - h->_m[0]); // set up for next iteration h->_m[1] = h->_m[0]; h->_s[1] = h->_s[0]; } } void hist_reset(struct hist *h) { h->total = 0; h->higher = 0; h->lower = 0; h->highest = DBL_MIN; h->lowest = DBL_MAX; if (h->data) memset(h->data, 0, h->length * sizeof(unsigned)); } double hist_mean(struct hist *h) { return (h->total > 0) ? h->_m[0] : NAN; } double hist_var(struct hist *h) { return (h->total > 1) ? h->_s[0] / (h->total - 1) : NAN; } double hist_stddev(struct hist *h) { return sqrt(hist_var(h)); } void hist_print(struct hist *h, int details) { if (h->total > 0) { hist_cnt_t missed = h->total - h->higher - h->lower; stats("Counted values: %ju (%ju between %f and %f)", h->total, missed, h->low, h->high); stats("Highest: %g", h->highest); stats("Lowest: %g", h->lowest); stats("Mu: %g", hist_mean(h)); stats("1/Mu: %g", 1.0/hist_mean(h)); stats("Variance: %g", hist_var(h)); stats("Stddev: %g", hist_stddev(h)); if (details > 0 && h->total - h->higher - h->lower > 0) { char *buf = hist_dump(h); stats("Matlab: %s", buf); free(buf); hist_plot(h); } } else stats("Counted values: %ju", h->total); } void hist_plot(struct hist *h) { hist_cnt_t max = 1; /* Get highest bar */ for (int i = 0; i < h->length; i++) { if (h->data[i] > max) max = h->data[i]; } struct table_column cols[] = { { -9, "Value", "%+9.3g", NULL, TABLE_ALIGN_RIGHT }, { -6, "Count", "%6ju", NULL, TABLE_ALIGN_RIGHT }, { 0, "Plot", "%s", "occurences", TABLE_ALIGN_LEFT } }; struct table table = { .ncols = ARRAY_LEN(cols), .cols = cols }; /* Print plot */ table_header(&table); for (int i = 0; i < h->length; i++) { double value = VAL(h, i); hist_cnt_t cnt = h->data[i]; int bar = cols[2]._width * ((double) cnt / max); char *buf = strf("%s", ""); for (int i = 0; i < bar; i++) buf = strcatf(&buf, "\u2588"); table_row(&table, value, cnt, buf); free(buf); } table_footer(&table); } char * hist_dump(struct hist *h) { char *buf = alloc(128); strcatf(&buf, "[ "); for (int i = 0; i < h->length; i++) strcatf(&buf, "%ju ", h->data[i]); strcatf(&buf, "]"); return buf; } json_t * hist_json(struct hist *h) { json_t *json_buckets, *json_hist; json_hist = json_pack("{ s: f, s: f, s: i }", "low", h->low, "high", h->high, "total", h->total ); if (h->total > 0) { json_object_update(json_hist, json_pack("{ s: i, s: i, s: f, s: f, s: f, s: f, s: f }", "higher", h->higher, "lower", h->lower, "highest", h->highest, "lowest", h->lowest, "mean", hist_mean(h), "variance", hist_var(h), "stddev", hist_stddev(h) )); } if (h->total - h->lower - h->higher > 0) { json_buckets = json_array(); for (int i = 0; i < h->length; i++) json_array_append(json_buckets, json_integer(h->data[i])); json_object_set(json_hist, "buckets", json_buckets); } return json_hist; } int hist_dump_json(struct hist *h, FILE *f) { json_t *j = hist_json(h); int ret = json_dumpf(j, f, 0); json_decref(j); return ret; } int hist_dump_matlab(struct hist *h, FILE *f) { fprintf(f, "struct("); fprintf(f, "'low', %f, ", h->low); fprintf(f, "'high', %f, ", h->high); fprintf(f, "'total', %ju, ", h->total); fprintf(f, "'higher', %ju, ", h->higher); fprintf(f, "'lower', %ju, ", h->lower); fprintf(f, "'highest', %f, ", h->highest); fprintf(f, "'lowest', %f, ", h->lowest); fprintf(f, "'mean', %f, ", hist_mean(h)); fprintf(f, "'variance', %f, ", hist_var(h)); fprintf(f, "'stddev', %f, ", hist_stddev(h)); if (h->total - h->lower - h->higher > 0) { char *buf = hist_dump(h); fprintf(f, "'buckets', %s", buf); free(buf); } else fprintf(f, "'buckets', zeros(1, %d)", h->length); fprintf(f, ")"); return 0; }