Actual source code: hists.c

petsc-3.6.1 2015-07-22
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 over 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;
 56:   PetscBool      isnull;

 62:   PetscHeaderCreate(h, PETSC_DRAWHG_CLASSID,  "PetscDrawHG", "Histogram", "Draw", PetscObjectComm((PetscObject)draw), PetscDrawHGDestroy, NULL);

 64:   h->view        = NULL;
 65:   h->destroy     = NULL;
 66:   h->win         = draw;

 68:   PetscObjectReference((PetscObject) draw);

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

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

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

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

 89:   if (!isnull) {
 90:     PetscDrawAxisCreate(draw, &h->axis);
 91:     PetscLogObjectParent((PetscObject)h, (PetscObject)h->axis);
 92:   } else h->axis = NULL;
 93:   *hist = h;
 94:   return(0);
 95: }

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

102:    Not Collective (ignored except on processor 0 of PetscDrawHG)

104:    Input Parameter:
105: +  hist - The histogram context.
106: -  dim  - The number of curves.

108:    Level: intermediate

110:    Concepts: histogram^setting number of bins

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

119:   if (hist->maxBins < bins) {
120:     PetscFree(hist->bins);
121:     PetscMalloc1(bins, &hist->bins);
122:     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:   Not Collective (ignored except on processor 0 of PetscDrawHG)

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

140:   Level: intermediate

142:   Concepts: histogram^resetting
143: @*/
144: PetscErrorCode  PetscDrawHGReset(PetscDrawHG hist)
145: {
148:   hist->xmin      = PETSC_MAX_REAL;
149:   hist->xmax      = PETSC_MIN_REAL;
150:   hist->ymin      = 0.0;
151:   hist->ymax      = 0.0;
152:   hist->numValues = 0;
153:   return(0);
154: }

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

161:   Collective over PetscDrawHG

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

166:   Level: intermediate

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

175:   if (!*hist) return(0);

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

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

192:   Not Collective (ignored except on processor 0 of PetscDrawHG)

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

198:   Level: intermediate

200:   Concepts: histogram^adding values

202: .seealso: PetscDrawHGAddValues()
203: @*/
204: PetscErrorCode  PetscDrawHGAddValue(PetscDrawHG hist, PetscReal value)
205: {
208:   /* Allocate more memory if necessary */
209:   if (hist->numValues >= hist->maxValues) {
210:     PetscReal      *tmp;

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

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

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

252: /*@
253:   PetscDrawHGDraw - Redraws a histogram.

255:   Not Collective (ignored except on processor 0 of PetscDrawHG)

257:   Input Parameter:
258: . hist - The histogram context

260:   Level: intermediate

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

275:   PetscObjectTypeCompare((PetscObject) draw, PETSC_DRAW_NULL, &isnull);
276:   if (isnull) return(0);
277:   if ((hist->xmin >= hist->xmax) || (hist->ymin >= hist->ymax)) return(0);
278:   if (hist->numValues < 1) return(0);

280: #if 0
281:   MPI_Comm_rank(PetscObjectComm((PetscObject)hist),&rank);
282:   if (rank) return(0);
283: #endif

285:   color = hist->color;
286:   if (color == PETSC_DRAW_ROTATE) bcolor = 2;
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:   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:     /* Draw bins */
322:     binLeft  = xmin;
323:     binRight = xmax;
324:     PetscDrawRectangle(draw,binLeft,ymin,binRight,bins[0],bcolor,bcolor,bcolor,bcolor);

326:     if (color == PETSC_DRAW_ROTATE && bins[0] != 0.0) bcolor++;
327:     if (bcolor > 31) bcolor = 2;

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:   } else {
333:     numBins    = hist->numBins;
334:     numBinsOld = hist->numBins;
335:     if (hist->integerBins && (((int) xmax - xmin) + 1.0e-05 > xmax - xmin)) {
336:       initSize = (int) ((int) xmax - xmin)/numBins;
337:       while (initSize*numBins != (int) xmax - xmin) {
338:         initSize = PetscMax(initSize - 1, 1);
339:         numBins  = (int) ((int) xmax - xmin)/initSize;
340:         PetscDrawHGSetNumberBins(hist, numBins);
341:       }
342:     }
343:     binSize = (xmax - xmin)/numBins;
344:     bins    = hist->bins;

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

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

365:     PetscDrawAxisSetLimits(hist->axis, xmin, xmax, ymin, ymax);
366:     if (hist->calcStats) {
367:       mean /= numValues;
368:       if (numValues > 1) var = (var - numValues*mean*mean) / (numValues-1);
369:       else var = 0.0;

371:       PetscSNPrintf(title, 256,"Mean: %g  Var: %g", (double)mean, (double)var);
372:       PetscSNPrintf(xlabel,256, "Total: %D", numValues);
373:       PetscDrawAxisSetLabels(hist->axis, title, xlabel, NULL);
374:     }
375:     PetscDrawAxisDraw(hist->axis);
376:     /* Draw bins */
377:     for (i = 0; i < numBins; i++) {
378:       binLeft  = xmin + binSize*i;
379:       binRight = xmin + binSize*(i+1);
380:       PetscDrawRectangle(draw,binLeft,ymin,binRight,bins[i],bcolor,bcolor,bcolor,bcolor);
381:       if (color == PETSC_DRAW_ROTATE && bins[i]) bcolor++;
382:       if (bcolor > 31) bcolor = 2;
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:     }
387:     PetscDrawHGSetNumberBins(hist, numBinsOld);
388:   }
389:   PetscDrawSynchronizedFlush(draw);
390:   PetscDrawPause(draw);
391:   return(0);
392: }

396: /*@
397:   PetscDrawHGView - Prints the histogram information.

399:   Not collective

401:   Input Parameter:
402: . hist - The histogram context

404:   Level: beginner

406: .keywords:  draw, histogram
407: @*/
408: PetscErrorCode  PetscDrawHGView(PetscDrawHG hist,PetscViewer viewer)
409: {
410:   PetscReal      xmax,xmin,*bins,*values,binSize,binLeft,binRight,mean,var;
412:   PetscInt       numBins,numBinsOld,numValues,initSize,i,p;

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

419:   if (!viewer){
420:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)hist),&viewer);
421:   }
422:   PetscObjectPrintClassNamePrefixType((PetscObject)hist,viewer);
423:   xmax      = hist->xmax;
424:   xmin      = hist->xmin;
425:   numValues = hist->numValues;
426:   values    = hist->values;
427:   mean      = 0.0;
428:   var       = 0.0;
429:   if (xmax == xmin) {
430:     /* Calculate number of points in the bin */
431:     bins    = hist->bins;
432:     bins[0] = 0.;
433:     for (p = 0; p < numValues; p++) {
434:       if (values[p] == xmin) bins[0]++;
435:       mean += values[p];
436:       var  += values[p]*values[p];
437:     }
438:     /* Draw bins */
439:     PetscViewerASCIIPrintf(viewer, "Bin %2d (%6.2g - %6.2g): %.0g\n", 0, (double)xmin, (double)xmax, (double)bins[0]);
440:   } else {
441:     numBins    = hist->numBins;
442:     numBinsOld = hist->numBins;
443:     if (hist->integerBins && (((int) xmax - xmin) + 1.0e-05 > xmax - xmin)) {
444:       initSize = (int) ((int) xmax - xmin)/numBins;
445:       while (initSize*numBins != (int) xmax - xmin) {
446:         initSize = PetscMax(initSize - 1, 1);
447:         numBins  = (int) ((int) xmax - xmin)/initSize;
448:         PetscDrawHGSetNumberBins(hist, numBins);
449:       }
450:     }
451:     binSize = (xmax - xmin)/numBins;
452:     bins    = hist->bins;

454:     /* Calculate number of points in each bin */
455:     PetscMemzero(bins, numBins * sizeof(PetscReal));
456:     for (i = 0; i < numBins; i++) {
457:       binLeft  = xmin + binSize*i;
458:       binRight = xmin + binSize*(i+1);
459:       for (p = 0; p < numValues; p++) {
460:         if ((values[p] >= binLeft) && (values[p] < binRight)) bins[i]++;
461:         /* Handle last bin separately */
462:         if ((i == numBins-1) && (values[p] == binRight)) bins[i]++;
463:         if (!i) {
464:           mean += values[p];
465:           var  += values[p]*values[p];
466:         }
467:       }
468:     }
469:     /* Draw bins */
470:     for (i = 0; i < numBins; i++) {
471:       binLeft  = xmin + binSize*i;
472:       binRight = xmin + binSize*(i+1);
473:       PetscViewerASCIIPrintf(viewer, "Bin %2d (%6.2g - %6.2g): %.0g\n", (int)i, (double)binLeft, (double)binRight, (double)bins[i]);
474:     }
475:     PetscDrawHGSetNumberBins(hist, numBinsOld);
476:   }

