gsl_histogram_flat.c

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_errno.h>

#define MAXG 2048

static gsl_histogram *
ini_histogram(const int,
              const double,
              const double);


static short
histogram_is_flat(const gsl_histogram *);


int bins = 10;

int
main(void)
{
  gsl_histogram *m      = NULL;
  short         is_flat = 0;

  m = ini_histogram(bins, 0, 100);
  gsl_histogram_increment(m, 5);
  gsl_histogram_increment(m, 25);
  gsl_histogram_increment(m, 55);
  gsl_histogram_increment(m, 57);
  gsl_histogram_increment(m, 65);
  gsl_histogram_fprintf(stderr, m, "%6.3g", "%6g");

  if (histogram_is_flat(m))
    printf("Histogram is flat\n");
  else
    printf("Histogram is NOT flat\n");

  return 0;
}


/* ==== */
static short
histogram_is_flat(const gsl_histogram *z)
{
  double  val, avg, sum = 0.0, flat = 0.8;
  size_t  lbin, gbin;      /* lowest/greatest populated bin */
  int     i, b = 0, is_flat = 1;

  fprintf(stderr, "[[histogram_is_flat()]]\n");
  for (i = 0; i < bins; i++) {
    val = gsl_histogram_get(z, i);
    if (val > 0) {
      lbin = i;
      break;
    }
  }

  printf("lbin is %i\n", lbin);
  for (i = bins - 1; i >= 0; i--) {
    val = gsl_histogram_get(z, i);
    if (val > 0) {
      gbin = i;
      break;
    }
  }
  printf("gbin is %i\n", gbin);

  // compute average over histogram
  for (i = lbin; i <= gbin; i++) {
    if ((val = gsl_histogram_get(z, i)) != 0) {
      sum += val;
      b++;
    }
  }
  avg = sum / b;

  // evaluate if populated bins are all within accepted range
  for (i = lbin; i <= gbin; i++) {
    if ((val = gsl_histogram_get(z, i)) != 0) {
      if (val < flat * avg) {
        is_flat = 0;
        break;
      }
    }
  }


  return is_flat;
}


/* ==== */
static gsl_histogram *
ini_histogram(const int     n,
              const double  min,
              const double  max)
{
  gsl_histogram *a = gsl_histogram_alloc(n);

  gsl_histogram_set_ranges_uniform(a, min, max);
  return a;
}