Improve the way in which histogram bin ranges are chosen
[pspp] / src / math / histogram.c
index 0c641b56e127ddbdf9a180363629101fad746b6c..0b2369f63c7b6ecad3329d6c728e2ed5eebf7dcf 100644 (file)
@@ -1,5 +1,5 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 2004, 2008, 2009 Free Software Foundation, Inc.
+   Copyright (C) 2004, 2008, 2009, 2011 Free Software Foundation, Inc.
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
 
 #include <config.h>
-#include "histogram.h"
 
-#include <gl/xalloc.h>
-#include <libpspp/assertion.h>
-#include <libpspp/cast.h>
+#include "math/histogram.h"
 
 #include <gsl/gsl_histogram.h>
-#include "chart-geometry.h"
 #include <math.h>
 
+#include "libpspp/assertion.h"
+#include "libpspp/cast.h"
+#include "math/chart-geometry.h"
+
+#include "gl/xalloc.h"
 
 void
 histogram_add (struct histogram *h, double y, double c)
@@ -54,27 +55,39 @@ destroy (struct statistic *s)
 
 
 struct histogram *
-histogram_create (int bins, double min, double max)
+histogram_create (double bin_width, double min, double max)
 {
+  int bins;
   struct histogram *h = xmalloc (sizeof *h);
   struct statistic *stat = &h->parent;
+
   double upper_limit, lower_limit;
 
-  double bin_width = chart_rounded_tick ((max - min) / (double) bins);
-  double bin_width_2 = bin_width / 2.0;
+  assert (max >= min);
 
-  int n =  ceil (max / (bin_width_2) ) ;
+  if (max == min)
+     bin_width = 1;
 
-  assert (max >= min);
+  lower_limit = floor (2 * min / bin_width) - 1;
+  upper_limit = floor (2 * max / bin_width) + 1;
+  
+  /* The range must be an even number of half bin_widths */
+  if ( (int)(upper_limit - lower_limit) % 2)
+    {
+      /* Extend the range at the end which gives the least unused space */
+      if (remainder (min, bin_width / 2.0) > remainder (max, bin_width / 2.0))
+       lower_limit --;
+      else
+       upper_limit ++;
+    }
 
-  if ( ! (n % 2 ) ) n++;
-  upper_limit = n * bin_width_2;
+  bins = (upper_limit - lower_limit) / 2.0;
 
-  n =  floor (min / (bin_width_2) ) ;
-  if ( ! (n % 2 ) ) n--;
-  lower_limit = n * bin_width_2;
+  upper_limit *= bin_width / 2;
+  lower_limit *= bin_width / 2;  
 
   h->gsl_hist = gsl_histogram_alloc (bins);
+
   gsl_histogram_set_ranges_uniform (h->gsl_hist, lower_limit, upper_limit);
 
   stat->accumulate = acc;