Actual source code: hists.c
petsc-3.9.4 2018-09-11
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: The difference between a bar chart, PetscDrawBar, and a histogram, PetscDrawHG, is explained here http://stattrek.com/statistics/charts/histogram.aspx?Tutorial=AP
46: The histogram is only displayed when PetscDrawHGDraw() is called.
48: 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
49: zeroth MPI process in the communicator. All MPI processes in the communicator must call PetscDrawHGDraw() to display the updated graph.
51: Level: intermediate
53: Concepts: histogram^creating
55: .seealso: PetscDrawHGDestroy(), PetscDrawHG, PetscDrawBarCreate(), PetscDrawBar, PetscDrawLGCreate(), PetscDrawLG, PetscDrawSPCreate(), PetscDrawSP,
56: PetscDrawHGSetNumberBins(), PetscDrawHGReset(), PetscDrawHGAddValue(), PetscDrawHGDraw(), PetscDrawHGSave(), PetscDrawHGView(), PetscDrawHGSetColor(),
57: PetscDrawHGSetLimits(), PetscDrawHGCalcStats(), PetscDrawHGIntegerBins(), PetscDrawHGGetAxis(), PetscDrawAxis, PetscDrawHGGetDraw()
59: @*/
60: PetscErrorCode PetscDrawHGCreate(PetscDraw draw,int bins,PetscDrawHG *hist)
61: {
62: PetscDrawHG h;
70: PetscHeaderCreate(h,PETSC_DRAWHG_CLASSID,"DrawHG","Histogram","Draw",PetscObjectComm((PetscObject)draw),PetscDrawHGDestroy,NULL);
71: PetscLogObjectParent((PetscObject)draw,(PetscObject)h);
73: PetscObjectReference((PetscObject)draw);
74: h->win = draw;
76: h->view = NULL;
77: h->destroy = NULL;
78: h->color = PETSC_DRAW_GREEN;
79: h->xmin = PETSC_MAX_REAL;
80: h->xmax = PETSC_MIN_REAL;
81: h->ymin = 0.;
82: h->ymax = 1.;
83: h->numBins = bins;
84: h->maxBins = bins;
86: PetscMalloc1(h->maxBins,&h->bins);
88: h->numValues = 0;
89: h->maxValues = CHUNKSIZE;
90: h->calcStats = PETSC_FALSE;
91: h->integerBins = PETSC_FALSE;
93: PetscMalloc1(h->maxValues,&h->values);
94: PetscLogObjectMemory((PetscObject)h,(h->maxBins + h->maxValues)*sizeof(PetscReal));
96: PetscDrawAxisCreate(draw,&h->axis);
97: PetscLogObjectParent((PetscObject)h,(PetscObject)h->axis);
99: *hist = h;
100: return(0);
101: }
103: /*@
104: PetscDrawHGSetNumberBins - Change the number of bins that are to be drawn.
106: Logically Collective on PetscDrawHG
108: Input Parameter:
109: + hist - The histogram context.
110: - bins - The number of bins.
112: Level: intermediate
114: Concepts: histogram^setting number of bins
116: .seealso: PetscDrawHGCreate(), PetscDrawHG, PetscDrawHGDraw(), PetscDrawHGIntegerBins()
118: @*/
119: PetscErrorCode PetscDrawHGSetNumberBins(PetscDrawHG hist, int bins)
120: {
127: if (hist->maxBins < bins) {
128: PetscFree(hist->bins);
129: PetscMalloc1(bins, &hist->bins);
130: PetscLogObjectMemory((PetscObject)hist, (bins - hist->maxBins) * sizeof(PetscReal));
131: hist->maxBins = bins;
132: }
133: hist->numBins = bins;
134: return(0);
135: }
137: /*@
138: PetscDrawHGReset - Clears histogram to allow for reuse with new data.
140: Logically Collective on PetscDrawHG
142: Input Parameter:
143: . hist - The histogram context.
145: Level: intermediate
147: Concepts: histogram^resetting
149: .seealso: PetscDrawHGCreate(), PetscDrawHG, PetscDrawHGDraw(), PetscDrawHGAddValue()
151: @*/
152: PetscErrorCode PetscDrawHGReset(PetscDrawHG hist)
153: {
157: hist->xmin = PETSC_MAX_REAL;
158: hist->xmax = PETSC_MIN_REAL;
159: hist->ymin = 0.0;
160: hist->ymax = 0.0;
161: hist->numValues = 0;
162: return(0);
163: }
165: /*@C
166: PetscDrawHGDestroy - Frees all space taken up by histogram data structure.
168: Collective on PetscDrawHG
170: Input Parameter:
171: . hist - The histogram context
173: Level: intermediate
175: .seealso: PetscDrawHGCreate(), PetscDrawHG
176: @*/
177: PetscErrorCode PetscDrawHGDestroy(PetscDrawHG *hist)
178: {
182: if (!*hist) return(0);
184: if (--((PetscObject)(*hist))->refct > 0) {*hist = NULL; return(0);}
186: PetscFree((*hist)->bins);
187: PetscFree((*hist)->values);
188: PetscDrawAxisDestroy(&(*hist)->axis);
189: PetscDrawDestroy(&(*hist)->win);
190: PetscHeaderDestroy(hist);
191: return(0);
192: }
194: /*@
195: PetscDrawHGAddValue - Adds another value to the histogram.
197: Logically Collective on PetscDrawHG
199: Input Parameters:
200: + hist - The histogram
201: - value - The value
203: Level: intermediate
205: Concepts: histogram^adding values
207: .seealso: PetscDrawHGCreate(), PetscDrawHG, PetscDrawHGDraw(), PetscDrawHGAddValue(), PetscDrawHGReset()
208: @*/
209: PetscErrorCode PetscDrawHGAddValue(PetscDrawHG hist, PetscReal value)
210: {
214: /* Allocate more memory if necessary */
215: if (hist->numValues >= hist->maxValues) {
216: PetscReal *tmp;
219: PetscMalloc1(hist->maxValues+CHUNKSIZE, &tmp);
220: PetscLogObjectMemory((PetscObject)hist, CHUNKSIZE * sizeof(PetscReal));
221: PetscMemcpy(tmp, hist->values, hist->maxValues * sizeof(PetscReal));
222: PetscFree(hist->values);
224: hist->values = tmp;
225: hist->maxValues += CHUNKSIZE;
226: }
227: /* I disagree with the original Petsc implementation here. There should be no overshoot, but rather the
228: stated convention of using half-open intervals (always the way to go) */
229: if (!hist->numValues) {
230: hist->xmin = value;
231: hist->xmax = value;
232: #if 1
233: } else {
234: /* Update limits */
235: if (value > hist->xmax) hist->xmax = value;
236: if (value < hist->xmin) hist->xmin = value;
237: #else
238: } else if (hist->numValues == 1) {
239: /* Update limits -- We need to overshoot the largest value somewhat */
240: if (value > hist->xmax) hist->xmax = value + 0.001*(value - hist->xmin)/hist->numBins;
241: if (value < hist->xmin) {
242: hist->xmin = value;
243: hist->xmax = hist->xmax + 0.001*(hist->xmax - hist->xmin)/hist->numBins;
244: }
245: } else {
246: /* Update limits -- We need to overshoot the largest value somewhat */
247: if (value > hist->xmax) hist->xmax = value + 0.001*(hist->xmax - hist->xmin)/hist->numBins;
248: if (value < hist->xmin) hist->xmin = value;
249: #endif
250: }
252: hist->values[hist->numValues++] = value;
253: return(0);
254: }
256: /*@
257: PetscDrawHGDraw - Redraws a histogram.
259: Collective on PetscDrawHG
261: Input Parameter:
262: . hist - The histogram context
264: Level: intermediate
266: .seealso: PetscDrawHGCreate(), PetscDrawHG, PetscDrawHGDraw(), PetscDrawHGAddValue(), PetscDrawHGReset()
268: @*/
269: PetscErrorCode PetscDrawHGDraw(PetscDrawHG hist)
270: {
271: PetscDraw draw;
272: PetscBool isnull;
273: PetscReal xmin,xmax,ymin,ymax,*bins,*values,binSize,binLeft,binRight,maxHeight,mean,var;
274: char title[256];
275: char xlabel[256];
276: PetscInt numBins,numBinsOld,numValues,initSize,i,p,bcolor,color;
277: PetscMPIInt rank;
282: PetscDrawIsNull(hist->win,&isnull);
283: if (isnull) return(0);
284: MPI_Comm_rank(PetscObjectComm((PetscObject)hist),&rank);
286: if ((hist->xmin >= hist->xmax) || (hist->ymin >= hist->ymax)) return(0);
287: if (hist->numValues < 1) return(0);
289: color = hist->color;
290: if (color == PETSC_DRAW_ROTATE) bcolor = PETSC_DRAW_BLACK+1;
291: else bcolor = color;
293: xmin = hist->xmin;
294: xmax = hist->xmax;
295: ymin = hist->ymin;
296: ymax = hist->ymax;
297: numValues = hist->numValues;
298: values = hist->values;
299: mean = 0.0;
300: var = 0.0;
302: draw = hist->win;
303: PetscDrawCheckResizedWindow(draw);
304: PetscDrawClear(draw);
306: if (xmin == xmax) {
307: /* Calculate number of points in each bin */
308: bins = hist->bins;
309: bins[0] = 0.;
310: for (p = 0; p < numValues; p++) {
311: if (values[p] == xmin) bins[0]++;
312: mean += values[p];
313: var += values[p]*values[p];
314: }
315: maxHeight = bins[0];
316: if (maxHeight > ymax) ymax = hist->ymax = maxHeight;
317: xmax = xmin + 1;
318: PetscDrawAxisSetLimits(hist->axis, xmin, xmax, ymin, ymax);
319: if (hist->calcStats) {
320: mean /= numValues;
321: if (numValues > 1) var = (var - numValues*mean*mean) / (numValues-1);
322: else var = 0.0;
323: PetscSNPrintf(title, 256, "Mean: %g Var: %g", (double)mean, (double)var);
324: PetscSNPrintf(xlabel,256, "Total: %D", numValues);
325: PetscDrawAxisSetLabels(hist->axis, title, xlabel, NULL);
326: }
327: PetscDrawAxisDraw(hist->axis);
328: PetscDrawCollectiveBegin(draw);
329: if (!rank) { /* Draw bins */
330: binLeft = xmin;
331: binRight = xmax;
332: PetscDrawRectangle(draw,binLeft,ymin,binRight,bins[0],bcolor,bcolor,bcolor,bcolor);
333: PetscDrawLine(draw,binLeft,ymin,binLeft,bins[0],PETSC_DRAW_BLACK);
334: PetscDrawLine(draw,binRight,ymin,binRight,bins[0],PETSC_DRAW_BLACK);
335: PetscDrawLine(draw,binLeft,bins[0],binRight,bins[0],PETSC_DRAW_BLACK);
336: }
337: PetscDrawCollectiveEnd(draw);
338: } else {
339: numBins = hist->numBins;
340: numBinsOld = hist->numBins;
341: if (hist->integerBins && (((int) xmax - xmin) + 1.0e-05 > xmax - xmin)) {
342: initSize = (int) ((int) xmax - xmin)/numBins;
343: while (initSize*numBins != (int) xmax - xmin) {
344: initSize = PetscMax(initSize - 1, 1);
345: numBins = (int) ((int) xmax - xmin)/initSize;
346: PetscDrawHGSetNumberBins(hist, numBins);
347: }
348: }
349: binSize = (xmax - xmin)/numBins;
350: bins = hist->bins;
352: PetscMemzero(bins, numBins * sizeof(PetscReal));
354: maxHeight = 0.0;
355: for (i = 0; i < numBins; i++) {
356: binLeft = xmin + binSize*i;
357: binRight = xmin + binSize*(i+1);
358: for (p = 0; p < numValues; p++) {
359: if ((values[p] >= binLeft) && (values[p] < binRight)) bins[i]++;
360: /* Handle last bin separately */
361: if ((i == numBins-1) && (values[p] == binRight)) bins[i]++;
362: if (!i) {
363: mean += values[p];
364: var += values[p]*values[p];
365: }
366: }
367: maxHeight = PetscMax(maxHeight, bins[i]);
368: }
369: if (maxHeight > ymax) ymax = hist->ymax = maxHeight;
371: PetscDrawAxisSetLimits(hist->axis, xmin, xmax, ymin, ymax);
372: if (hist->calcStats) {
373: mean /= numValues;
374: if (numValues > 1) var = (var - numValues*mean*mean) / (numValues-1);
375: else var = 0.0;
376: PetscSNPrintf(title, 256,"Mean: %g Var: %g", (double)mean, (double)var);
377: PetscSNPrintf(xlabel,256, "Total: %D", numValues);
378: PetscDrawAxisSetLabels(hist->axis, title, xlabel, NULL);
379: }
380: PetscDrawAxisDraw(hist->axis);
381: PetscDrawCollectiveBegin(draw);
382: if (!rank) { /* Draw bins */
383: for (i = 0; i < numBins; i++) {
384: binLeft = xmin + binSize*i;
385: binRight = xmin + binSize*(i+1);
386: PetscDrawRectangle(draw,binLeft,ymin,binRight,bins[i],bcolor,bcolor,bcolor,bcolor);
387: PetscDrawLine(draw,binLeft,ymin,binLeft,bins[i],PETSC_DRAW_BLACK);
388: PetscDrawLine(draw,binRight,ymin,binRight,bins[i],PETSC_DRAW_BLACK);
389: PetscDrawLine(draw,binLeft,bins[i],binRight,bins[i],PETSC_DRAW_BLACK);
390: if (color == PETSC_DRAW_ROTATE && bins[i]) bcolor++;
391: if (bcolor > PETSC_DRAW_BASIC_COLORS-1) bcolor = PETSC_DRAW_BLACK+1;
392: }
393: }
394: PetscDrawCollectiveEnd(draw);
395: PetscDrawHGSetNumberBins(hist,numBinsOld);
396: }
398: PetscDrawFlush(draw);
399: PetscDrawPause(draw);
400: return(0);
401: }
403: /*@
404: PetscDrawHGSave - Saves a drawn image
406: Collective on PetscDrawHG
408: Input Parameter:
409: . hist - The histogram context
411: Level: intermediate
413: Concepts: histogram^saving
415: .seealso: PetscDrawHGCreate(), PetscDrawHGGetDraw(), PetscDrawSetSave(), PetscDrawSave(), PetscDrawHGDraw()
416: @*/
417: PetscErrorCode PetscDrawHGSave(PetscDrawHG hg)
418: {
423: PetscDrawSave(hg->win);
424: return(0);
425: }
427: /*@
428: PetscDrawHGView - Prints the histogram information.
430: Not collective
432: Input Parameter:
433: . hist - The histogram context
435: Level: beginner
437: .seealso: PetscDrawHGCreate(), PetscDrawHGGetDraw(), PetscDrawSetSave(), PetscDrawSave(), PetscDrawHGDraw()
439: .keywords: draw, histogram
440: @*/
441: PetscErrorCode PetscDrawHGView(PetscDrawHG hist,PetscViewer viewer)
442: {
443: PetscReal xmax,xmin,*bins,*values,binSize,binLeft,binRight,mean,var;
445: PetscInt numBins,numBinsOld,numValues,initSize,i,p;
450: if ((hist->xmin > hist->xmax) || (hist->ymin >= hist->ymax)) return(0);
451: if (hist->numValues < 1) return(0);
453: if (!viewer){
454: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)hist),&viewer);
455: }
456: PetscObjectPrintClassNamePrefixType((PetscObject)hist,viewer);
457: xmax = hist->xmax;
458: xmin = hist->xmin;
459: numValues = hist->numValues;
460: values = hist->values;
461: mean = 0.0;
462: var = 0.0;
463: if (xmax == xmin) {
464: /* Calculate number of points in the bin */
465: bins = hist->bins;
466: bins[0] = 0.;
467: for (p = 0; p < numValues; p++) {
468: if (values[p] == xmin) bins[0]++;
469: mean += values[p];
470: var += values[p]*values[p];
471: }
472: /* Draw bins */
473: PetscViewerASCIIPrintf(viewer, "Bin %2d (%6.2g - %6.2g): %.0g\n", 0, (double)xmin, (double)xmax, (double)bins[0]);
474: } else {
475: numBins = hist->numBins;
476: numBinsOld = hist->numBins;
477: if (hist->integerBins && (((int) xmax - xmin) + 1.0e-05 > xmax - xmin)) {
478: initSize = (int) ((int) xmax - xmin)/numBins;
479: while (initSize*numBins != (int) xmax - xmin) {
480: initSize = PetscMax(initSize - 1, 1);
481: numBins = (int) ((int) xmax - xmin)/initSize;
482: PetscDrawHGSetNumberBins(hist, numBins);
483: }
484: }
485: binSize = (xmax - xmin)/numBins;
486: bins = hist->bins;
488: /* Calculate number of points in each bin */
489: PetscMemzero(bins, numBins * sizeof(PetscReal));
490: for (i = 0; i < numBins; i++) {
491: binLeft = xmin + binSize*i;
492: binRight = xmin + binSize*(i+1);
493: for (p = 0; p < numValues; p++) {
494: if ((values[p] >= binLeft) && (values[p] < binRight)) bins[i]++;
495: /* Handle last bin separately */
496: if ((i == numBins-1) && (values[p] == binRight)) bins[i]++;
497: if (!i) {
498: mean += values[p];
499: var += values[p]*values[p];
500: }
501: }
502: }
503: /* Draw bins */
504: for (i = 0; i < numBins; i++) {
505: binLeft = xmin + binSize*i;
506: binRight = xmin + binSize*(i+1);
507: PetscViewerASCIIPrintf(viewer, "Bin %2d (%6.2g - %6.2g): %.0g\n", (int)i, (double)binLeft, (double)binRight, (double)bins[i]);
508: }
509: PetscDrawHGSetNumberBins(hist, numBinsOld);
510: }
512: if (hist->calcStats) {
513: mean /= numValues;
514: if (numValues > 1) var = (var - numValues*mean*mean) / (numValues-1);
515: else var = 0.0;
516: PetscViewerASCIIPrintf(viewer, "Mean: %g Var: %g\n", (double)mean, (double)var);
517: PetscViewerASCIIPrintf(viewer, "Total: %D\n", numValues);
518: }
519: return(0);
520: }
522: /*@
523: PetscDrawHGSetColor - Sets the color the bars will be drawn with.
525: Logically Collective on PetscDrawHG
527: Input Parameters:
528: + hist - The histogram context
529: - color - one of the colors defined in petscdraw.h or PETSC_DRAW_ROTATE to make each bar a
530: different color
532: Level: intermediate
534: .seealso: PetscDrawHGCreate(), PetscDrawHGGetDraw(), PetscDrawSetSave(), PetscDrawSave(), PetscDrawHGDraw(), PetscDrawHGGetAxis()
536: @*/
537: PetscErrorCode PetscDrawHGSetColor(PetscDrawHG hist,int color)
538: {
542: hist->color = color;
543: return(0);
544: }
546: /*@
547: PetscDrawHGSetLimits - Sets the axis limits for a histogram. If more
548: points are added after this call, the limits will be adjusted to
549: include those additional points.
551: Logically Collective on PetscDrawHG
553: Input Parameters:
554: + hist - The histogram context
555: - x_min,x_max,y_min,y_max - The limits
557: Level: intermediate
559: Concepts: histogram^setting axis
561: .seealso: PetscDrawHGCreate(), PetscDrawHGGetDraw(), PetscDrawSetSave(), PetscDrawSave(), PetscDrawHGDraw(), PetscDrawHGGetAxis()
563: @*/
564: PetscErrorCode PetscDrawHGSetLimits(PetscDrawHG hist, PetscReal x_min, PetscReal x_max, int y_min, int y_max)
565: {
569: hist->xmin = x_min;
570: hist->xmax = x_max;
571: hist->ymin = y_min;
572: hist->ymax = y_max;
573: return(0);
574: }
576: /*@
577: PetscDrawHGCalcStats - Turns on calculation of descriptive statistics
579: Not collective
581: Input Parameters:
582: + hist - The histogram context
583: - calc - Flag for calculation
585: Level: intermediate
587: .keywords: draw, histogram, statistics
589: .seealso: PetscDrawHGCreate(), PetscDrawHGAddValue(), PetscDrawHGView(), PetscDrawHGDraw()
591: @*/
592: PetscErrorCode PetscDrawHGCalcStats(PetscDrawHG hist, PetscBool calc)
593: {
597: hist->calcStats = calc;
598: return(0);
599: }
601: /*@
602: PetscDrawHGIntegerBins - Turns on integer width bins
604: Not collective
606: Input Parameters:
607: + hist - The histogram context
608: - ints - Flag for integer width bins
610: Level: intermediate
612: .keywords: draw, histogram, statistics
614: .seealso: PetscDrawHGCreate(), PetscDrawHGAddValue(), PetscDrawHGView(), PetscDrawHGDraw(), PetscDrawHGSetColor()
616: @*/
617: PetscErrorCode PetscDrawHGIntegerBins(PetscDrawHG hist, PetscBool ints)
618: {
622: hist->integerBins = ints;
623: return(0);
624: }
626: /*@C
627: PetscDrawHGGetAxis - Gets the axis context associated with a histogram.
628: This is useful if one wants to change some axis property, such as
629: labels, color, etc. The axis context should not be destroyed by the
630: application code.
632: Not Collective, PetscDrawAxis is parallel if PetscDrawHG is parallel
634: Input Parameter:
635: . hist - The histogram context
637: Output Parameter:
638: . axis - The axis context
640: Level: intermediate
642: .seealso: PetscDrawHGCreate(), PetscDrawHGAddValue(), PetscDrawHGView(), PetscDrawHGDraw(), PetscDrawHGSetColor(), PetscDrawAxis, PetscDrawHGSetLimits()
644: @*/
645: PetscErrorCode PetscDrawHGGetAxis(PetscDrawHG hist,PetscDrawAxis *axis)
646: {
650: *axis = hist->axis;
651: return(0);
652: }
654: /*@C
655: PetscDrawHGGetDraw - Gets the draw context associated with a histogram.
657: Not Collective, PetscDraw is parallel if PetscDrawHG is parallel
659: Input Parameter:
660: . hist - The histogram context
662: Output Parameter:
663: . draw - The draw context
665: Level: intermediate
667: .seealso: PetscDrawHGCreate(), PetscDrawHGAddValue(), PetscDrawHGView(), PetscDrawHGDraw(), PetscDrawHGSetColor(), PetscDrawAxis, PetscDrawHGSetLimits()
669: @*/
670: PetscErrorCode PetscDrawHGGetDraw(PetscDrawHG hist,PetscDraw *draw)
671: {
675: *draw = hist->win;
676: return(0);
677: }