pintos-os.org Git - pspp/blob - src/math/histogram.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2004, 2008, 2009, 2011, 2012 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18
  19 #include "math/histogram.h"
  20
  21 #include <gsl/gsl_histogram.h>
  22 #include <math.h>
  23
  24 #include "libpspp/message.h"
  25 #include "libpspp/assertion.h"
  26 #include "libpspp/cast.h"
  27 #include "math/chart-geometry.h"
  28
  29 #include "gettext.h"
  30 #define _(msgid) gettext (msgid)
  31 #define N_(msgid) msgid
  32
  33
  34 #include "gl/xalloc.h"
  35
  36 void
  37 histogram_add (struct histogram *h, double y, double c)
  38 {
  39   struct statistic *stat = &h->parent;
  40   stat->accumulate (stat, NULL, c, 0, y);
  41 }
  42
  43 static void
  44 acc (struct statistic *s, const struct ccase *cx UNUSED, double c, double cc UNUSED, double y)
  45 {
  46   struct histogram *hist = UP_CAST (s, struct histogram, parent);
  47
  48   gsl_histogram_accumulate (hist->gsl_hist, y, c);
  49 }
  50
  51 static void
  52 destroy (struct statistic *s)
  53 {
  54   struct histogram *h = UP_CAST (s, struct histogram, parent);
  55   gsl_histogram_free (h->gsl_hist);
  56   free (s);
  57 }
  58
  59
  60 static
  61 double get_slack (double limit, double half_bin_width, int *n_half_bins)
  62 {
  63   double ipart, remainder;
  64
  65   assert (half_bin_width > 0);
  66
  67   remainder =  modf (limit / half_bin_width, &ipart);
  68
  69   /* In C modf and % behave in an unexpected (to me at any rate) manner
  70      when presented with a negative value
  71
  72      For example, modf (-7.0 / 3.0) returns -2.0 R -0.3333
  73   */
  74
  75
  76   *n_half_bins = ipart;
  77
  78   return remainder * half_bin_width;
  79 }
  80
  81
  82 /* This functions adjusts the upper and lower range of the histogram to make them fit BIN_WIDTH
  83    MIN and MAX are the lowest and highest data to be plotted in the histogram.
  84    ADJ_MIN and ADJ_MAX are locations of the adjusted values of MIN and MAX (the range will
  85    always be  equal or slightly larger).
  86    Returns the number of bins.
  87  */
  88 static int
  89 adjust_bin_ranges (double bin_width, double min, double max, double *adj_min, double *adj_max)
  90 {
  91   const double half_bin_width = bin_width / 2.0;
  92
  93   /* The lower and upper limits of the histogram, in units of half
  94      bin widths */
  95   int lower_limit, upper_limit;
  96
  97   double lower_slack =  get_slack (min, half_bin_width, &lower_limit);
  98   double upper_slack = -get_slack (max, half_bin_width, &upper_limit);
  99
 100   assert (max > min);
 101
 102   /* If min is negative, then lower_slack may be less than zero.
 103      In this case, the lower bound must be extended in the negative direction
 104      so that it is less than OR EQUAL to min.
 105    */
 106   if (lower_slack < 0)
 107     {
 108       lower_limit--;
 109       lower_slack += half_bin_width;
 110     }
 111   assert (lower_limit * half_bin_width <= min);
 112
 113   /* However, the upper bound must be extended regardless, because histogram bins
 114      span the range [lower, upper). In other words, the upper bound must be
 115      greater than max.
 116   */
 117   upper_limit++;;
 118   upper_slack += half_bin_width;
 119   assert (upper_limit * half_bin_width > max);
 120
 121   /* The range must be an EVEN number of half bin_widths */
 122   if ( (upper_limit - lower_limit) % 2)
 123     {
 124       /* Extend the range at the end which gives the least unused space */
 125       if (upper_slack > lower_slack)
 126         {
 127           lower_limit--;
 128           lower_slack += half_bin_width;
 129         }
 130       else
 131         {
 132           upper_limit++;
 133           upper_slack += half_bin_width;
 134         }
 135     }
 136
 137   /* But the range should be aligned to an ODD number of
 138      half bin widths, so that the labels are aesthetically pleasing ones.
 139      Otherwise we are likely to get labels such as -3 -1 1 3 instead of -2 0 2 4
 140   */
 141   if ( lower_limit % 2 == 0)
 142     {
 143       /* If there is not enough slack at either end to perform a shift,
 144          then we must extend the range so that there is.  We must extend
 145          by two half bin widths in order to preserve the EVEN condition
 146          established above.  Also, we extend on the end with the least
 147          slack, in order to keep things as balanced as possible. */
 148       if ( upper_slack > lower_slack && upper_slack <= half_bin_width)
 149         {
 150           lower_limit -= 2;
 151           lower_slack += 2 * half_bin_width;
 152         }
 153
 154       if (lower_slack > upper_slack && lower_slack < half_bin_width)
 155         {
 156           upper_limit += 2;
 157           upper_slack += 2 * half_bin_width;
 158         }
 159
 160       if (upper_slack > lower_slack)
 161         {
 162           assert (upper_slack > half_bin_width);
 163
 164           /* Adjust the range to the left */
 165           lower_limit --;
 166           upper_limit --;
 167           upper_slack -= half_bin_width;
 168           lower_slack += half_bin_width;
 169         }
 170       else
 171         {
 172           assert (lower_slack >= half_bin_width);
 173
 174           /* Adjust the range to the right */
 175           lower_limit ++;
 176           upper_limit ++;
 177           lower_slack -= half_bin_width;
 178           upper_slack += half_bin_width;
 179         }
 180     }
 181
 182   /* If there are any completely empty bins, then remove them,
 183      since empty bins don't really add much information to the histogram.
 184    */
 185   if (upper_slack > 2 * half_bin_width)
 186     {
 187       upper_slack -= 2 * half_bin_width;
 188       upper_limit -=2;
 189     }
 190
 191   if (lower_slack >= 2 * half_bin_width)
 192     {
 193       lower_slack -= 2 * half_bin_width;
 194       lower_limit +=2;
 195     }
 196
 197   *adj_min = lower_limit * half_bin_width;
 198   *adj_max = upper_limit * half_bin_width;
 199
 200   assert (*adj_max > max);
 201   assert (*adj_min <= min);
 202
 203   return (upper_limit - lower_limit) / 2.0;
 204 }
 205
 206
 207
 208 struct histogram *
 209 histogram_create (double bin_width, double min, double max)
 210 {
 211   const int MAX_BINS = 25;
 212   struct histogram *h;
 213   struct statistic *stat;
 214   int bins;
 215   double adjusted_min, adjusted_max;
 216
 217   if (max == min)
 218     {
 219       msg (MW, _("Not creating histogram because the data contains less than 2 distinct values"));
 220       return NULL;
 221     }
 222
 223   assert (bin_width > 0);
 224
 225   bins = adjust_bin_ranges (bin_width, min, max, &adjusted_min, &adjusted_max);
 226
 227   /* Force the number of bins to lie in a sensible range. */
 228   if (bins > MAX_BINS)
 229     {
 230       bins = adjust_bin_ranges ((max - min) / (double) (MAX_BINS - 1),
 231                                 min, max, &adjusted_min, &adjusted_max);
 232     }
 233
 234   /* Can this ever happen? */
 235   if (bins < 1)
 236     bins = 1;
 237
 238   h = xmalloc (sizeof *h);
 239
 240   h->gsl_hist = gsl_histogram_alloc (bins);
 241
 242   gsl_histogram_set_ranges_uniform (h->gsl_hist, adjusted_min, adjusted_max);
 243
 244   stat = &h->parent;
 245   stat->accumulate = acc;
 246   stat->destroy = destroy;
 247
 248   return h;
 249 }
 250