Actual source code: hists.c
2: /*
3: Contains the data structure for plotting a histogram in a window with an axis.
4: */
5: #include <petscdraw.h>
6: #include <petsc/private/petscimpl.h>
7: #include <petscviewer.h>
9: PetscClassId PETSC_DRAWHG_CLASSID = 0;
11: struct _p_PetscDrawHG {
12: PETSCHEADER(int);
13: PetscErrorCode (*destroy)(PetscDrawSP);
14: PetscErrorCode (*view)(PetscDrawSP, PetscViewer);
15: PetscDraw win;
16: PetscDrawAxis axis;
17: PetscReal xmin, xmax;
18: PetscReal ymin, ymax;
19: int numBins;
20: int maxBins;
21: PetscReal *bins;
22: int numValues;
23: int maxValues;
24: PetscReal *values;
25: int color;
26: PetscBool calcStats;
27: PetscBool integerBins;
28: };
30: #define CHUNKSIZE 100
32: /*@C
33: PetscDrawHGCreate - Creates a histogram data structure.
35: Collective
37: Input Parameters:
38: + draw - The window where the graph will be made
39: - bins - The number of bins to use
41: Output Parameters:
42: . hist - The histogram context
44: Notes:
45: The difference between a bar chart, `PetscDrawBar`, and a histogram, `PetscDrawHG`, is explained here https://stattrek.com/statistics/charts/histogram.aspx?Tutorial=AP
47: The histogram is only displayed when `PetscDrawHGDraw()` is called.
49: The MPI communicator that owns the `PetscDraw` owns this `PetscDrawHG`, but the calls to set options and add data are ignored on all processes except the
50: zeroth MPI process in the communicator. All MPI ranks in the communicator must call `PetscDrawHGDraw()` to display the updated graph.
52: Level: intermediate
54: .seealso: `PetscDrawHGDestroy()`, `PetscDrawHG`, `PetscDrawBarCreate()`, `PetscDrawBar`, `PetscDrawLGCreate()`, `PetscDrawLG`, `PetscDrawSPCreate()`, `PetscDrawSP`,
55: `PetscDrawHGSetNumberBins()`, `PetscDrawHGReset()`, `PetscDrawHGAddValue()`, `PetscDrawHGDraw()`, `PetscDrawHGSave()`, `PetscDrawHGView()`, `PetscDrawHGSetColor()`,
56: `PetscDrawHGSetLimits()`, `PetscDrawHGCalcStats()`, `PetscDrawHGIntegerBins()`, `PetscDrawHGGetAxis()`, `PetscDrawAxis`, `PetscDrawHGGetDraw()`
57: @*/
58: PetscErrorCode PetscDrawHGCreate(PetscDraw draw, int bins, PetscDrawHG *hist)
59: {
60: PetscDrawHG h;
66: PetscHeaderCreate(h, PETSC_DRAWHG_CLASSID, "DrawHG", "Histogram", "Draw", PetscObjectComm((PetscObject)draw), PetscDrawHGDestroy, NULL);
68: PetscObjectReference((PetscObject)draw);
69: h->win = draw;
71: h->view = NULL;
72: h->destroy = NULL;
73: h->color = PETSC_DRAW_GREEN;
74: h->xmin = PETSC_MAX_REAL;
75: h->xmax = PETSC_MIN_REAL;
76: h->ymin = 0.;
77: h->ymax = 1.;
78: h->numBins = bins;
79: h->maxBins = bins;
81: PetscMalloc1(h->maxBins, &h->bins);
83: h->numValues = 0;
84: h->maxValues = CHUNKSIZE;
85: h->calcStats = PETSC_FALSE;
86: h->integerBins = PETSC_FALSE;
88: PetscMalloc1(h->maxValues, &h->values);
89: PetscDrawAxisCreate(draw, &h->axis);
91: *hist = h;
92: return 0;
93: }
95: /*@
96: PetscDrawHGSetNumberBins - Change the number of bins that are to be drawn in the histogram
98: Logically Collective
100: Input Parameters:
101: + hist - The histogram context.
102: - bins - The number of bins.
104: Level: intermediate
106: .seealso: `PetscDrawHGCreate()`, `PetscDrawHG`, `PetscDrawHGDraw()`, `PetscDrawHGIntegerBins()`
107: @*/
108: PetscErrorCode PetscDrawHGSetNumberBins(PetscDrawHG hist, int bins)
109: {
113: if (hist->maxBins < bins) {
114: PetscFree(hist->bins);
115: PetscMalloc1(bins, &hist->bins);
116: hist->maxBins = bins;
117: }
118: hist->numBins = bins;
119: return 0;
120: }
122: /*@
123: PetscDrawHGReset - Clears histogram to allow for reuse with new data.
125: Logically Collective
127: Input Parameter:
128: . hist - The histogram context.
130: Level: intermediate
132: .seealso: `PetscDrawHGCreate()`, `PetscDrawHG`, `PetscDrawHGDraw()`, `PetscDrawHGAddValue()`
133: @*/
134: PetscErrorCode PetscDrawHGReset(PetscDrawHG hist)
135: {
138: hist->xmin = PETSC_MAX_REAL;
139: hist->xmax = PETSC_MIN_REAL;
140: hist->ymin = 0.0;
141: hist->ymax = 0.0;
142: hist->numValues = 0;
143: return 0;
144: }
146: /*@C
147: PetscDrawHGDestroy - Frees all space taken up by histogram data structure.
149: Collective
151: Input Parameter:
152: . hist - The histogram context
154: Level: intermediate
156: .seealso: `PetscDrawHGCreate()`, `PetscDrawHG`
157: @*/
158: PetscErrorCode PetscDrawHGDestroy(PetscDrawHG *hist)
159: {
160: if (!*hist) return 0;
162: if (--((PetscObject)(*hist))->refct > 0) {
163: *hist = NULL;
164: return 0;
165: }
167: PetscFree((*hist)->bins);
168: PetscFree((*hist)->values);
169: PetscDrawAxisDestroy(&(*hist)->axis);
170: PetscDrawDestroy(&(*hist)->win);
171: PetscHeaderDestroy(hist);
172: return 0;
173: }
175: /*@
176: PetscDrawHGAddValue - Adds another value to the histogram.
178: Logically Collective
180: Input Parameters:
181: + hist - The histogram
182: - value - The value
184: Level: intermediate
186: .seealso: `PetscDrawHGCreate()`, `PetscDrawHG`, `PetscDrawHGDraw()`, `PetscDrawHGAddValue()`, `PetscDrawHGReset()`
187: @*/
188: PetscErrorCode PetscDrawHGAddValue(PetscDrawHG hist, PetscReal value)
189: {
192: /* Allocate more memory if necessary */
193: if (hist->numValues >= hist->maxValues) {
194: PetscReal *tmp;
196: PetscMalloc1(hist->maxValues + CHUNKSIZE, &tmp);
197: PetscArraycpy(tmp, hist->values, hist->maxValues);
198: PetscFree(hist->values);
200: hist->values = tmp;
201: hist->maxValues += CHUNKSIZE;
202: }
203: /* I disagree with the original Petsc implementation here. There should be no overshoot, but rather the
204: stated convention of using half-open intervals (always the way to go) */
205: if (!hist->numValues && (hist->xmin == PETSC_MAX_REAL) && (hist->xmax == PETSC_MIN_REAL)) {
206: hist->xmin = value;
207: hist->xmax = value;
208: #if 1
209: } else {
210: /* Update limits */
211: if (value > hist->xmax) hist->xmax = value;
212: if (value < hist->xmin) hist->xmin = value;
213: #else
214: } else if (hist->numValues == 1) {
215: /* Update limits -- We need to overshoot the largest value somewhat */
216: if (value > hist->xmax) hist->xmax = value + 0.001 * (value - hist->xmin) / hist->numBins;
217: if (value < hist->xmin) {
218: hist->xmin = value;
219: hist->xmax = hist->xmax + 0.001 * (hist->xmax - hist->xmin) / hist->numBins;
220: }
221: } else {
222: /* Update limits -- We need to overshoot the largest value somewhat */
223: if (value > hist->xmax) hist->xmax = value + 0.001 * (hist->xmax - hist->xmin) / hist->numBins;
224: if (value < hist->xmin) hist->xmin = value;
225: #endif
226: }
228: hist->values[hist->numValues++] = value;
229: return 0;
230: }
232: /*@
233: PetscDrawHGDraw - Redraws a histogram.
235: Collective
237: Input Parameter:
238: . hist - The histogram context
240: Level: intermediate
242: .seealso: `PetscDrawHGCreate()`, `PetscDrawHG`, `PetscDrawHGDraw()`, `PetscDrawHGAddValue()`, `PetscDrawHGReset()`
243: @*/
244: PetscErrorCode PetscDrawHGDraw(PetscDrawHG hist)
245: {
246: PetscDraw draw;
247: PetscBool isnull;
248: PetscReal xmin, xmax, ymin, ymax, *bins, *values, binSize, binLeft, binRight, maxHeight, mean, var;
249: char title[256];
250: char xlabel[256];
251: PetscInt numBins, numBinsOld, numValues, initSize, i, p, bcolor, color;
252: PetscMPIInt rank;
255: PetscDrawIsNull(hist->win, &isnull);
256: if (isnull) return 0;
257: MPI_Comm_rank(PetscObjectComm((PetscObject)hist), &rank);
259: if ((hist->xmin >= hist->xmax) || (hist->ymin >= hist->ymax)) return 0;
260: if (hist->numValues < 1) return 0;
262: color = hist->color;
263: if (color == PETSC_DRAW_ROTATE) bcolor = PETSC_DRAW_BLACK + 1;
264: else bcolor = color;
266: xmin = hist->xmin;
267: xmax = hist->xmax;
268: ymin = hist->ymin;
269: ymax = hist->ymax;
270: numValues = hist->numValues;
271: values = hist->values;
272: mean = 0.0;
273: var = 0.0;
275: draw = hist->win;
276: PetscDrawCheckResizedWindow(draw);
277: PetscDrawClear(draw);
279: if (xmin == xmax) {
280: /* Calculate number of points in each bin */
281: bins = hist->bins;
282: bins[0] = 0.;
283: for (p = 0; p < numValues; p++) {
284: if (values[p] == xmin) bins[0]++;
285: mean += values[p];
286: var += values[p] * values[p];
287: }
288: maxHeight = bins[0];
289: if (maxHeight > ymax) ymax = hist->ymax = maxHeight;
290: xmax = xmin + 1;
291: PetscDrawAxisSetLimits(hist->axis, xmin, xmax, ymin, ymax);
292: if (hist->calcStats) {
293: mean /= numValues;
294: if (numValues > 1) var = (var - numValues * mean * mean) / (numValues - 1);
295: else var = 0.0;
296: PetscSNPrintf(title, 256, "Mean: %g Var: %g", (double)mean, (double)var);
297: PetscSNPrintf(xlabel, 256, "Total: %" PetscInt_FMT, numValues);
298: PetscDrawAxisSetLabels(hist->axis, title, xlabel, NULL);
299: }
300: PetscDrawAxisDraw(hist->axis);
301: PetscDrawCollectiveBegin(draw);
302: if (rank == 0) { /* Draw bins */
303: binLeft = xmin;
304: binRight = xmax;
305: PetscDrawRectangle(draw, binLeft, ymin, binRight, bins[0], bcolor, bcolor, bcolor, bcolor);
306: PetscDrawLine(draw, binLeft, ymin, binLeft, bins[0], PETSC_DRAW_BLACK);
307: PetscDrawLine(draw, binRight, ymin, binRight, bins[0], PETSC_DRAW_BLACK);
308: PetscDrawLine(draw, binLeft, bins[0], binRight, bins[0], PETSC_DRAW_BLACK);
309: }
310: PetscDrawCollectiveEnd(draw);
311: } else {
312: numBins = hist->numBins;
313: numBinsOld = hist->numBins;
314: if (hist->integerBins && (((int)xmax - xmin) + 1.0e-05 > xmax - xmin)) {
315: initSize = (int)((int)xmax - xmin) / numBins;
316: while (initSize * numBins != (int)xmax - xmin) {
317: initSize = PetscMax(initSize - 1, 1);
318: numBins = (int)((int)xmax - xmin) / initSize;
319: PetscDrawHGSetNumberBins(hist, numBins);
320: }
321: }
322: binSize = (xmax - xmin) / numBins;
323: bins = hist->bins;
325: PetscArrayzero(bins, numBins);
327: maxHeight = 0.0;
328: for (i = 0; i < numBins; i++) {
329: binLeft = xmin + binSize * i;
330: binRight = xmin + binSize * (i + 1);
331: for (p = 0; p < numValues; p++) {
332: if ((values[p] >= binLeft) && (values[p] < binRight)) bins[i]++;
333: /* Handle last bin separately */
334: if ((i == numBins - 1) && (values[p] == binRight)) bins[i]++;
335: if (!i) {
336: mean += values[p];
337: var += values[p] * values[p];
338: }
339: }
340: maxHeight = PetscMax(maxHeight, bins[i]);
341: }
342: if (maxHeight > ymax) ymax = hist->ymax = maxHeight;
344: PetscDrawAxisSetLimits(hist->axis, xmin, xmax, ymin, ymax);
345: if (hist->calcStats) {
346: mean /= numValues;
347: if (numValues > 1) var = (var - numValues * mean * mean) / (numValues - 1);
348: else var = 0.0;
349: PetscSNPrintf(title, 256, "Mean: %g Var: %g", (double)mean, (double)var);
350: PetscSNPrintf(xlabel, 256, "Total: %" PetscInt_FMT, numValues);
351: PetscDrawAxisSetLabels(hist->axis, title, xlabel, NULL);
352: }
353: PetscDrawAxisDraw(hist->axis);
354: PetscDrawCollectiveBegin(draw);
355: if (rank == 0) { /* Draw bins */
356: for (i = 0; i < numBins; i++) {
357: binLeft = xmin + binSize * i;
358: binRight = xmin + binSize * (i + 1);
359: PetscDrawRectangle(draw, binLeft, ymin, binRight, bins[i], bcolor, bcolor, bcolor, bcolor);
360: PetscDrawLine(draw, binLeft, ymin, binLeft, bins[i], PETSC_DRAW_BLACK);
361: PetscDrawLine(draw, binRight, ymin, binRight, bins[i], PETSC_DRAW_BLACK);
362: PetscDrawLine(draw, binLeft, bins[i], binRight, bins[i], PETSC_DRAW_BLACK);
363: if (color == PETSC_DRAW_ROTATE && bins[i]) bcolor++;
364: if (bcolor > PETSC_DRAW_BASIC_COLORS - 1) bcolor = PETSC_DRAW_BLACK + 1;
365: }
366: }
367: PetscDrawCollectiveEnd(draw);
368: PetscDrawHGSetNumberBins(hist, numBinsOld);
369: }
371: PetscDrawFlush(draw);
372: PetscDrawPause(draw);
373: return 0;
374: }
376: /*@
377: PetscDrawHGSave - Saves a drawn image
379: Collective
381: Input Parameter:
382: . hist - The histogram context
384: Level: intermediate
386: .seealso: `PetscDrawSave()`, `PetscDrawHGCreate()`, `PetscDrawHGGetDraw()`, `PetscDrawSetSave()`, `PetscDrawSave()`, `PetscDrawHGDraw()`
387: @*/
388: PetscErrorCode PetscDrawHGSave(PetscDrawHG hg)
389: {
391: PetscDrawSave(hg->win);
392: return 0;
393: }
395: /*@
396: PetscDrawHGView - Prints the histogram information to a viewer
398: Not collective
400: Input Parameter:
401: . hist - The histogram context
403: Level: beginner
405: .seealso: `PetscDrawHG`, `PetscViewer`, `PetscDrawHGCreate()`, `PetscDrawHGGetDraw()`, `PetscDrawSetSave()`, `PetscDrawSave()`, `PetscDrawHGDraw()`
406: @*/
407: PetscErrorCode PetscDrawHGView(PetscDrawHG hist, PetscViewer viewer)
408: {
409: PetscReal xmax, xmin, *bins, *values, binSize, binLeft, binRight, mean, var;
410: PetscInt numBins, numBinsOld, numValues, initSize, i, p;
414: if ((hist->xmin > hist->xmax) || (hist->ymin >= hist->ymax)) return 0;
415: if (hist->numValues < 1) return 0;
417: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)hist), &viewer);
418: PetscObjectPrintClassNamePrefixType((PetscObject)hist, viewer);
419: xmax = hist->xmax;
420: xmin = hist->xmin;
421: numValues = hist->numValues;
422: values = hist->values;
423: mean = 0.0;
424: var = 0.0;
425: if (xmax == xmin) {
426: /* Calculate number of points in the bin */
427: bins = hist->bins;
428: bins[0] = 0.;
429: for (p = 0; p < numValues; p++) {
430: if (values[p] == xmin) bins[0]++;
431: mean += values[p];
432: var += values[p] * values[p];
433: }
434: /* Draw bins */
435: PetscViewerASCIIPrintf(viewer, "Bin %2d (%6.2g - %6.2g): %.0g\n", 0, (double)xmin, (double)xmax, (double)bins[0]);
436: } else {
437: numBins = hist->numBins;
438: numBinsOld = hist->numBins;
439: if (hist->integerBins && (((int)xmax - xmin) + 1.0e-05 > xmax - xmin)) {
440: initSize = (int)((int)xmax - xmin) / numBins;
441: while (initSize * numBins != (int)xmax - xmin) {
442: initSize = PetscMax(initSize - 1, 1);
443: numBins = (int)((int)xmax - xmin) / initSize;
444: PetscDrawHGSetNumberBins(hist, numBins);
445: }
446: }
447: binSize = (xmax - xmin) / numBins;
448: bins = hist->bins;
450: /* Calculate number of points in each bin */
451: PetscArrayzero(bins, numBins);
452: for (i = 0; i < numBins; i++) {
453: binLeft = xmin + binSize * i;
454: binRight = xmin + binSize * (i + 1);
455: for (p = 0; p < numValues; p++) {
456: if ((values[p] >= binLeft) && (values[p] < binRight)) bins[i]++;
457: /* Handle last bin separately */
458: if ((i == numBins - 1) && (values[p] == binRight)) bins[i]++;
459: if (!i) {
460: mean += values[p];
461: var += values[p] * values[p];
462: }
463: }
464: }
465: /* Draw bins */
466: for (i = 0; i < numBins; i++) {
467: binLeft = xmin + binSize * i;
468: binRight = xmin + binSize * (i + 1);
469: PetscViewerASCIIPrintf(viewer, "Bin %2d (%6.2g - %6.2g): %.0g\n", (int)i, (double)binLeft, (double)binRight, (double)bins[i]);
470: }
471: PetscDrawHGSetNumberBins(hist, numBinsOld);
472: }
474: if (hist->calcStats) {
475: mean /= numValues;
476: if (numValues > 1) var = (var - numValues * mean * mean) / (numValues - 1);
477: else var = 0.0;
478: PetscViewerASCIIPrintf(viewer, "Mean: %g Var: %g\n", (double)mean, (double)var);
479: PetscViewerASCIIPrintf(viewer, "Total: %" PetscInt_FMT "\n", numValues);
480: }
481: return 0;
482: }
484: /*@
485: PetscDrawHGSetColor - Sets the color the bars will be drawn with.
487: Logically Collective
489: Input Parameters:
490: + hist - The histogram context
491: - color - one of the colors defined in petscdraw.h or `PETSC_DRAW_ROTATE` to make each bar a
492: different color
494: Level: intermediate
496: .seealso: `PetscDrawHG`, `PetscDrawHGCreate()`, `PetscDrawHGGetDraw()`, `PetscDrawSetSave()`, `PetscDrawSave()`, `PetscDrawHGDraw()`, `PetscDrawHGGetAxis()`
497: @*/
498: PetscErrorCode PetscDrawHGSetColor(PetscDrawHG hist, int color)
499: {
502: hist->color = color;
503: return 0;
504: }
506: /*@
507: PetscDrawHGSetLimits - Sets the axis limits for a histogram. If more
508: points are added after this call, the limits will be adjusted to
509: include those additional points.
511: Logically Collective
513: Input Parameters:
514: + hist - The histogram context
515: - x_min,x_max,y_min,y_max - The limits
517: Level: intermediate
519: .seealso: `PetscDrawHG`, `PetscDrawHGCreate()`, `PetscDrawHGGetDraw()`, `PetscDrawSetSave()`, `PetscDrawSave()`, `PetscDrawHGDraw()`, `PetscDrawHGGetAxis()`
520: @*/
521: PetscErrorCode PetscDrawHGSetLimits(PetscDrawHG hist, PetscReal x_min, PetscReal x_max, int y_min, int y_max)
522: {
525: hist->xmin = x_min;
526: hist->xmax = x_max;
527: hist->ymin = y_min;
528: hist->ymax = y_max;
529: return 0;
530: }
532: /*@
533: PetscDrawHGCalcStats - Turns on calculation of descriptive statistics associated with the histogram
535: Not collective
537: Input Parameters:
538: + hist - The histogram context
539: - calc - Flag for calculation
541: Level: intermediate
543: .seealso: `PetscDrawHG`, `PetscDrawHGCreate()`, `PetscDrawHGAddValue()`, `PetscDrawHGView()`, `PetscDrawHGDraw()`
544: @*/
545: PetscErrorCode PetscDrawHGCalcStats(PetscDrawHG hist, PetscBool calc)
546: {
549: hist->calcStats = calc;
550: return 0;
551: }
553: /*@
554: PetscDrawHGIntegerBins - Turns on integer width bins
556: Not collective
558: Input Parameters:
559: + hist - The histogram context
560: - ints - Flag for integer width bins
562: Level: intermediate
564: .seealso: `PetscDrawHG`, `PetscDrawHGCreate()`, `PetscDrawHGAddValue()`, `PetscDrawHGView()`, `PetscDrawHGDraw()`, `PetscDrawHGSetColor()`
565: @*/
566: PetscErrorCode PetscDrawHGIntegerBins(PetscDrawHG hist, PetscBool ints)
567: {
570: hist->integerBins = ints;
571: return 0;
572: }
574: /*@C
575: PetscDrawHGGetAxis - Gets the axis context associated with a histogram.
576: This is useful if one wants to change some axis property, such as
577: labels, color, etc. The axis context should not be destroyed by the
578: application code.
580: Not Collective, axis is parallel if hist is parallel
582: Input Parameter:
583: . hist - The histogram context
585: Output Parameter:
586: . axis - The axis context
588: Level: intermediate
590: .seealso: `PetscDrawHG`, `PetscDrawAxis`, `PetscDrawHGCreate()`, `PetscDrawHGAddValue()`, `PetscDrawHGView()`, `PetscDrawHGDraw()`, `PetscDrawHGSetColor()`, `PetscDrawAxis`, `PetscDrawHGSetLimits()`
591: @*/
592: PetscErrorCode PetscDrawHGGetAxis(PetscDrawHG hist, PetscDrawAxis *axis)
593: {
596: *axis = hist->axis;
597: return 0;
598: }
600: /*@C
601: PetscDrawHGGetDraw - Gets the draw context associated with a histogram.
603: Not Collective, draw is parallel if hist is parallel
605: Input Parameter:
606: . hist - The histogram context
608: Output Parameter:
609: . draw - The draw context
611: Level: intermediate
613: .seealso: `PetscDraw`, `PetscDrawHG`, `PetscDrawHGCreate()`, `PetscDrawHGAddValue()`, `PetscDrawHGView()`, `PetscDrawHGDraw()`, `PetscDrawHGSetColor()`, `PetscDrawAxis`, `PetscDrawHGSetLimits()`
614: @*/
615: PetscErrorCode PetscDrawHGGetDraw(PetscDrawHG hist, PetscDraw *draw)
616: {
619: *draw = hist->win;
620: return 0;
621: }