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 "data/settings.h"
  25 #include "libpspp/message.h"
  26 #include "libpspp/assertion.h"
  27 #include "libpspp/cast.h"
  28 #include "math/chart-geometry.h"
  29
  30 #include "gettext.h"
  31 #define _(msgid) gettext (msgid)
  32 #define N_(msgid) msgid
  33
  34
  35 #include "gl/xalloc.h"
  36
  37 void
  38 histogram_add (struct histogram *h, double y, double c)
  39 {
  40   struct statistic *stat = &h->parent;
  41   stat->accumulate (stat, NULL, c, 0, y);
  42 }
  43
  44 static void
  45 acc (struct statistic *s, const struct ccase *cx UNUSED, double c, double cc UNUSED, double y)
  46 {
  47   struct histogram *hist = UP_CAST (s, struct histogram, parent);
  48   gsl_histogram *gslh = hist->gsl_hist;
  49
  50   /* Include cases which are just on the boundary */
  51   if (y == gsl_histogram_max (gslh))
  52     {
  53       double lower, upper;
  54       gsl_histogram_get_range (gslh, gsl_histogram_bins (gslh)-1, &lower, &upper);
  55       gsl_histogram_accumulate (gslh, lower + (upper - lower)/2.0, c);
  56     }
  57   else
  58     gsl_histogram_accumulate (hist->gsl_hist, y, c);
  59 }
  60
  61 static void
  62 destroy (struct statistic *s)
  63 {
  64   struct histogram *h = UP_CAST (s, struct histogram, parent);
  65   gsl_histogram_free (h->gsl_hist);
  66   free (s);
  67 }
  68
  69
  70 /* Find a bin width which is adapted to the scaling of the x axis
  71 In the example here, the binwidth is half of the tick interval.
  72
  73         binwidth
  74          >   <
  75      |....+....+....+.      .+....|
  76    LOWER  1    2    3     N_TICKS
  77         ^LOWDBL                 ^HIGHDBL
  78
  79 This only works, when the min and max value for the histogram are adapted
  80 such that (max-min) is a multiple of the binwidth. Then the location of the
  81 first bin has to be aligned to the ticks.
  82 */
  83 static int
  84 hist_find_pretty_no_of_bins(double bin_width_in, double min, double max,
  85                             double *adjusted_min, double *adjusted_max)
  86 {
  87   double lower, interval;
  88   int n_ticks;
  89   double binwidth;
  90   int nbins;
  91
  92   chart_get_scale (max, min, &lower, &interval, &n_ticks);
  93
  94   if (bin_width_in >= 2 * interval)
  95     {
  96       binwidth = floor(bin_width_in/interval) * interval;
  97       *adjusted_min = lower;
  98     }
  99   else if (bin_width_in >= 1.5 * interval)
 100     {
 101       binwidth = 1.5 * interval;
 102       if (min < (lower + 0.5 * interval))
 103         *adjusted_min = lower;
 104       else
 105         *adjusted_min = lower + 0.5 * interval;
 106     }
 107   else if (bin_width_in >= interval)
 108     {
 109       binwidth = interval;
 110       *adjusted_min = lower;
 111     }
 112   else if (bin_width_in >= (2.0/3.0 * interval))
 113     {
 114       binwidth = (2.0/3.0 * interval);
 115       if (min >= lower + binwidth)
 116         *adjusted_min = lower + binwidth;
 117       else
 118         *adjusted_min = lower;
 119     }
 120   else
 121     {
 122       int i;
 123       for(i = 2; bin_width_in < interval/i; i++);
 124       binwidth = interval/i;
 125       *adjusted_min = floor((min - lower)/binwidth)*binwidth + lower;
 126     }
 127
 128   nbins = ceil((max-*adjusted_min)/binwidth);
 129   *adjusted_max = nbins*binwidth + *adjusted_min;
 130
 131   return nbins;
 132 }
 133
 134
 135 struct histogram *
 136 histogram_create (double bin_width_in, double min, double max)
 137 {
 138   struct histogram *h;
 139   struct statistic *stat;
 140   int bins;
 141   double adjusted_min, adjusted_max;
 142
 143   if (max == min)
 144     {
 145       msg (MW, _("Not creating histogram because the data contains less than 2 distinct values"));
 146       return NULL;
 147     }
 148
 149   assert (bin_width_in > 0);
 150
 151   bins = hist_find_pretty_no_of_bins(bin_width_in, min, max, &adjusted_min, &adjusted_max);
 152
 153   h = xmalloc (sizeof *h);
 154
 155   h->gsl_hist = gsl_histogram_alloc (bins);
 156
 157   /* The bin ranges could be computed with gsl_histogram_set_ranges_uniform,
 158      but the number of bins is adapted to the ticks of the axis such that for example
 159      data ranging from 1.0 to 7.0 with 6 bins will have bin limits at
 160      2.0,3.0,4.0 and 5.0. Due to numerical accuracy the computed bin limits might
 161      be 4.99999999 for a value which is expected to be 5.0. But the limits of
 162      the histogram bins should be that what is expected from the displayed ticks.
 163      Therefore the bin limits are computed from the rounded values which is similar
 164      to the procedure at the chart_get_ticks_format. Actual bin limits should be
 165      exactly what is displayed at the ticks.
 166      But... I cannot reproduce the problem that I see with gsl_histogram_set_ranges_uniform
 167      with the code below without rounding...
 168   */
 169   {
 170     double *ranges = xmalloc (sizeof (double) * (bins + 1));
 171     double interval = (adjusted_max - adjusted_min) / bins;
 172     for (int i = 0; i < bins; i++)
 173       ranges[i] = adjusted_min + interval * i;
 174     ranges[bins] = adjusted_max;
 175     gsl_histogram_set_ranges (h->gsl_hist, ranges, bins + 1);
 176     free (ranges);
 177   }
 178
 179   stat = &h->parent;
 180   stat->accumulate = acc;
 181   stat->destroy = destroy;
 182
 183   return h;
 184 }
 185