Actual source code: hists.c

petsc-3.7.5 2017-01-01
Report Typos and Errors
  2: /*
  3:   Contains the data structure for plotting a histogram in a window with an axis.
  4: */
  5: #include <petscdraw.h>         /*I "petscdraw.h" I*/
  6: #include <petsc/private/petscimpl.h>         /*I "petscsys.h" I*/
  7: #include <petscviewer.h>         /*I "petscviewer.h" I*/

  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

 34: /*@C
 35:    PetscDrawHGCreate - Creates a histogram data structure.

 37:    Collective on PetscDraw

 39:    Input Parameters:
 40: +  draw  - The window where the graph will be made
 41: -  bins - The number of bins to use

 43:    Output Parameters:
 44: .  hist - The histogram context

 46:    Level: intermediate

 48:    Concepts: histogram^creating

 50: .seealso: PetscDrawHGDestroy()

 52: @*/
 53: PetscErrorCode  PetscDrawHGCreate(PetscDraw draw,int bins,PetscDrawHG *hist)
 54: {
 55:   PetscDrawHG    h;


 63:   PetscHeaderCreate(h,PETSC_DRAWHG_CLASSID,"DrawHG","Histogram","Draw",PetscObjectComm((PetscObject)draw),PetscDrawHGDestroy,NULL);
 64:   PetscLogObjectParent((PetscObject)draw,(PetscObject)h);

 66:   PetscObjectReference((PetscObject)draw);
 67:   h->win = draw;

 69:   h->view        = NULL;
 70:   h->destroy     = NULL;
 71:   h->color       = PETSC_DRAW_GREEN;
 72:   h->xmin        = PETSC_MAX_REAL;
 73:   h->xmax        = PETSC_MIN_REAL;
 74:   h->ymin        = 0.;
 75:   h->ymax        = 1.;
 76:   h->numBins     = bins;
 77:   h->maxBins     = bins;

 79:   PetscMalloc1(h->maxBins,&h->bins);

 81:   h->numValues   = 0;
 82:   h->maxValues   = CHUNKSIZE;
 83:   h->calcStats   = PETSC_FALSE;
 84:   h->integerBins = PETSC_FALSE;

 86:   PetscMalloc1(h->maxValues,&h->values);
 87:   PetscLogObjectMemory((PetscObject)h,(h->maxBins + h->maxValues)*sizeof(PetscReal));

 89:   PetscDrawAxisCreate(draw,&h->axis);
 90:   PetscLogObjectParent((PetscObject)h,(PetscObject)h->axis);

 92:   *hist = h;
 93:   return(0);
 94: }

 98: /*@
 99:    PetscDrawHGSetNumberBins - Change the number of bins that are to be drawn.

101:    Logically Collective on PetscDrawHG

103:    Input Parameter:
104: +  hist - The histogram context.
105: -  bins  - The number of bins.

107:    Level: intermediate

109:    Concepts: histogram^setting number of bins

111: @*/
112: PetscErrorCode  PetscDrawHGSetNumberBins(PetscDrawHG hist, int bins)
113: {


120:   if (hist->maxBins < bins) {
121:     PetscFree(hist->bins);
122:     PetscMalloc1(bins, &hist->bins);
123:     PetscLogObjectMemory((PetscObject)hist, (bins - hist->maxBins) * sizeof(PetscReal));
124:     hist->maxBins = bins;
125:   }
126:   hist->numBins = bins;
127:   return(0);
128: }

132: /*@
133:   PetscDrawHGReset - Clears histogram to allow for reuse with new data.

135:   Logically Collective on PetscDrawHG

137:   Input Parameter:
138: . hist - The histogram context.

140:   Level: intermediate

142:   Concepts: histogram^resetting
143: @*/
144: PetscErrorCode  PetscDrawHGReset(PetscDrawHG hist)
145: {

149:   hist->xmin      = PETSC_MAX_REAL;
150:   hist->xmax      = PETSC_MIN_REAL;
151:   hist->ymin      = 0.0;
152:   hist->ymax      = 0.0;
153:   hist->numValues = 0;
154:   return(0);
155: }