478:   if (hist->calcStats) {
479:     mean /= numValues;
480:     if (numValues > 1) var = (var - numValues*mean*mean) / (numValues-1);
481:     else var = 0.0;
482:     PetscViewerASCIIPrintf(viewer, "Mean: %g  Var: %g\n", (double)mean, (double)var);
483:     PetscViewerASCIIPrintf(viewer, "Total: %D\n", numValues);
484:   }
485:   return(0);
486: }

490: /*@
491:   PetscDrawHGSetColor - Sets the color the bars will be drawn with.

493:   Not Collective (ignored except on processor 0 of PetscDrawHG)

495:   Input Parameters:
496: + hist - The histogram context
497: - color - one of the colors defined in petscdraw.h or PETSC_DRAW_ROTATE to make each bar a
498:           different color

500:   Level: intermediate

502: @*/
503: PetscErrorCode  PetscDrawHGSetColor(PetscDrawHG hist, int color)
504: {
507:   hist->color = color;
508:   return(0);
509: }

513: /*@
514:   PetscDrawHGSetLimits - Sets the axis limits for a histogram. If more
515:   points are added after this call, the limits will be adjusted to
516:   include those additional points.

518:   Not Collective (ignored except on processor 0 of PetscDrawHG)

520:   Input Parameters:
521: + hist - The histogram context
522: - x_min,x_max,y_min,y_max - The limits

524:   Level: intermediate

526:   Concepts: histogram^setting axis
527: @*/
528: PetscErrorCode  PetscDrawHGSetLimits(PetscDrawHG hist, PetscReal x_min, PetscReal x_max, int y_min, int y_max)
529: {
532:   hist->xmin = x_min;
533:   hist->xmax = x_max;
534:   hist->ymin = y_min;
535:   hist->ymax = y_max;
536:   return(0);
537: }

