X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Fhistogram.c;h=b00067210f08676736d589660c9eb85f75216831;hb=e2da62d735c597afeef2e0e9b36e5a4a83d7da94;hp=d14eb45e3423a55077c84d8b93a370d70b6c383b;hpb=478954c1ed028744ab870f0e6a5befcc8f643758;p=pspp

diff --git a/src/math/histogram.c b/src/math/histogram.c
index d14eb45e34..b00067210f 100644
--- a/src/math/histogram.c
+++ b/src/math/histogram.c
@@ -21,6 +21,7 @@
 #include <gsl/gsl_histogram.h>
 #include <math.h>
 
+#include "data/settings.h"
 #include "libpspp/message.h"
 #include "libpspp/assertion.h"
 #include "libpspp/cast.h"
@@ -57,79 +58,197 @@ destroy (struct statistic *s)
 }
 
 
-struct histogram *
-histogram_create (double bin_width, double min, double max)
+static
+double get_slack (double limit, double half_bin_width, int *n_half_bins)
 {
-  int bins;
-  struct histogram *h = xmalloc (sizeof *h);
-  struct statistic *stat = &h->parent;
+  double ipart, remainder;
+
+  assert (half_bin_width > 0);
+
+  remainder =  modf (limit / half_bin_width, &ipart);
+
+  /* In C modf and % behave in an unexpected (to me at any rate) manner
+     when presented with a negative value
+
+     For example, modf (-7.0 / 3.0) returns -2.0 R -0.3333
+  */
+
+  
+  *n_half_bins = ipart;
+
+  return remainder * half_bin_width;
+}
+
 
