Actual source code: hists.c
petsc-3.14.6 2021-03-30
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 on PetscDraw
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 processes 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()
58: @*/
59: PetscErrorCode PetscDrawHGCreate(PetscDraw draw,int bins,PetscDrawHG *hist)
60: {
61: PetscDrawHG h;
69: PetscHeaderCreate(h,PETSC_DRAWHG_CLASSID,"DrawHG","Histogram","Draw",PetscObjectComm((PetscObject)draw),PetscDrawHGDestroy,NULL);
70: PetscLogObjectParent((PetscObject)draw,(PetscObject)h);
72: PetscObjectReference((PetscObject)draw);
73: h->win = draw;
75: h->view = NULL;
76: h->destroy = NULL;
77: h->color = PETSC_DRAW_GREEN;
78: h->xmin = PETSC_MAX_REAL;
79: h->xmax = PETSC_MIN_REAL;
80: h->ymin = 0.;
81: h->ymax = 1.;
82: h->numBins = bins;
83: h->maxBins = bins;
85: PetscMalloc1(h->maxBins,&h->bins);
87: h->numValues = 0;
88: h->maxValues = CHUNKSIZE;
89: h->calcStats = PETSC_FALSE;
90: h->integerBins = PETSC_FALSE;
92: PetscMalloc1(h->maxValues,&h->values);
93: PetscLogObjectMemory((PetscObject)h,(h->maxBins + h->maxValues)*sizeof(PetscReal));
95: PetscDrawAxisCreate(draw,&h->axis);
96: PetscLogObjectParent((PetscObject)h,(PetscObject)h->axis);
98: *hist = h;
99: return(0);
100: }
102: /*@
103: PetscDrawHGSetNumberBins - Change the number of bins that are to be drawn.
105: Logically Collective on PetscDrawHG
107: Input Parameter:
108: + hist - The histogram context.
109: - bins - The number of bins.
111: Level: intermediate
113: .seealso: PetscDrawHGCreate(), PetscDrawHG, PetscDrawHGDraw(), PetscDrawHGIntegerBins()
115: @*/
116: PetscErrorCode PetscDrawHGSetNumberBins(PetscDrawHG hist, int bins)
117: {
124: if (hist->maxBins < bins) {
125: PetscFree(hist->bins);
126: PetscMalloc1(bins, &hist->bins);
127: PetscLogObjectMemory((PetscObject)hist, (bins - hist->maxBins) * sizeof(PetscReal));
128: hist->maxBins = bins;
129: }
130: hist->numBins = bins;
131: return(0);
132: }
134: /*@
135: PetscDrawHGReset - Clears histogram to allow for reuse with new data.
137: Logically Collective on PetscDrawHG
139: Input Parameter:
140: . hist - The histogram context.
142: Level: intermediate
144: .seealso: PetscDrawHGCreate(), PetscDrawHG, PetscDrawHGDraw(), PetscDrawHGAddValue()
146: @*/
147: PetscErrorCode PetscDrawHGReset(PetscDrawHG hist)
148: {
152: hist->xmin = PETSC_MAX_REAL;
153: hist->xmax = PETSC_MIN_REAL;
154: hist->ymin = 0.0;
155: hist->ymax = 0.0;
156: hist->numValues = 0;
157: return(0);
158: }
160: /*@C
161: PetscDrawHGDestroy - Frees all space taken up by histogram data structure.
163: Collective on PetscDrawHG
165: Input Parameter:
166: . hist - The histogram context
168: Level: intermediate
170: .seealso: PetscDrawHGCreate(), PetscDrawHG
171: @*/
172: PetscErrorCode PetscDrawHGDestroy(PetscDrawHG *hist)
173: {
177: if (!*hist) return(0);
179: if (--((PetscObject)(*hist))->refct > 0) {*hist = NULL; return(0);}
181: PetscFree((*hist)->bins);
182: PetscFree((*hist)->values);
183: PetscDrawAxisDestroy(&(*hist)->axis);
184: PetscDrawDestroy(&(*hist)->win);
185: PetscHeaderDestroy(hist);
186: return(0);
187: }
189: /*@
190: PetscDrawHGAddValue - Adds another value to the histogram.
192: Logically Collective on PetscDrawHG
194: Input Parameters:
195: + hist - The histogram
196: - value - The value
198: Level: intermediate
200: .seealso: PetscDrawHGCreate(), PetscDrawHG, PetscDrawHGDraw(), PetscDrawHGAddValue(), PetscDrawHGReset()
201: @*/
202: PetscErrorCode PetscDrawHGAddValue(PetscDrawHG hist, PetscReal value)
203: {
207: /* Allocate more memory if necessary */
208: if (hist->numValues >= hist->maxValues) {
209: PetscReal *tmp;
212: PetscMalloc1(hist->maxValues+CHUNKSIZE, &tmp);
213: PetscLogObjectMemory((PetscObject)hist, CHUNKSIZE * sizeof(PetscReal));
214: PetscArraycpy(tmp, hist->values, hist->maxValues);
215: PetscFree(hist->values);
217: hist->values = tmp;
218: hist->maxValues += CHUNKSIZE;
219: }
220: /* I disagree with the original Petsc implementation here. There should be no overshoot, but rather the
221: stated convention of using half-open intervals (always the way to go) */
222: if (!hist->numValues) {
223: hist->xmin = value;
224: hist->xmax = value;
225: #if 1
226: } else {
227: /* Update limits */
228: if (value > hist->xmax) hist->xmax = value;
229: if (value < hist->xmin) hist->xmin = value;
230: #else
231: } else if (hist->numValues == 1) {
232: /* Update limits -- We need to overshoot the largest value somewhat */
233: if (value > hist->xmax) hist->xmax = value + 0.001*(value - hist->xmin)/hist->numBins;
234: if (value < hist->xmin) {
235: hist->xmin = value;
236: hist->xmax = hist->xmax + 0.001*(hist->xmax - hist->xmin)/hist->numBins;
237: }
238: } else {
239: /* Update limits -- We need to overshoot the largest value somewhat */
240: if (value > hist->xmax) hist->xmax = value + 0.001*(hist->xmax - hist->xmin)/hist->numBins;
241: if (value < hist->xmin) hist->xmin = value;
242: #endif
243: }
245: hist->values[hist->numValues++] = value;
246: return(0);
247: }
249: /*@
250: PetscDrawHGDraw - Redraws a histogram.
252: Collective on PetscDrawHG
254: Input Parameter:
255: . hist - The histogram context
257: Level: intermediate
259: .seealso: PetscDrawHGCreate(), PetscDrawHG, PetscDrawHGDraw(), PetscDrawHGAddValue(), PetscDrawHGReset()
261: @*/
262: PetscErrorCode PetscDrawHGDraw(PetscDrawHG hist)
263: {
264: PetscDraw draw;
265: PetscBool isnull;
266: PetscReal xmin,xmax,ymin,ymax,*bins,*values,binSize,binLeft,binRight,maxHeight,mean,var;
267: char title[256];
268: char xlabel[256];
269: PetscInt numBins,numBinsOld,numValues,initSize,i,p,bcolor,color;
270: PetscMPIInt rank;
275: PetscDrawIsNull(hist->win,&isnull);
276: if (isnull) return(0);
277: MPI_Comm_rank(PetscObjectComm((PetscObject)hist),&rank);
279: if ((hist->xmin >= hist->xmax) || (hist->ymin >= hist->ymax)) return(0);
280: if (hist->numValues < 1) return(0);
282: color = hist->color;
283: if (color == PETSC_DRAW_ROTATE) bcolor = PETSC_DRAW_BLACK+1;
284: else bcolor = color;
286: xmin = hist->xmin;
287: xmax = hist->xmax;
288: ymin = hist->ymin;
289: ymax = hist->ymax;
290: numValues = hist->numValues;
291: values = hist->values;
292: mean = 0.0;
293: var = 0.0;
295: draw = hist->win;
296: PetscDrawCheckResizedWindow(draw);
297: PetscDrawClear(draw);
299: if (xmin == xmax) {
300: /* Calculate number of points in each bin */
301: bins = hist->bins;
302: bins[0] = 0.;
303: for (p = 0; p < numValues; p++) {
304: if (values[p] == xmin) bins[0]++;
305: mean += values[p];
306: var += values[p]*values[p];
307: }
308: maxHeight = bins[0];
309: if (maxHeight > ymax) ymax = hist->ymax = maxHeight;
310: xmax = xmin + 1;
311: PetscDrawAxisSetLimits(hist->axis, xmin, xmax, ymin, ymax);
312: if (hist->calcStats) {
313: mean /= numValues;
314: if (numValues > 1) var = (var - numValues*mean*mean) / (numValues-1);
315: else var = 0.0;
316: PetscSNPrintf(title, 256, "Mean: %g Var: %g", (double)mean, (double)var);
317: PetscSNPrintf(xlabel,256, "Total: %D", numValues);
318: PetscDrawAxisSetLabels(hist->axis, title, xlabel, NULL);
319: }
320: PetscDrawAxisDraw(hist->axis);
321: PetscDrawCollectiveBegin(draw);
322: if (!rank) { /* Draw bins */
323: binLeft = xmin;
324: binRight = xmax;
325: PetscDrawRectangle(draw,binLeft,ymin,binRight,bins[0],bcolor,bcolor,bcolor,bcolor);
326: PetscDrawLine(draw,binLeft,ymin,binLeft,bins[0],PETSC_DRAW_BLACK);
327: PetscDrawLine(draw,binRight,ymin,binRight,bins[0],PETSC_DRAW_BLACK);
328: PetscDrawLine(draw,binLeft,bins[0],binRight,bins[0],PETSC_DRAW_BLACK);
329: }
330: PetscDrawCollectiveEnd(draw);
331: } else {
332: numBins = hist->numBins;
333: numBinsOld = hist->numBins;
334: if (hist->integerBins && (((int) xmax - xmin) + 1.0e-05 > xmax - xmin)) {
335: initSize = (int) ((int) xmax - xmin)/numBins;
336: while (initSize*numBins != (int) xmax - xmin) {
337: initSize = PetscMax(initSize - 1, 1);
338: numBins = (int) ((int) xmax - xmin)/initSize;
339: PetscDrawHGSetNumberBins(hist, numBins);
340: }
341: }
342: binSize = (xmax - xmin)/numBins;
343: bins = hist->bins;
345: PetscArrayzero(bins, numBins);
347: maxHeight = 0.0;
348: for (i = 0; i < numBins; i++) {
349: binLeft = xmin + binSize*i;
350: binRight = xmin + binSize*(i+1);
351: for (p = 0; p < numValues; p++) {
352: if ((values[p] >= binLeft) && (values[p] < binRight)) bins[i]++;
353: /* Handle last bin separately */
354: if ((i == numBins-1) && (values[p] == binRight)) bins[i]++;
355: if (!i) {
356: mean += values[p];
357: var += values[p]*values[p];
358: }
359: }
360: maxHeight = PetscMax(maxHeight, bins[i]);
361: }
362: if (maxHeight > ymax) ymax = hist->ymax = maxHeight;
364: PetscDrawAxisSetLimits(hist->axis, xmin, xmax, ymin, ymax);
365: if (hist->calcStats) {
366: mean /= numValues;
367: if (numValues > 1) var = (var - numValues*mean*mean) / (numValues-1);
368: else var = 0.0;
369: PetscSNPrintf(title, 256,"Mean: %g Var: %g", (double)mean, (double)var);
370: PetscSNPrintf(xlabel,256, "Total: %D", numValues);
371: PetscDrawAxisSetLabels(hist->axis, title, xlabel, NULL);
372: }
373: PetscDrawAxisDraw(hist->axis);
374: PetscDrawCollectiveBegin(draw);
375: if (!rank) { /* Draw bins */
376: for (i = 0; i < numBins; i++) {
377: binLeft = xmin + binSize*i;
378: binRight = xmin + binSize*(i+1);
379: PetscDrawRectangle(draw,binLeft,ymin,binRight,bins[i],bcolor,bcolor,bcolor,bcolor);
380: PetscDrawLine(draw,binLeft,ymin,binLeft,bins[i],PETSC_DRAW_BLACK);
381: PetscDrawLine(draw,binRight,ymin,binRight,bins[i],PETSC_DRAW_BLACK);
382: PetscDrawLine(draw,binLeft,bins[i],binRight,bins[i],PETSC_DRAW_BLACK);
383: if (color == PETSC_DRAW_ROTATE && bins[i]) bcolor++;
384: if (bcolor > PETSC_DRAW_BASIC_COLORS-1) bcolor = PETSC_DRAW_BLACK+1;
385: }
386: }
387: PetscDrawCollectiveEnd(draw);
388: PetscDrawHGSetNumberBins(hist,numBinsOld);
389: }
391: PetscDrawFlush(draw);
392: PetscDrawPause(draw);
393: return(0);
394: }
396: /*@
397: PetscDrawHGSave - Saves a drawn image
399: Collective on PetscDrawHG
401: Input Parameter:
402: . hist - The histogram context
404: Level: intermediate
406: .seealso: PetscDrawHGCreate(), PetscDrawHGGetDraw(), PetscDrawSetSave(), PetscDrawSave(), PetscDrawHGDraw()
407: @*/
408: PetscErrorCode PetscDrawHGSave(PetscDrawHG hg)
409: {
414: PetscDrawSave(hg->win);
415: return(0);
416: }
418: /*@
419: PetscDrawHGView - Prints the histogram information.
421: Not collective
423: Input Parameter:
424: . hist - The histogram context
426: Level: beginner
428: .seealso: PetscDrawHGCreate(), PetscDrawHGGetDraw(), PetscDrawSetSave(), PetscDrawSave(), PetscDrawHGDraw()
430: @*/
431: PetscErrorCode PetscDrawHGView(PetscDrawHG hist,PetscViewer viewer)
432: {
433: PetscReal xmax,xmin,*bins,*values,binSize,binLeft,binRight,mean,var;
435: PetscInt numBins,numBinsOld,numValues,initSize,i,p;
440: if ((hist->xmin > hist->xmax) || (hist->ymin >= hist->ymax)) return(0);
441: if (hist->numValues < 1) return(0);
443: if (!viewer){
444: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)hist),&viewer);
445: }
446: PetscObjectPrintClassNamePrefixType((PetscObject)hist,viewer);
447: xmax = hist->xmax;
448: xmin = hist->xmin;
449: numValues = hist->numValues;
450: values = hist->values;
451: mean = 0.0;
452: var = 0.0;
453: if (xmax == xmin) {
454: /* Calculate number of points in the bin */
455: bins = hist->bins;
456: bins[0] = 0.;
457: for (p = 0; p < numValues; p++) {
458: if (values[p] == xmin) bins[0]++;
459: mean += values[p];
460: var += values[p]*values[p];
461: }
462: /* Draw bins */
463: PetscViewerASCIIPrintf(viewer, "Bin %2d (%6.2g - %6.2g): %.0g\n", 0, (double)xmin, (double)xmax, (double)bins[0]);
464: } else {
465: numBins = hist->numBins;
466: numBinsOld = hist->numBins;
467: if (hist->integerBins && (((int) xmax - xmin) + 1.0e-05 > xmax - xmin)) {
468: initSize = (int) ((int) xmax - xmin)/numBins;
469: while (initSize*numBins != (int) xmax - xmin) {
470: initSize = PetscMax(initSize - 1, 1);
471: numBins = (int) ((int) xmax - xmin)/initSize;
472: PetscDrawHGSetNumberBins(hist, numBins);
473: }
474: }
475: binSize = (xmax - xmin)/numBins;
476: bins = hist->bins;
478: /* Calculate number of points in each bin */
479: PetscArrayzero(bins, numBins);
480: for (i = 0; i < numBins; i++) {
481: binLeft = xmin + binSize*i;
482: binRight = xmin + binSize*(i+1);
483: for (p = 0; p < numValues; p++) {
484: if ((values[p] >= binLeft) && (values[p] < binRight)) bins[i]++;
485: /* Handle last bin separately */
486: if ((i == numBins-1) && (values[p] == binRight)) bins[i]++;
487: if (!i) {
488: mean += values[p];
489: var += values[p]*values[p];
490: }
491: }
492: }
493: /* Draw bins */
494: for (i = 0; i < numBins; i++) {
495: binLeft = xmin + binSize*i;
496: binRight = xmin + binSize*(i+1);
497: PetscViewerASCIIPrintf(viewer, "Bin %2d (%6.2g - %6.2g): %.0g\n", (int)i, (double)binLeft, (double)binRight, (double)bins[i]);
498: }
499: PetscDrawHGSetNumberBins(hist, numBinsOld);
500: }
502: if (hist->calcStats) {
503: mean /= numValues;
504: if (numValues > 1) var = (var - numValues*mean*mean) / (numValues-1);
505: else var = 0.0;
506: PetscViewerASCIIPrintf(viewer, "Mean: %g Var: %g\n", (double)mean, (double)var);
507: PetscViewerASCIIPrintf(viewer, "Total: %D\n", numValues);
508: }
509: return(0);
510: }
512: /*@
513: PetscDrawHGSetColor - Sets the color the bars will be drawn with.
515: Logically Collective on PetscDrawHG
517: Input Parameters:
518: + hist - The histogram context
519: - color - one of the colors defined in petscdraw.h or PETSC_DRAW_ROTATE to make each bar a
520: different color
522: Level: intermediate
524: .seealso: PetscDrawHGCreate(), PetscDrawHGGetDraw(), PetscDrawSetSave(), PetscDrawSave(), PetscDrawHGDraw(), PetscDrawHGGetAxis()
526: @*/
527: PetscErrorCode PetscDrawHGSetColor(PetscDrawHG hist,int color)
528: {
532: hist->color = color;
533: return(0);
534: }
536: /*@
537: PetscDrawHGSetLimits - Sets the axis limits for a histogram. If more
538: points are added after this call, the limits will be adjusted to
539: include those additional points.
541: Logically Collective on PetscDrawHG
543: Input Parameters:
544: + hist - The histogram context
545: - x_min,x_max,y_min,y_max - The limits
547: Level: intermediate
549: .seealso: PetscDrawHGCreate(), PetscDrawHGGetDraw(), PetscDrawSetSave(), PetscDrawSave(), PetscDrawHGDraw(), PetscDrawHGGetAxis()
551: @*/
552: PetscErrorCode PetscDrawHGSetLimits(PetscDrawHG hist, PetscReal x_min, PetscReal x_max, int y_min, int y_max)
553: {
557: hist->xmin = x_min;
558: hist->xmax = x_max;
559: hist->ymin = y_min;
560: hist->ymax = y_max;
561: return(0);
562: }
564: /*@
565: PetscDrawHGCalcStats - Turns on calculation of descriptive statistics
567: Not collective
569: Input Parameters:
570: + hist - The histogram context
571: - calc - Flag for calculation
573: Level: intermediate
575: .seealso: PetscDrawHGCreate(), PetscDrawHGAddValue(), PetscDrawHGView(), PetscDrawHGDraw()
577: @*/
578: PetscErrorCode PetscDrawHGCalcStats(PetscDrawHG hist, PetscBool calc)
579: {
583: hist->calcStats = calc;
584: return(0);
585: }
587: /*@
588: PetscDrawHGIntegerBins - Turns on integer width bins
590: Not collective
592: Input Parameters:
593: + hist - The histogram context
594: - ints - Flag for integer width bins
596: Level: intermediate
598: .seealso: PetscDrawHGCreate(), PetscDrawHGAddValue(), PetscDrawHGView(), PetscDrawHGDraw(), PetscDrawHGSetColor()
600: @*/
601: PetscErrorCode PetscDrawHGIntegerBins(PetscDrawHG hist, PetscBool ints)
602: {
606: hist->integerBins = ints;
607: return(0);
608: }
610: /*@C
611: PetscDrawHGGetAxis - Gets the axis context associated with a histogram.
612: This is useful if one wants to change some axis property, such as
613: labels, color, etc. The axis context should not be destroyed by the
614: application code.
616: Not Collective, PetscDrawAxis is parallel if PetscDrawHG is parallel
618: Input Parameter:
619: . hist - The histogram context
621: Output Parameter:
622: . axis - The axis context
624: Level: intermediate
626: .seealso: PetscDrawHGCreate(), PetscDrawHGAddValue(), PetscDrawHGView(), PetscDrawHGDraw(), PetscDrawHGSetColor(), PetscDrawAxis, PetscDrawHGSetLimits()
628: @*/
629: PetscErrorCode PetscDrawHGGetAxis(PetscDrawHG hist,PetscDrawAxis *axis)
630: {
634: *axis = hist->axis;
635: return(0);
636: }
638: /*@C
639: PetscDrawHGGetDraw - Gets the draw context associated with a histogram.
641: Not Collective, PetscDraw is parallel if PetscDrawHG is parallel
643: Input Parameter:
644: . hist - The histogram context
646: Output Parameter:
647: . draw - The draw context
649: Level: intermediate
651: .seealso: PetscDrawHGCreate(), PetscDrawHGAddValue(), PetscDrawHGView(), PetscDrawHGDraw(), PetscDrawHGSetColor(), PetscDrawAxis, PetscDrawHGSetLimits()
653: @*/
654: PetscErrorCode PetscDrawHGGetDraw(PetscDrawHG hist,PetscDraw *draw)
655: {
659: *draw = hist->win;
660: return(0);
661: }