Actual source code: hists.c

      1: #define PETSC_DLL
      2: /*
      3:   Contains the data structure for plotting a histogram in a window with an axis.
      4: */

     6:  #include petsc.h

      8: PetscCookie DRAWHG_COOKIE = 0;

     10: struct _p_DrawHG {
     11:   PETSCHEADER(int);
     12:   PetscErrorCode (*destroy)(PetscDrawSP);
     13:   PetscErrorCode (*view)(PetscDrawSP,PetscViewer);
     14:   PetscDraw     win;
     15:   PetscDrawAxis axis;
     16:   PetscReal     xmin,xmax;
     17:   PetscReal     ymin,ymax;
     18:   int           numBins;
     19:   int           maxBins;
     20:   PetscReal     *bins;
     21:   int           numValues;
     22:   int           maxValues;
     23:   PetscReal     *values;
     24:   int           color;
     25:   PetscTruth    calcStats;
     26:   PetscTruth    integerBins;
     27: };

     29: #define CHUNKSIZE 100

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

     36:    Collective over PetscDraw

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

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

     45:    Level: intermediate

     47:    Contributed by: Matthew Knepley

     49:    Concepts: histogram^creating

     51: .seealso: PetscDrawHGDestroy()

     53: @*/
     54: PetscErrorCode PETSC_DLLEXPORT PetscDrawHGCreate(PetscDraw draw, int bins, PetscDrawHG *hist) {
     55:   PetscDrawHG h;
     56:   MPI_Comm    comm;
     57:   PetscTruth  isnull;

     63:   PetscObjectGetComm((PetscObject) draw, &comm);
     64:   PetscHeaderCreate(h, _p_DrawHG, int, DRAWHG_COOKIE, 0, "PetscDrawHG", comm, PetscDrawHGDestroy, PETSC_NULL);
     65:   h->view        = PETSC_NULL;
     66:   h->destroy     = PETSC_NULL;
     67:   h->win         = draw;
     68:   PetscObjectReference((PetscObject) draw);
     69:   h->color       = PETSC_DRAW_GREEN;
     70:   h->xmin        = PETSC_MAX;
     71:   h->xmax        = PETSC_MIN;
     72:   h->ymin        = 0.;
     73:   h->ymax        = 1.;
     74:   h->numBins     = bins;
     75:   h->maxBins     = bins;
     76:   PetscMalloc(h->maxBins * sizeof(PetscReal), &h->bins);
     77:   h->numValues   = 0;
     78:   h->maxValues   = CHUNKSIZE;
     79:   h->calcStats   = PETSC_FALSE;
     80:   h->integerBins = PETSC_FALSE;
     81:   PetscMalloc(h->maxValues * sizeof(PetscReal), &h->values);
     82:   PetscLogObjectMemory(h, (h->maxBins + h->maxValues)*sizeof(PetscReal));
     83:   PetscTypeCompare((PetscObject) draw, PETSC_DRAW_NULL, &isnull);
     84:   if (!isnull) {
     85:     PetscDrawAxisCreate(draw, &h->axis);
     86:     PetscLogObjectParent(h, h->axis);
     87:   } else {
     88:     h->axis = PETSC_NULL;
     89:   }
     90:   *hist = h;
     91:   return(0);
     92: }

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

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

    101:    Input Parameter:
    102: +  hist - The histogram context.
    103: -  dim  - The number of curves.

    105:    Level: intermediate

    107:   Contributed by: Matthew Knepley

    109:    Concepts: histogram^setting number of bins

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

    117:   if (hist->maxBins < bins) {
    118:     PetscFree(hist->bins);
    119:     PetscMalloc(bins * sizeof(PetscReal), &hist->bins);
    120:     PetscLogObjectMemory(hist, (bins - hist->maxBins) * sizeof(PetscReal));
    121:     hist->maxBins = bins;
    122:   }
    123:   hist->numBins = bins;
    124:   return(0);
    125: }

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

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

    134:   Input Parameter:
    135: . hist - The histogram context.

    137:   Level: intermediate

    139:   Contributed by: Matthew Knepley

    141:   Concepts: histogram^resetting
    142: @*/
    143: PetscErrorCode PETSC_DLLEXPORT PetscDrawHGReset(PetscDrawHG hist)
    144: {
    147:   hist->xmin      = PETSC_MAX;
    148:   hist->xmax      = PETSC_MIN;
    149:   hist->ymin      = 0;
    150:   hist->ymax      = 0;
    151:   hist->numValues = 0;
    152:   return(0);
    153: }

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

    160:   Collective over PetscDrawHG

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

    165:   Level: intermediate

    167:   Contributed by: Matthew Knepley

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

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

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

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

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

    200:   Level: intermediate

    202:   Contributed by: Matthew Knepley

    204:   Concepts: histogram^adding values

    206: .seealso: PetscDrawHGAddValues()
    207: @*/
    208: PetscErrorCode PETSC_DLLEXPORT PetscDrawHGAddValue(PetscDrawHG hist, PetscReal value)
    209: {
    212:   /* Allocate more memory if necessary */
    213:   if (hist->numValues >= hist->maxValues) {
    214:     PetscReal *tmp;

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

    257:   hist->values[hist->numValues++] = value;
    258:   return(0);
    259: }

    263: /*@
    264:   PetscDrawHGDraw - Redraws a histogram.

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

    268:   Input Parameter:
    269: . hist - The histogram context

    271:   Level: intermediate

    273:   Contributed by: Matthew Knepley
    274: @*/
    275: PetscErrorCode PETSC_DLLEXPORT PetscDrawHGDraw(PetscDrawHG hist)
    276: {
    277:   PetscDraw  draw = hist->win;
    278:   PetscTruth isnull;
    279:   PetscReal  xmin,xmax,ymin,ymax,*bins,*values,binSize,binLeft,binRight,maxHeight,mean,var;
    280:   char       title[256];
    281:   char       xlabel[256];
    282:   int        numBins,numBinsOld,numValues,initSize,i,p,bcolor,color;

    287:   PetscTypeCompare((PetscObject) draw, PETSC_DRAW_NULL, &isnull);
    288:   if (isnull) return(0);
    289:   if ((hist->xmin >= hist->xmax) || (hist->ymin >= hist->ymax)) return(0);
    290:   if (hist->numValues < 1) return(0);

    292: #if 0
    293:   MPI_Comm_rank(hist->comm,&rank);
    294:   if (rank) return(0);
    295: #endif

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

    356:     PetscMemzero(bins, numBins * sizeof(PetscReal));
    357:     maxHeight = 0;
    358:     for (i = 0; i < numBins; i++) {
    359:       binLeft   = xmin + binSize*i;
    360:       binRight  = xmin + binSize*(i+1);
    361:       for(p = 0; p < numValues; p++) {
    362:         if ((values[p] >= binLeft) && (values[p] < binRight)) bins[i]++;
    363:         /* Handle last bin separately */
    364:         if ((i == numBins-1) && (values[p] == binRight)) bins[i]++;
    365:         if (!i) {
    366:           mean += values[p];
    367:           var  += values[p]*values[p];
    368:         }
    369:       }
    370:       maxHeight = PetscMax(maxHeight, bins[i]);
    371:     }
    372:     if (maxHeight > ymax) ymax = hist->ymax = maxHeight;

    374:     PetscDrawAxisSetLimits(hist->axis, xmin, xmax, ymin, ymax);
    375:     if (hist->calcStats) {
    376:       mean /= numValues;
    377:       if (numValues > 1) {
    378:         var = (var - numValues*mean*mean) / (numValues-1);
    379:       } else {
    380:         var = 0.0;
    381:       }
    382:       sprintf(title, "Mean: %g  Var: %g", mean, var);
    383:       sprintf(xlabel, "Total: %d", numValues);
    384:       PetscDrawAxisSetLabels(hist->axis, title, xlabel, PETSC_NULL);
    385:     }
    386:     PetscDrawAxisDraw(hist->axis);
    387:     /* Draw bins */
    388:     for (i = 0; i < numBins; i++) {
    389:       binLeft   = xmin + binSize*i;
    390:       binRight  = xmin + binSize*(i+1);
    391:       PetscDrawRectangle(draw,binLeft,ymin,binRight,bins[i],bcolor,bcolor,bcolor,bcolor);
    392:       if (color == PETSC_DRAW_ROTATE && bins[i]) bcolor++; if (bcolor > 31) bcolor = 2;
    393:       PetscDrawLine(draw,binLeft,ymin,binLeft,bins[i],PETSC_DRAW_BLACK);
    394:       PetscDrawLine(draw,binRight,ymin,binRight,bins[i],PETSC_DRAW_BLACK);
    395:       PetscDrawLine(draw,binLeft,bins[i],binRight,bins[i],PETSC_DRAW_BLACK);
    396:     }
    397:     PetscDrawHGSetNumberBins(hist, numBinsOld);
    398:   }
    399:   PetscDrawSynchronizedFlush(draw);
    400:   PetscDrawPause(draw);
    401:   return(0);
    402: }

    406: /*@
    407:   PetscDrawHGPrint - Prints the histogram information.

    409:   Not collective

    411:   Input Parameter:
    412: . hist - The histogram context

    414:   Level: beginner

    416:   Contributed by: Matthew Knepley

    418: .keywords:  draw, histogram
    419: @*/
    420: PetscErrorCode PETSC_DLLEXPORT PetscDrawHGPrint(PetscDrawHG hist)
    421: {
    422:   PetscReal xmax,xmin,*bins,*values,binSize,binLeft,binRight,mean,var;
    424:   int       numBins,numBinsOld,numValues,initSize,i,p;

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

    431:   xmax      = hist->xmax;
    432:   xmin      = hist->xmin;
    433:   numValues = hist->numValues;
    434:   values    = hist->values;
    435:   mean      = 0.0;
    436:   var       = 0.0;
    437:   if (xmax == xmin) {
    438:     /* Calculate number of points in the bin */
    439:     bins    = hist->bins;
    440:     bins[0] = 0;
    441:     for(p = 0; p < numValues; p++) {
    442:       if (values[p] == xmin) bins[0]++;
    443:       mean += values[p];
    444:       var  += values[p]*values[p];
    445:     }
    446:     /* Draw bins */
    447:     PetscPrintf(hist->comm, "Bin %2d (%6.2g - %6.2g): %.0g\n", 0, xmin, xmax, bins[0]);
    448:   } else {
    449:     numBins    = hist->numBins;
    450:     numBinsOld = hist->numBins;
    451:     if (hist->integerBins && (((int) xmax - xmin) + 1.0e-05 > xmax - xmin)) {
    452:       initSize = (int) ((int) xmax - xmin)/numBins;
    453:       while (initSize*numBins != (int) xmax - xmin) {
    454:         initSize = PetscMax(initSize - 1, 1);
    455:         numBins  = (int) ((int) xmax - xmin)/initSize;
    456:         PetscDrawHGSetNumberBins(hist, numBins);
    457:       }
    458:     }
    459:     binSize = (xmax - xmin)/numBins;
    460:     bins    = hist->bins;

    462:     /* Calculate number of points in each bin */
    463:     PetscMemzero(bins, numBins * sizeof(PetscReal));
    464:     for (i = 0; i < numBins; i++) {
    465:       binLeft   = xmin + binSize*i;
    466:       binRight  = xmin + binSize*(i+1);
    467:       for(p = 0; p < numValues; p++) {
    468:         if ((values[p] >= binLeft) && (values[p] < binRight)) bins[i]++;
    469:         /* Handle last bin separately */
    470:         if ((i == numBins-1) && (values[p] == binRight)) bins[i]++;
    471:         if (!i) {
    472:           mean += values[p];
    473:           var  += values[p]*values[p];
    474:         }
    475:       }
    476:     }
    477:     /* Draw bins */
    478:     for (i = 0; i < numBins; i++) {
    479:       binLeft   = xmin + binSize*i;
    480:       binRight  = xmin + binSize*(i+1);
    481:       PetscPrintf(hist->comm, "Bin %2d (%6.2g - %6.2g): %.0g\n", i, binLeft, binRight, bins[i]);
    482:     }
    483:     PetscDrawHGSetNumberBins(hist, numBinsOld);
    484:   }

    486:   if (hist->calcStats) {
    487:     mean /= numValues;
    488:     if (numValues > 1) {
    489:       var = (var - numValues*mean*mean) / (numValues-1);
    490:     } else {
    491:       var = 0.0;
    492:     }
    493:     PetscPrintf(hist->comm, "Mean: %g  Var: %g\n", mean, var);
    494:     PetscPrintf(hist->comm, "Total: %d\n", numValues);
    495:   }
    496:   return(0);
    497: }
    498:
    501: /*@
    502:   PetscDrawHGSetColor - Sets the color the bars will be drawn with.

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

    506:   Input Parameters:
    507: + hist - The histogram context
    508: - color - one of the colors defined in petscdraw.h or PETSC_DRAW_ROTATE to make each bar a
    509:           different color

    511:   Level: intermediate

    513: @*/
    514: PetscErrorCode PETSC_DLLEXPORT PetscDrawHGSetColor(PetscDrawHG hist, int color)
    515: {
    518:   hist->color = color;
    519:   return(0);
    520: }

    524: /*@
    525:   PetscDrawHGSetLimits - Sets the axis limits for a histogram. If more
    526:   points are added after this call, the limits will be adjusted to
    527:   include those additional points.

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

    531:   Input Parameters:
    532: + hist - The histogram context
    533: - x_min,x_max,y_min,y_max - The limits

    535:   Level: intermediate

    537:   Contributed by: Matthew Knepley

    539:   Concepts: histogram^setting axis
    540: @*/
    541: PetscErrorCode PETSC_DLLEXPORT PetscDrawHGSetLimits(PetscDrawHG hist, PetscReal x_min, PetscReal x_max, int y_min, int y_max)
    542: {
    545:   hist->xmin = x_min;
    546:   hist->xmax = x_max;
    547:   hist->ymin = y_min;
    548:   hist->ymax = y_max;
    549:   return(0);
    550: }

    554: /*@
    555:   PetscDrawHGCalcStats - Turns on calculation of descriptive statistics

    557:   Not collective

    559:   Input Parameters:
    560: + hist - The histogram context
    561: - calc - Flag for calculation

    563:   Level: intermediate

    565:   Contributed by: Matthew Knepley

    567: .keywords:  draw, histogram, statistics

    569: @*/
    570: PetscErrorCode PETSC_DLLEXPORT PetscDrawHGCalcStats(PetscDrawHG hist, PetscTruth calc)
    571: {
    574:   hist->calcStats = calc;
    575:   return(0);
    576: }

    580: /*@
    581:   PetscDrawHGIntegerBins - Turns on integer width bins

    583:   Not collective

    585:   Input Parameters:
    586: + hist - The histogram context
    587: - ints - Flag for integer width bins

    589:   Level: intermediate

    591:   Contributed by: Matthew Knepley

    593: .keywords:  draw, histogram, statistics
    594: @*/
    595: PetscErrorCode PETSC_DLLEXPORT PetscDrawHGIntegerBins(PetscDrawHG hist, PetscTruth ints)
    596: {
    599:   hist->integerBins = ints;
    600:   return(0);
    601: }

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

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

    613:   Input Parameter:
    614: . hist - The histogram context

    616:   Output Parameter:
    617: . axis - The axis context

    619:   Level: intermediate

    621:   Contributed by: Matthew Knepley
    622: @*/
    623: PetscErrorCode PETSC_DLLEXPORT PetscDrawHGGetAxis(PetscDrawHG hist, PetscDrawAxis *axis)
    624: {
    628:   *axis = hist->axis;
    629:   return(0);
    630: }

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

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

    639:   Input Parameter:
    640: . hist - The histogram context

    642:   Output Parameter:
    643: . win  - The draw context

    645:   Level: intermediate

    647:   Contributed by: Matthew Knepley
    648: @*/
    649: PetscErrorCode PETSC_DLLEXPORT PetscDrawHGGetDraw(PetscDrawHG hist, PetscDraw *win)
    650: {
    654:   *win = hist->win;
    655:   return(0);
    656: }