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   gsl_histogram_set_ranges_uniform (h->gsl_hist, adjusted_min, adjusted_max);
 158
 159   stat = &h->parent;
 160   stat->accumulate = acc;
 161   stat->destroy = destroy;
 162
 163   return h;
 164 }
 165