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 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;
67: PetscHeaderCreate(h,PETSC_DRAWHG_CLASSID,"DrawHG","Histogram","Draw",PetscObjectComm((PetscObject)draw),PetscDrawHGDestroy,NULL);
68: PetscLogObjectParent((PetscObject)draw,(PetscObject)h);
70: PetscObjectReference((PetscObject)draw);
71: h->win = draw;
73: h->view = NULL;
74: h->destroy = NULL;
75: h->color = PETSC_DRAW_GREEN;
76: h->xmin = PETSC_MAX_REAL;
77: h->xmax = PETSC_MIN_REAL;
78: h->ymin = 0.;
79: h->ymax = 1.;
80: h->numBins = bins;
81: h->maxBins = bins;
83: PetscMalloc1(h->maxBins,&h->bins);
85: h->numValues = 0;
86: h->maxValues = CHUNKSIZE;
87: h->calcStats = PETSC_FALSE;
88: h->integerBins = PETSC_FALSE;
90: PetscMalloc1(h->maxValues,&h->values);
91: PetscLogObjectMemory((PetscObject)h,(h->maxBins + h->maxValues)*sizeof(PetscReal));
93: PetscDrawAxisCreate(draw,&h->axis);
94: PetscLogObjectParent((PetscObject)h,(PetscObject)h->axis);
96: *hist = h;
97: return 0;
98: }
100: /*@
101: PetscDrawHGSetNumberBins - Change the number of bins that are to be drawn.
103: Logically Collective on PetscDrawHG
105: Input Parameters:
106: + hist - The histogram context.
107: - bins - The number of bins.
109: Level: intermediate
111: .seealso: PetscDrawHGCreate(), PetscDrawHG, PetscDrawHGDraw(), PetscDrawHGIntegerBins()
113: @*/
114: PetscErrorCode PetscDrawHGSetNumberBins(PetscDrawHG hist, int bins)
115: {
119: if (hist->maxBins < bins) {
120: PetscFree(hist->bins);
121: PetscMalloc1(bins, &hist->bins);
122: PetscLogObjectMemory((PetscObject)hist, (bins - hist->maxBins) * sizeof(PetscReal));
123: hist->maxBins = bins;
124: }
125: hist->numBins = bins;
126: return 0;
127: }
129: /*@
130: PetscDrawHGReset - Clears histogram to allow for reuse with new data.
132: Logically Collective on PetscDrawHG
134: Input Parameter:
135: . hist - The histogram context.
137: Level: intermediate
139: .seealso: PetscDrawHGCreate(), PetscDrawHG, PetscDrawHGDraw(), PetscDrawHGAddValue()
141: @*/
142: PetscErrorCode PetscDrawHGReset(PetscDrawHG hist)
143: {
146: hist->xmin = PETSC_MAX_REAL;
147: hist->xmax = PETSC_MIN_REAL;
148: hist->ymin = 0.0;
149: hist->ymax = 0.0;
150: hist->numValues = 0;
151: return 0;
152: }
154: /*@C
155: PetscDrawHGDestroy - Frees all space taken up by histogram data structure.
157: Collective on PetscDrawHG
159: Input Parameter:
160: . hist - The histogram context
162: Level: intermediate
164: .seealso: PetscDrawHGCreate(), PetscDrawHG
165: @*/
166: PetscErrorCode PetscDrawHGDestroy(PetscDrawHG *hist)
167: {
168: if (!*hist) return 0;
170: if (--((PetscObject)(*hist))->refct > 0) {*hist = NULL; return 0;}
172: PetscFree((*hist)->bins);
173: PetscFree((*hist)->values);
174: PetscDrawAxisDestroy(&(*hist)->axis);
175: PetscDrawDestroy(&(*hist)->win);
176: PetscHeaderDestroy(hist);
177: return 0;
178: }
180: /*@
181: PetscDrawHGAddValue - Adds another value to the histogram.
183: Logically Collective on PetscDrawHG
185: Input Parameters:
186: + hist - The histogram
187: - value - The value
189: Level: intermediate
191: .seealso: PetscDrawHGCreate(), PetscDrawHG, PetscDrawHGDraw(), PetscDrawHGAddValue(), PetscDrawHGReset()
192: @*/
193: PetscErrorCode PetscDrawHGAddValue(PetscDrawHG hist, PetscReal value)
194: {
197: /* Allocate more memory if necessary */
198: if (hist->numValues >= hist->maxValues) {
199: PetscReal *tmp;
201: PetscMalloc1(hist->maxValues+CHUNKSIZE, &tmp);
202: PetscLogObjectMemory((PetscObject)hist, CHUNKSIZE * sizeof(PetscReal));
203: PetscArraycpy(tmp, hist->values, hist->maxValues);
204: PetscFree(hist->values);
206: hist->values = tmp;
207: hist->maxValues += CHUNKSIZE;
208: }
209: /* I disagree with the original Petsc implementation here. There should be no overshoot, but rather the
210: stated convention of using half-open intervals (always the way to go) */
211: if (!hist->numValues && (hist->xmin == PETSC_MAX_REAL) && (hist->xmax == PETSC_MIN_REAL)) {
212: hist->xmin = value;
213: hist->xmax = value;
214: #if 1
215: } else {
216: /* Update limits */
217: if (value > hist->xmax) hist->xmax = value;
218: if (value < hist->xmin) hist->xmin = value;
219: #else
220: } else if (hist->numValues == 1) {
221: /* Update limits -- We need to overshoot the largest value somewhat */
222: if (value > hist->xmax) hist->xmax = value + 0.001*(value - hist->xmin)/hist->numBins;
223: if (value < hist->xmin) {
224: hist->xmin = value;
225: hist->xmax = hist->xmax + 0.001*(hist->xmax - hist->xmin)/hist->numBins;
226: }
227: } else {
228: /* Update limits -- We need to overshoot the largest value somewhat */
229: if (value > hist->xmax) hist->xmax = value + 0.001*(hist->xmax - hist->xmin)/hist->numBins;
230: if (value < hist->xmin) hist->xmin = value;
231: #endif
232: }
234: hist->values[hist->numValues++] = value;
235: return 0;
236: }
238: /*@
239: PetscDrawHGDraw - Redraws a histogram.
241: Collective on PetscDrawHG
243: Input Parameter:
244: . hist - The histogram context
246: Level: intermediate
248: .seealso: PetscDrawHGCreate(), PetscDrawHG, PetscDrawHGDraw(), PetscDrawHGAddValue(), PetscDrawHGReset()
250: @*/
251: PetscErrorCode PetscDrawHGDraw(PetscDrawHG hist)
252: {
253: PetscDraw draw;
254: PetscBool isnull;
255: PetscReal xmin,xmax,ymin,ymax,*bins,*values,binSize,binLeft,binRight,maxHeight,mean,var;
256: char title[256];
257: char xlabel[256];
258: PetscInt numBins,numBinsOld,numValues,initSize,i,p,bcolor,color;
259: PetscMPIInt rank;
263: PetscDrawIsNull(hist->win,&isnull);
264: if (isnull) return 0;
265: MPI_Comm_rank(PetscObjectComm((PetscObject)hist),&rank);
267: if ((hist->xmin >= hist->xmax) || (hist->ymin >= hist->ymax)) return 0;
268: if (hist->numValues < 1) return 0;
270: color = hist->color;
271: if (color == PETSC_DRAW_ROTATE) bcolor = PETSC_DRAW_BLACK+1;
272: else bcolor = color;
274: xmin = hist->xmin;
275: xmax = hist->xmax;
276: ymin = hist->ymin;
277: ymax = hist->ymax;
278: numValues = hist->numValues;
279: values = hist->values;
280: mean = 0.0;
281: var = 0.0;
283: draw = hist->win;
284: PetscDrawCheckResizedWindow(draw);
285: PetscDrawClear(draw);
287: if (xmin == xmax) {
288: /* Calculate number of points in each bin */
289: bins = hist->bins;
290: bins[0] = 0.;
291: for (p = 0; p < numValues; p++) {
292: if (values[p] == xmin) bins[0]++;
293: mean += values[p];
294: var += values[p]*values[p];
295: }
296: maxHeight = bins[0];
297: if (maxHeight > ymax) ymax = hist->ymax = maxHeight;
298: xmax = xmin + 1;
299: PetscDrawAxisSetLimits(hist->axis, xmin, xmax, ymin, ymax);
300: if (hist->calcStats) {
301: mean /= numValues;
302: if (numValues > 1) var = (var - numValues*mean*mean) / (numValues-1);
303: else var = 0.0;
304: PetscSNPrintf(title, 256, "Mean: %g Var: %g", (double)mean, (double)var);
305: PetscSNPrintf(xlabel,256, "Total: %" PetscInt_FMT, numValues);
306: PetscDrawAxisSetLabels(hist->axis, title, xlabel, NULL);
307: }
308: PetscDrawAxisDraw(hist->axis);
309: PetscDrawCollectiveBegin(draw);
310: if (rank == 0) { /* Draw bins */
311: binLeft = xmin;
312: binRight = xmax;
313: PetscDrawRectangle(draw,binLeft,ymin,binRight,bins[0],bcolor,bcolor,bcolor,bcolor);
314: PetscDrawLine(draw,binLeft,ymin,binLeft,bins[0],PETSC_DRAW_BLACK);
315: PetscDrawLine(draw,binRight,ymin,binRight,bins[0],PETSC_DRAW_BLACK);
316: PetscDrawLine(draw,binLeft,bins[0],binRight,bins[0],PETSC_DRAW_BLACK);
317: }
318: PetscDrawCollectiveEnd(draw);
319: } else {
320: numBins = hist->numBins;
321: numBinsOld = hist->numBins;
322: if (hist->integerBins && (((int) xmax - xmin) + 1.0e-05 > xmax - xmin)) {
323: initSize = (int) ((int) xmax - xmin)/numBins;
324: while (initSize*numBins != (int) xmax - xmin) {
325: initSize = PetscMax(initSize - 1, 1);
326: numBins = (int) ((int) xmax - xmin)/initSize;
327: PetscDrawHGSetNumberBins(hist, numBins);
328: }
329: }
330: binSize = (xmax - xmin)/numBins;
331: bins = hist->bins;
333: PetscArrayzero(bins, numBins);
335: maxHeight = 0.0;
336: for (i = 0; i < numBins; i++) {
337: binLeft = xmin + binSize*i;
338: binRight = xmin + binSize*(i+1);
339: for (p = 0; p < numValues; p++) {
340: if ((values[p] >= binLeft) && (values[p] < binRight)) bins[i]++;
341: /* Handle last bin separately */
342: if ((i == numBins-1) && (values[p] == binRight)) bins[i]++;
343: if (!i) {
344: mean += values[p];
345: var += values[p]*values[p];
346: }
347: }
348: maxHeight = PetscMax(maxHeight, bins[i]);
349: }
350: if (maxHeight > ymax) ymax = hist->ymax = maxHeight;
352: PetscDrawAxisSetLimits(hist->axis, xmin, xmax, ymin, ymax);
353: if (hist->calcStats) {
354: mean /= numValues;
355: if (numValues > 1) var = (var - numValues*mean*mean) / (numValues-1);
356: else var = 0.0;
357: PetscSNPrintf(title, 256,"Mean: %g Var: %g", (double)mean, (double)var);
358: PetscSNPrintf(xlabel,256, "Total: %" PetscInt_FMT, numValues);
359: PetscDrawAxisSetLabels(hist->axis, title, xlabel, NULL);
360: }
361: PetscDrawAxisDraw(hist->axis);
362: PetscDrawCollectiveBegin(draw);
363: if (rank == 0) { /* Draw bins */
364: for (i = 0; i < numBins; i++) {
365: binLeft = xmin + binSize*i;
366: binRight = xmin + binSize*(i+1);
367: PetscDrawRectangle(draw,binLeft,ymin,binRight,bins[i],bcolor,bcolor,bcolor,bcolor);
368: PetscDrawLine(draw,binLeft,ymin,binLeft,bins[i],PETSC_DRAW_BLACK);
369: PetscDrawLine(draw,binRight,ymin,binRight,bins[i],PETSC_DRAW_BLACK);
370: PetscDrawLine(draw,binLeft,bins[i],binRight,bins[i],PETSC_DRAW_BLACK);
371: if (color == PETSC_DRAW_ROTATE && bins[i]) bcolor++;
372: if (bcolor > PETSC_DRAW_BASIC_COLORS-1) bcolor = PETSC_DRAW_BLACK+1;
373: }
374: }
375: PetscDrawCollectiveEnd(draw);
376: PetscDrawHGSetNumberBins(hist,numBinsOld);
377: }
379: PetscDrawFlush(draw);
380: PetscDrawPause(draw);
381: return 0;
382: }
384: /*@
385: PetscDrawHGSave - Saves a drawn image
387: Collective on PetscDrawHG
389: Input Parameter:
390: . hist - The histogram context
392: Level: intermediate
394: .seealso: PetscDrawHGCreate(), PetscDrawHGGetDraw(), PetscDrawSetSave(), PetscDrawSave(), PetscDrawHGDraw()
395: @*/
396: PetscErrorCode PetscDrawHGSave(PetscDrawHG hg)
397: {
399: PetscDrawSave(hg->win);
400: return 0;
401: }
403: /*@
404: PetscDrawHGView - Prints the histogram information.
406: Not collective
408: Input Parameter:
409: . hist - The histogram context
411: Level: beginner
413: .seealso: PetscDrawHGCreate(), PetscDrawHGGetDraw(), PetscDrawSetSave(), PetscDrawSave(), PetscDrawHGDraw()
415: @*/
416: PetscErrorCode PetscDrawHGView(PetscDrawHG hist,PetscViewer viewer)
417: {
418: PetscReal xmax,xmin,*bins,*values,binSize,binLeft,binRight,mean,var;
419: PetscInt numBins,numBinsOld,numValues,initSize,i,p;
423: if ((hist->xmin > hist->xmax) || (hist->ymin >= hist->ymax)) return 0;
424: if (hist->numValues < 1) return 0;
426: if (!viewer) {
427: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)hist),&viewer);
428: }
429: PetscObjectPrintClassNamePrefixType((PetscObject)hist,viewer);
430: xmax = hist->xmax;
431: xmin = hist->xmin;
432: numValues = hist->numValues;
433: values = hist->values;
434: mean = 0.0;
435: var = 0.0;
436: if (xmax == xmin) {
437: /* Calculate number of points in the bin */
438: bins = hist->bins;
439: bins[0] = 0.;
440: for (p = 0; p < numValues; p++) {
441: if (values[p] == xmin) bins[0]++;
442: mean += values[p];
443: var += values[p]*values[p];
444: }
445: /* Draw bins */
446: PetscViewerASCIIPrintf(viewer, "Bin %2d (%6.2g - %6.2g): %.0g\n", 0, (double)xmin, (double)xmax, (double)bins[0]);
447: } else {
448: numBins = hist->numBins;
449: numBinsOld = hist->numBins;
450: if (hist->integerBins && (((int) xmax - xmin) + 1.0e-05 > xmax - xmin)) {
451: initSize = (int) ((int) xmax - xmin)/numBins;
452: while (initSize*numBins != (int) xmax - xmin) {
453: initSize = PetscMax(initSize - 1, 1);
454: numBins = (int) ((int) xmax - xmin)/initSize;
455: PetscDrawHGSetNumberBins(hist, numBins);
456: }
457: }
458: binSize = (xmax - xmin)/numBins;
459: bins = hist->bins;
461: /* Calculate number of points in each bin */
462: PetscArrayzero(bins, numBins);
463: for (i = 0; i < numBins; i++) {
464: binLeft = xmin + binSize*i;
465: binRight = xmin + binSize*(i+1);
466: for (p = 0; p < numValues; p++) {
467: if ((values[p] >= binLeft) && (values[p] < binRight)) bins[i]++;
468: /* Handle last bin separately */
469: if ((i == numBins-1) && (values[p] == binRight)) bins[i]++;
470: if (!i) {
471: mean += values[p];
472: var += values[p]*values[p];
473: }
474: }
475: }
476: /* Draw bins */
477: for (i = 0; i < numBins; i++) {
478: binLeft = xmin + binSize*i;
479: binRight = xmin + binSize*(i+1);
480: PetscViewerASCIIPrintf(viewer, "Bin %2d (%6.2g - %6.2g): %.0g\n", (int)i, (double)binLeft, (double)binRight, (double)bins[i]);
481: }
482: PetscDrawHGSetNumberBins(hist, numBinsOld);
483: }
485: if (hist->calcStats) {
486: mean /= numValues;
487: if (numValues > 1) var = (var - numValues*mean*mean) / (numValues-1);
488: else var = 0.0;
489: PetscViewerASCIIPrintf(viewer, "Mean: %g Var: %g\n", (double)mean, (double)var);
490: PetscViewerASCIIPrintf(viewer, "Total: %" PetscInt_FMT "\n", numValues);
491: }
492: return 0;
493: }
495: /*@
496: PetscDrawHGSetColor - Sets the color the bars will be drawn with.
498: Logically Collective on PetscDrawHG
500: Input Parameters:
501: + hist - The histogram context
502: - color - one of the colors defined in petscdraw.h or PETSC_DRAW_ROTATE to make each bar a
503: different color
505: Level: intermediate
507: .seealso: PetscDrawHGCreate(), PetscDrawHGGetDraw(), PetscDrawSetSave(), PetscDrawSave(), PetscDrawHGDraw(), PetscDrawHGGetAxis()
509: @*/
510: PetscErrorCode PetscDrawHGSetColor(PetscDrawHG hist,int color)
511: {
514: hist->color = color;
515: return 0;
516: }
518: /*@
519: PetscDrawHGSetLimits - Sets the axis limits for a histogram. If more
520: points are added after this call, the limits will be adjusted to
521: include those additional points.
523: Logically Collective on PetscDrawHG
525: Input Parameters:
526: + hist - The histogram context
527: - x_min,x_max,y_min,y_max - The limits
529: Level: intermediate
531: .seealso: PetscDrawHGCreate(), PetscDrawHGGetDraw(), PetscDrawSetSave(), PetscDrawSave(), PetscDrawHGDraw(), PetscDrawHGGetAxis()
533: @*/
534: PetscErrorCode PetscDrawHGSetLimits(PetscDrawHG hist, PetscReal x_min, PetscReal x_max, int y_min, int y_max)
535: {
538: hist->xmin = x_min;
539: hist->xmax = x_max;
540: hist->ymin = y_min;
541: hist->ymax = y_max;
542: return 0;
543: }
545: /*@
546: PetscDrawHGCalcStats - Turns on calculation of descriptive statistics
548: Not collective
550: Input Parameters:
551: + hist - The histogram context
552: - calc - Flag for calculation
554: Level: intermediate
556: .seealso: PetscDrawHGCreate(), PetscDrawHGAddValue(), PetscDrawHGView(), PetscDrawHGDraw()
558: @*/
559: PetscErrorCode PetscDrawHGCalcStats(PetscDrawHG hist, PetscBool calc)
560: {
563: hist->calcStats = calc;
564: return 0;
565: }
567: /*@
568: PetscDrawHGIntegerBins - Turns on integer width bins
570: Not collective
572: Input Parameters:
573: + hist - The histogram context
574: - ints - Flag for integer width bins
576: Level: intermediate
578: .seealso: PetscDrawHGCreate(), PetscDrawHGAddValue(), PetscDrawHGView(), PetscDrawHGDraw(), PetscDrawHGSetColor()
580: @*/
581: PetscErrorCode PetscDrawHGIntegerBins(PetscDrawHG hist, PetscBool ints)
582: {
585: hist->integerBins = ints;
586: return 0;
587: }
589: /*@C
590: PetscDrawHGGetAxis - Gets the axis context associated with a histogram.
591: This is useful if one wants to change some axis property, such as
592: labels, color, etc. The axis context should not be destroyed by the
593: application code.
595: Not Collective, PetscDrawAxis is parallel if PetscDrawHG is parallel
597: Input Parameter:
598: . hist - The histogram context
600: Output Parameter:
601: . axis - The axis context
603: Level: intermediate
605: .seealso: PetscDrawHGCreate(), PetscDrawHGAddValue(), PetscDrawHGView(), PetscDrawHGDraw(), PetscDrawHGSetColor(), PetscDrawAxis, PetscDrawHGSetLimits()
607: @*/
608: PetscErrorCode PetscDrawHGGetAxis(PetscDrawHG hist,PetscDrawAxis *axis)
609: {
612: *axis = hist->axis;
613: return 0;
614: }
616: /*@C
617: PetscDrawHGGetDraw - Gets the draw context associated with a histogram.
619: Not Collective, PetscDraw is parallel if PetscDrawHG is parallel
621: Input Parameter:
622: . hist - The histogram context
624: Output Parameter:
625: . draw - The draw context
627: Level: intermediate
629: .seealso: PetscDrawHGCreate(), PetscDrawHGAddValue(), PetscDrawHGView(), PetscDrawHGDraw(), PetscDrawHGSetColor(), PetscDrawAxis, PetscDrawHGSetLimits()
631: @*/
632: PetscErrorCode PetscDrawHGGetDraw(PetscDrawHG hist,PetscDraw *draw)
633: {
636: *draw = hist->win;
637: return 0;
638: }