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 struct histogram *
  61 histogram_create (double bin_width, double min, double max)
  62 {
  63   int bins;
  64   struct histogram *h = xmalloc (sizeof *h);
  65   struct statistic *stat = &h->parent;
  66
  67   const double half_bin_width = bin_width / 2.0;
  68
  69   /* The lower and upper limits of the histogram, in units of half
  70      bin widths */
  71   int lower_limit, upper_limit;
  72
  73   /* -1 if the lower end of the range contains more unused space
  74      than the upper end.
  75      +1 otherwise.  */
  76   short sparse_end = 0;
  77
  78   double ul, ll;
  79   double lower_remainder = fabs (modf (min / half_bin_width, &ll));
  80   double upper_remainder = fabs (modf (max / half_bin_width, &ul));
  81
  82
  83   if (max == min)
  84     {
  85       msg (MW, _("Not creating histogram because the data contains less than 2 distinct values"));
  86       free (h);
  87       return NULL;
  88     }
  89
  90   assert (max > min);
  91
  92   lower_limit = ll;
  93
  94   /* If the minimum value is not aligned on a half bin width,
  95      then the lower bound must be extended so that the histogram range includes it. */
  96   if (lower_remainder > 0)
  97     lower_limit--;
  98
  99   /* However, the upper bound must be extended regardless, because histogram bins
 100      span the range [lower, upper) */
 101   upper_limit = ul + 1;
 102
 103   /* So, in the case of the maximum value coinciding with a half bin width,
 104      the upper end will be the sparse end (because is got extended by a complete
 105      half bin width).   In other cases, it depends which got the bigger extension. */
 106   if (upper_remainder == 0)
 107     sparse_end = +1;
 108   else
 109     sparse_end = lower_remainder < upper_remainder ? -1 : +1;
 110
 111   /* The range must be an EVEN number of half bin_widths */
 112   if ( (upper_limit - lower_limit) % 2)
 113     {
 114       /* Extend the range at the end which gives the least unused space */
 115       if (sparse_end == +1)
 116         lower_limit--;
 117       else
 118         upper_limit++;
 119
 120       /* Now the other end has more space */
 121       sparse_end *= -1;
 122     }
 123
 124   /* But the range should be aligned to an ODD number of
 125      half bin widths, so that the labels are aesthetically pleasing ones.
 126      Otherwise we are likely to get labels such as -3 -1 1 3 instead of -2 0 2 4
 127   */
 128   if ( lower_limit % 2 == 0)
 129     {
 130       /* Shift the range away from the sparse end, EXCEPT if that is the upper end,
 131          and it was extended to prevent the maximum value from getting lost */
 132       if (sparse_end == +1 && upper_remainder > 0)
 133         {
 134           lower_limit --;
 135           upper_limit --;
 136         }
 137       else
 138         {
 139           lower_limit ++;
 140           upper_limit ++;
 141         }
 142     }
 143   bins = (upper_limit - lower_limit) / 2.0;
 144
 145   /* Force the number of bins to lie in a sensible range.
 146      FIXME: We should redo the above calculations if this happens! */
 147   if (bins > 25)
 148     bins = 25;
 149
 150   if (bins < 1)
 151     bins = 1;
 152
 153   h->gsl_hist = gsl_histogram_alloc (bins);
 154
 155   gsl_histogram_set_ranges_uniform (h->gsl_hist,
 156                                     lower_limit * half_bin_width,
 157                                     upper_limit * half_bin_width);
 158
 159   stat->accumulate = acc;
 160   stat->destroy = destroy;
 161
 162   return h;
 163 }
 164