541: /*@
542:   PetscDrawHGCalcStats - Turns on calculation of descriptive statistics

544:   Not collective

546:   Input Parameters:
547: + hist - The histogram context
548: - calc - Flag for calculation

550:   Level: intermediate

552: .keywords:  draw, histogram, statistics

554: @*/
555: PetscErrorCode  PetscDrawHGCalcStats(PetscDrawHG hist, PetscBool calc)
556: {
559:   hist->calcStats = calc;
560:   return(0);
561: }

565: /*@
566:   PetscDrawHGIntegerBins - Turns on integer width bins

568:   Not collective

570:   Input Parameters:
571: + hist - The histogram context
572: - ints - Flag for integer width bins

574:   Level: intermediate

576: .keywords:  draw, histogram, statistics
577: @*/
578: PetscErrorCode  PetscDrawHGIntegerBins(PetscDrawHG hist, PetscBool ints)
579: {
582:   hist->integerBins = ints;
583:   return(0);
584: }

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

594:   Not Collective (ignored except on processor 0 of PetscDrawHG)

596:   Input Parameter:
597: . hist - The histogram context

599:   Output Parameter:
600: . axis - The axis context

602:   Level: intermediate

604: @*/
605: PetscErrorCode  PetscDrawHGGetAxis(PetscDrawHG hist, PetscDrawAxis *axis)
606: {
610:   *axis = hist->axis;
611:   return(0);
612: }

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: . win  - The draw context

627:   Level: intermediate

629: @*/
630: PetscErrorCode  PetscDrawHGGetDraw(PetscDrawHG hist, PetscDraw *win)
631: {
635:   *win = hist->win;
636:   return(0);
637: }