159: /*@C
160:   PetscDrawHGDestroy - Frees all space taken up by histogram data structure.

162:   Collective on PetscDrawHG

164:   Input Parameter:
165: . hist - The histogram context

167:   Level: intermediate

169: .seealso:  PetscDrawHGCreate()
170: @*/
171: PetscErrorCode  PetscDrawHGDestroy(PetscDrawHG *hist)
172: {

176:   if (!*hist) return(0);
178:   if (--((PetscObject)(*hist))->refct > 0) {*hist = NULL; return(0);}

180:   PetscFree((*hist)->bins);
181:   PetscFree((*hist)->values);
182:   PetscDrawAxisDestroy(&(*hist)->axis);
183:   PetscDrawDestroy(&(*hist)->win);
184:   PetscHeaderDestroy(hist);
185:   return(0);
186: }

190: /*@
191:   PetscDrawHGAddValue - Adds another value to the histogram.

193:   Logically Collective on PetscDrawHG

195:   Input Parameters:
196: + hist  - The histogram
197: - value - The value

199:   Level: intermediate

201:   Concepts: histogram^adding values

203: .seealso: PetscDrawHGAddValues()
204: @*/
205: PetscErrorCode  PetscDrawHGAddValue(PetscDrawHG hist, PetscReal value)
206: {

210:   /* Allocate more memory if necessary */
211:   if (hist->numValues >= hist->maxValues) {
212:     PetscReal      *tmp;

215:     PetscMalloc1(hist->maxValues+CHUNKSIZE, &tmp);
216:     PetscLogObjectMemory((PetscObject)hist, CHUNKSIZE * sizeof(PetscReal));
217:     PetscMemcpy(tmp, hist->values, hist->maxValues * sizeof(PetscReal));
218:     PetscFree(hist->values);

220:     hist->values     = tmp;
221:     hist->maxValues += CHUNKSIZE;
222:   }
223:   /* I disagree with the original Petsc implementation here. There should be no overshoot, but rather the
224:      stated convention of using half-open intervals (always the way to go) */
225:   if (!hist->numValues) {
226:     hist->xmin = value;
227:     hist->xmax = value;
228: #if 1
229:   } else {
230:     /* Update limits */
231:     if (value > hist->xmax) hist->xmax = value;
232:     if (value < hist->xmin) hist->xmin = value;
233: #else
234:   } else if (hist->numValues == 1) {
235:     /* Update limits -- We need to overshoot the largest value somewhat */
236:     if (value > hist->xmax) hist->xmax = value + 0.001*(value - hist->xmin)/hist->numBins;
237:     if (value < hist->xmin) {
238:       hist->xmin = value;
239:       hist->xmax = hist->xmax + 0.001*(hist->xmax - hist->xmin)/hist->numBins;
240:     }
241:   } else {
242:     /* Update limits -- We need to overshoot the largest value somewhat */
243:     if (value > hist->xmax) hist->xmax = value + 0.001*(hist->xmax - hist->xmin)/hist->numBins;
244:     if (value < hist->xmin) hist->xmin = value;
245: #endif
246:   }

248:   hist->values[hist->numValues++] = value;
249:   return(0);
250: }