+/* This functions adjusts the upper and lower range of the histogram to make them fit BIN_WIDTH
+   MIN and MAX are the lowest and highest data to be plotted in the histogram.
+   ADJ_MIN and ADJ_MAX are locations of the adjusted values of MIN and MAX (the range will
+   always be  equal or slightly larger).
+   Returns the number of bins.
+
+   The "testing_assert" expressions in this function should be algebraically correct.
+   However, due to floating point rounding they could fail, especially when small numbers
+   are involved.  In normal use, therefore, testing_assert does nothing.
+ */
+static int
+adjust_bin_ranges (double bin_width, double min, double max, double *adj_min, double *adj_max)
+{
   const double half_bin_width = bin_width / 2.0;
 
   /* The lower and upper limits of the histogram, in units of half
      bin widths */
   int lower_limit, upper_limit;
 
-  /* -1 if the lower end of the range contains more unused space
-     than the upper end.
-     +1 otherwise.  */
-  short sparse_end = 0;
+  double lower_slack =  get_slack (min, half_bin_width, &lower_limit);
+  double upper_slack = -get_slack (max, half_bin_width, &upper_limit);
 
-  if (max == min)
+  testing_assert (max > min);
+
+  /* If min is negative, then lower_slack may be less than zero.
+     In this case, the lower bound must be extended in the negative direction
+     so that it is less than OR EQUAL to min.
+   */
+  if (lower_slack < 0)
     {
-      msg (MW, _("Not creating histogram because the data contains less than 2 distinct values"));
-      free (h);
-      return NULL;
+      lower_limit--;
+      lower_slack += half_bin_width;
     }
+  testing_assert (lower_limit * half_bin_width <= min);
 
-  assert (max > min);
-
-  {
-    double ul, ll;
-    double lower_tail = modf (min / half_bin_width, &ll);
-    double upper_tail = modf (max / half_bin_width, &ul);
-    lower_limit = ll - 1;
-    upper_limit = ul + 1;
-    
-    sparse_end = lower_tail < upper_tail ? -1 : +1;
-  }
+  /* However, the upper bound must be extended regardless, because histogram bins
+     span the range [lower, upper). In other words, the upper bound must be
+     greater than max.
+  */
+  upper_limit++;;
+  upper_slack += half_bin_width;
+  testing_assert (upper_limit * half_bin_width > max);
 
   /* The range must be an EVEN number of half bin_widths */
   if ( (upper_limit - lower_limit) % 2)
     {
       /* Extend the range at the end which gives the least unused space */
-      if (sparse_end == +1)
-	lower_limit--;
+      if (upper_slack > lower_slack)
+        {
+          lower_limit--;
+          lower_slack += half_bin_width;
+        }
       else
-        upper_limit++;
-      
-      /* Now the other end has more space */
-      sparse_end *= -1;
+        {
+          upper_limit++;
+          upper_slack += half_bin_width;
+        }
     }
 
   /* But the range should be aligned to an ODD number of
-     half bin widths, so that the labels are aesthetically pleasing ones. */
+     half bin widths, so that the labels are aesthetically pleasing ones.
+     Otherwise we are likely to get labels such as -3 -1 1 3 instead of -2 0 2 4
+  */
   if ( lower_limit % 2 == 0)
     {
-      lower_limit += -sparse_end ;
-      upper_limit += -sparse_end ;
+      /* If there is not enough slack at either end to perform a shift,
+         then we must extend the range so that there is.  We must extend
+         by two half bin widths in order to preserve the EVEN condition
+         established above.  Also, we extend on the end with the least
+         slack, in order to keep things as balanced as possible. */
+      if ( upper_slack > lower_slack && upper_slack <= half_bin_width)
+        {
+          lower_limit -= 2;
+          lower_slack += 2 * half_bin_width;
+        }
+           
+      if (lower_slack > upper_slack && lower_slack < half_bin_width)
+        {
+          upper_limit += 2;
+          upper_slack += 2 * half_bin_width;
+        }
+
+      if (upper_slack > lower_slack)
+        {
+          testing_assert (upper_slack > half_bin_width);
+
+          /* Adjust the range to the left */
+          lower_limit --;
+          upper_limit --;
+          upper_slack -= half_bin_width;
+          lower_slack += half_bin_width;
+        }
+      else
+        {
+          testing_assert (lower_slack >= half_bin_width);
+
+          /* Adjust the range to the right */
+          lower_limit ++;
+          upper_limit ++;
+          lower_slack -= half_bin_width;
+          upper_slack += half_bin_width;
+        }
     }
 
-  bins = (upper_limit - lower_limit) / 2.0;
+  /* If there are any completely empty bins, then remove them,
+     since empty bins don't really add much information to the histogram.
+   */
+  if (upper_slack > 2 * half_bin_width)
+    {
+      upper_slack -= 2 * half_bin_width;
+      upper_limit -=2;
+    }
 
-  /* Force the number of bins to lie in a sensible range */
-  if (bins > 25) 
-    bins = 25;
+  if (lower_slack >= 2 * half_bin_width)
+    {
+      lower_slack -= 2 * half_bin_width;
+      lower_limit +=2;
+    }
 
+  *adj_min = lower_limit * half_bin_width;
+  *adj_max = upper_limit * half_bin_width;
+
+  testing_assert (*adj_max > max);
+  testing_assert (*adj_min <= min);
+
+  return (upper_limit - lower_limit) / 2.0;
+}
+
+
+
+struct histogram *
+histogram_create (double bin_width, double min, double max)
+{
+  const int MAX_BINS = 25;
+  struct histogram *h;
+  struct statistic *stat;
+  int bins;
+  double adjusted_min, adjusted_max;
+
+  if (max == min)
+    {
+      msg (MW, _("Not creating histogram because the data contains less than 2 distinct values"));
+      return NULL;
+    }
+
+  assert (bin_width > 0);
+
+  bin_width = chart_rounded_tick (bin_width);
+  bins = adjust_bin_ranges (bin_width, min, max, &adjusted_min, &adjusted_max);
+
+  /* Force the number of bins to lie in a sensible range. */
+  if (bins > MAX_BINS) 
+    {
+      bin_width = chart_rounded_tick ((max - min) / (double) (MAX_BINS - 1));
+      bins = adjust_bin_ranges (bin_width,
+                                min, max, &adjusted_min, &adjusted_max);
+    }
+
+  /* Can this ever happen? */
   if (bins < 1)
     bins = 1;
 
+  h = xmalloc (sizeof *h);
+
   h->gsl_hist = gsl_histogram_alloc (bins);
 
-  gsl_histogram_set_ranges_uniform (h->gsl_hist,
-                                    lower_limit * half_bin_width,
-                                    upper_limit * half_bin_width);
+  gsl_histogram_set_ranges_uniform (h->gsl_hist, adjusted_min, adjusted_max);
 
+  stat = &h->parent;
   stat->accumulate = acc;
   stat->destroy = destroy;