254: /*@
255:   PetscDrawHGDraw - Redraws a histogram.

257:   Collective on PetscDrawHG

259:   Input Parameter:
260: . hist - The histogram context

262:   Level: intermediate

264: @*/
265: PetscErrorCode  PetscDrawHGDraw(PetscDrawHG hist)
266: {
267:   PetscDraw      draw;
268:   PetscBool      isnull;
269:   PetscReal      xmin,xmax,ymin,ymax,*bins,*values,binSize,binLeft,binRight,maxHeight,mean,var;
270:   char           title[256];
271:   char           xlabel[256];
272:   PetscInt       numBins,numBinsOld,numValues,initSize,i,p,bcolor,color;
273:   PetscMPIInt    rank;

278:   PetscDrawIsNull(hist->win,&isnull);
279:   if (isnull) return(0);
280:   MPI_Comm_rank(PetscObjectComm((PetscObject)hist),&rank);

282:   if ((hist->xmin >= hist->xmax) || (hist->ymin >= hist->ymax)) return(0);
283:   if (hist->numValues < 1) return(0);

285:   color = hist->color;
286:   if (color == PETSC_DRAW_ROTATE) bcolor = PETSC_DRAW_BLACK+1;
287:   else bcolor = color;

289:   xmin      = hist->xmin;
290:   xmax      = hist->xmax;
291:   ymin      = hist->ymin;
292:   ymax      = hist->ymax;
293:   numValues = hist->numValues;
294:   values    = hist->values;
295:   mean      = 0.0;
296:   var       = 0.0;

298:   draw = hist->win;
299:   PetscDrawCheckResizedWindow(draw);
300:   PetscDrawClear(draw);

302:   if (xmin == xmax) {
303:     /* Calculate number of points in each bin */
304:     bins    = hist->bins;
305:     bins[0] = 0.;
306:     for (p = 0; p < numValues; p++) {
307:       if (values[p] == xmin) bins[0]++;
308:       mean += values[p];
309:       var  += values[p]*values[p];
310:     }
311:     maxHeight = bins[0];
312:     if (maxHeight > ymax) ymax = hist->ymax = maxHeight;
313:     xmax = xmin + 1;
314:     PetscDrawAxisSetLimits(hist->axis, xmin, xmax, ymin, ymax);
315:     if (hist->calcStats) {
316:       mean /= numValues;
317:       if (numValues > 1) var = (var - numValues*mean*mean) / (numValues-1);
318:       else var = 0.0;
319:       PetscSNPrintf(title, 256, "Mean: %g  Var: %g", (double)mean, (double)var);
320:       PetscSNPrintf(xlabel,256, "Total: %D", numValues);
321:       PetscDrawAxisSetLabels(hist->axis, title, xlabel, NULL);
322:     }
323:     PetscDrawAxisDraw(hist->axis);
324:     PetscDrawCollectiveBegin(draw);
325:     if (!rank) { /* Draw bins */
326:       binLeft  = xmin;
327:       binRight = xmax;
328:       PetscDrawRectangle(draw,binLeft,ymin,binRight,bins[0],bcolor,bcolor,bcolor,bcolor);
329:       PetscDrawLine(draw,binLeft,ymin,binLeft,bins[0],PETSC_DRAW_BLACK);
330:       PetscDrawLine(draw,binRight,ymin,binRight,bins[0],PETSC_DRAW_BLACK);
331:       PetscDrawLine(draw,binLeft,bins[0],binRight,bins[0],PETSC_DRAW_BLACK);
332:     }
333:     PetscDrawCollectiveEnd(draw);
334:   } else {
335:     numBins    = hist->numBins;
336:     numBinsOld = hist->numBins;
337:     if (hist->integerBins && (((int) xmax - xmin) + 1.0e-05 > xmax - xmin)) {
338:       initSize = (int) ((int) xmax - xmin)/numBins;
339:       while (initSize*numBins != (int) xmax - xmin) {
340:         initSize = PetscMax(initSize - 1, 1);
341:         numBins  = (int) ((int) xmax - xmin)/initSize;
342:         PetscDrawHGSetNumberBins(hist, numBins);
343:       }
344:     }
345:     binSize = (xmax - xmin)/numBins;
346:     bins    = hist->bins;

348:     PetscMemzero(bins, numBins * sizeof(PetscReal));

350:     maxHeight = 0.0;
351:     for (i = 0; i < numBins; i++) {
352:       binLeft  = xmin + binSize*i;
353:       binRight = xmin + binSize*(i+1);
354:       for (p = 0; p < numValues; p++) {
355:         if ((values[p] >= binLeft) && (values[p] < binRight)) bins[i]++;
356:         /* Handle last bin separately */
357:         if ((i == numBins-1) && (values[p] == binRight)) bins[i]++;
358:         if (!i) {
359:           mean += values[p];
360:           var  += values[p]*values[p];
361:         }
362:       }
363:       maxHeight = PetscMax(maxHeight, bins[i]);
364:     }
365:     if (maxHeight > ymax) ymax = hist->ymax = maxHeight;

367:     PetscDrawAxisSetLimits(hist->axis, xmin, xmax, ymin, ymax);
368:     if (hist->calcStats) {
369:       mean /= numValues;
370:       if (numValues > 1) var = (var - numValues*mean*mean) / (numValues-1);
371:       else var = 0.0;
372:       PetscSNPrintf(title, 256,"Mean: %g  Var: %g", (double)mean, (double)var);
373:       PetscSNPrintf(xlabel,256, "Total: %D", numValues);
374:       PetscDrawAxisSetLabels(hist->axis, title, xlabel, NULL);
375:     }
376:     PetscDrawAxisDraw(hist->axis);
377:     PetscDrawCollectiveBegin(draw);
378:     if (!rank) { /* Draw bins */
379:       for (i = 0; i < numBins; i++) {
380:         binLeft  = xmin + binSize*i;
381:         binRight = xmin + binSize*(i+1);
382:         PetscDrawRectangle(draw,binLeft,ymin,binRight,bins[i],bcolor,bcolor,bcolor,bcolor);
383:         PetscDrawLine(draw,binLeft,ymin,binLeft,bins[i],PETSC_DRAW_BLACK);
384:         PetscDrawLine(draw,binRight,ymin,binRight,bins[i],PETSC_DRAW_BLACK);
385:         PetscDrawLine(draw,binLeft,bins[i],binRight,bins[i],PETSC_DRAW_BLACK);
386:         if (color == PETSC_DRAW_ROTATE && bins[i]) bcolor++;
387:         if (bcolor > PETSC_DRAW_BASIC_COLORS-1) bcolor = PETSC_DRAW_BLACK+1;
388:       }
389:     }
390:     PetscDrawCollectiveEnd(draw);
391:     PetscDrawHGSetNumberBins(hist,numBinsOld);
392:   }

394:   PetscDrawFlush(draw);
395:   PetscDrawPause(draw);
396:   return(0);
397: }

401: /*@
402:   PetscDrawHGSave - Saves a drawn image

404:   Collective on PetscDrawHG

406:   Input Parameter:
407: . hist - The histogram context

409:   Level: intermediate

411:   Concepts: histogram^saving

413: .seealso:  PetscDrawHGCreate(), PetscDrawHGGetDraw(), PetscDrawSetSave(), PetscDrawSave()
414: @*/
415: PetscErrorCode  PetscDrawHGSave(PetscDrawHG hg)
416: {

421:   PetscDrawSave(hg->win);
422:   return(0);
423: }

427: /*@
428:   PetscDrawHGView - Prints the histogram information.

430:   Not collective

432:   Input Parameter:
433: . hist - The histogram context

435:   Level: beginner

437: .keywords:  draw, histogram
438: @*/
439: PetscErrorCode  PetscDrawHGView(PetscDrawHG hist,PetscViewer viewer)
440: {
441:   PetscReal      xmax,xmin,*bins,*values,binSize,binLeft,binRight,mean,var;
443:   PetscInt       numBins,numBinsOld,numValues,initSize,i,p;


448:   if ((hist->xmin > hist->xmax) || (hist->ymin >= hist->ymax)) return(0);
449:   if (hist->numValues < 1) return(0);

451:   if (!viewer){
452:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)hist),&viewer);
453:   }
454:   PetscObjectPrintClassNamePrefixType((PetscObject)hist,viewer);
455:   xmax      = hist->xmax;
456:   xmin      = hist->xmin;
457:   numValues = hist->numValues;
458:   values    = hist->values;
459:   mean      = 0.0;
460:   var       = 0.0;
461:   if (xmax == xmin) {
462:     /* Calculate number of points in the bin */
463:     bins    = hist->bins;
464:     bins[0] = 0.;
465:     for (p = 0; p < numValues; p++) {
466:       if (values[p] == xmin) bins[0]++;
467:       mean += values[p];
468:       var  += values[p]*values[p];
469:     }
470:     /* Draw bins */
471:     PetscViewerASCIIPrintf(viewer, "Bin %2d (%6.2g - %6.2g): %.0g\n", 0, (double)xmin, (double)xmax, (double)bins[0]);
472:   } else {
473:     numBins    = hist->numBins;
474:     numBinsOld = hist->numBins;
475:     if (hist->integerBins && (((int) xmax - xmin) + 1.0e-05 > xmax - xmin)) {
476:       initSize = (int) ((int) xmax - xmin)/numBins;
477:       while (initSize*numBins != (int) xmax - xmin) {
478:         initSize = PetscMax(initSize - 1, 1);
479:         numBins  = (int) ((int) xmax - xmin)/initSize;
480:         PetscDrawHGSetNumberBins(hist, numBins);
481:       }
482:     }
483:     binSize = (xmax - xmin)/numBins;
484:     bins    = hist->bins;

486:     /* Calculate number of points in each bin */
487:     PetscMemzero(bins, numBins * sizeof(PetscReal));
488:     for (i = 0; i < numBins; i++) {
489:       binLeft  = xmin + binSize*i;
490:       binRight = xmin + binSize*(i+1);
491:       for (p = 0; p < numValues; p++) {
492:         if ((values[p] >= binLeft) && (values[p] < binRight)) bins[i]++;
493:         /* Handle last bin separately */
494:         if ((i == numBins-1) && (values[p] == binRight)) bins[i]++;
495:         if (!i) {
496:           mean += values[p];
497:           var  += values[p]*values[p];
498:         }
499:       }
500:     }
501:     /* Draw bins */
502:     for (i = 0; i < numBins; i++) {
503:       binLeft  = xmin + binSize*i;
504:       binRight = xmin + binSize*(i+1);
505:       PetscViewerASCIIPrintf(viewer, "Bin %2d (%6.2g - %6.2g): %.0g\n", (int)i, (double)binLeft, (double)binRight, (double)bins[i]);
506:     }
507:     PetscDrawHGSetNumberBins(hist, numBinsOld);
508:   }

510:   if (hist->calcStats) {
511:     mean /= numValues;
512:     if (numValues > 1) var = (var - numValues*mean*mean) / (numValues-1);
513:     else var = 0.0;
514:     PetscViewerASCIIPrintf(viewer, "Mean: %g  Var: %g\n", (double)mean, (double)var);
515:     PetscViewerASCIIPrintf(viewer, "Total: %D\n", numValues);
516:   }
517:   return(0);
518: }

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: @*/
535: PetscErrorCode  PetscDrawHGSetColor(PetscDrawHG hist,int color)
536: {

540:   hist->color = color;
541:   return(0);
542: }

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
560: @*/
561: PetscErrorCode  PetscDrawHGSetLimits(PetscDrawHG hist, PetscReal x_min, PetscReal x_max, int y_min, int y_max)
562: {

566:   hist->xmin = x_min;
567:   hist->xmax = x_max;
568:   hist->ymin = y_min;
569:   hist->ymax = y_max;
570:   return(0);
571: }

575: /*@
576:   PetscDrawHGCalcStats - Turns on calculation of descriptive statistics

578:   Not collective

580:   Input Parameters:
581: + hist - The histogram context
582: - calc - Flag for calculation

584:   Level: intermediate

586: .keywords:  draw, histogram, statistics

588: @*/
589: PetscErrorCode  PetscDrawHGCalcStats(PetscDrawHG hist, PetscBool calc)
590: {

594:   hist->calcStats = calc;
595:   return(0);
596: }

600: /*@
601:   PetscDrawHGIntegerBins - Turns on integer width bins

603:   Not collective

605:   Input Parameters:
606: + hist - The histogram context
607: - ints - Flag for integer width bins

609:   Level: intermediate

611: .keywords:  draw, histogram, statistics
612: @*/
613: PetscErrorCode  PetscDrawHGIntegerBins(PetscDrawHG hist, PetscBool ints)
614: {

618:   hist->integerBins = ints;
619:   return(0);
620: }

624: /*@C
625:   PetscDrawHGGetAxis - Gets the axis context associated with a histogram.
626:   This is useful if one wants to change some axis property, such as
627:   labels, color, etc. The axis context should not be destroyed by the
628:   application code.

630:   Not Collective, PetscDrawAxis is parallel if PetscDrawHG is parallel

632:   Input Parameter:
633: . hist - The histogram context

635:   Output Parameter:
636: . axis - The axis context

638:   Level: intermediate

640: @*/
641: PetscErrorCode  PetscDrawHGGetAxis(PetscDrawHG hist,PetscDrawAxis *axis)
642: {
646:   *axis = hist->axis;
647:   return(0);
648: }

652: /*@C
653:   PetscDrawHGGetDraw - Gets the draw context associated with a histogram.

655:   Not Collective, PetscDraw is parallel if PetscDrawHG is parallel

657:   Input Parameter:
658: . hist - The histogram context

660:   Output Parameter:
661: . draw  - The draw context

663:   Level: intermediate

665: @*/
666: PetscErrorCode  PetscDrawHGGetDraw(PetscDrawHG hist,PetscDraw *draw)
667: {
671:   *draw = hist->win;
672:   return(0);
673: }