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
  49   gsl_histogram_accumulate (hist->gsl_hist, y, c);
  50 }
  51
  52 static void
  53 destroy (struct statistic *s)
  54 {
  55   struct histogram *h = UP_CAST (s, struct histogram, parent);
  56   gsl_histogram_free (h->gsl_hist);
  57   free (s);
  58 }
  59
  60
  61 /* Find a bin width which is adapted to the scaling of the x axis
  62 In the example here, the binwidth is half of the tick interval.
  63
  64         binwidth
  65          >   <
  66      |....+....+....+.      .+....|
  67    LOWER  1    2    3     N_TICKS
  68         ^LOWDBL                 ^HIGHDBL
  69
  70 This only works, when the min and max value for the histogram are adapted
  71 such that (max-min) is a multiple of the binwidth. Then the location of the
  72 first bin has to be aligned to the ticks.
  73 */
  74 static int
  75 hist_find_pretty_no_of_bins(double bin_width_in, double min, double max,
  76                             double *adjusted_min, double *adjusted_max)
  77 {
  78   double lower, interval;
  79   int n_ticks;
  80   double binwidth;
  81   int nbins;
  82
  83   chart_get_scale (max, min, &lower, &interval, &n_ticks);
  84
  85   if (bin_width_in >= 2 * interval)
  86     {
  87       binwidth = floor(bin_width_in/interval) * interval;
  88       *adjusted_min = lower;
  89     }
  90   else if (bin_width_in >= 1.5 * interval)
  91     {
  92       binwidth = 1.5 * interval;
  93       if (min < (lower + 0.5 * interval))
  94         *adjusted_min = lower;
  95       else
  96         *adjusted_min = lower + 0.5 * interval;
  97     }
  98   else if (bin_width_in >= interval)
  99     {
 100       binwidth = interval;
 101       *adjusted_min = lower;
 102     }
 103   else if (bin_width_in >= (2.0/3.0 * interval))
 104     {
 105       binwidth = (2.0/3.0 * interval);
 106       if (min >= lower + binwidth)
 107         *adjusted_min = lower + binwidth;
 108       else
 109         *adjusted_min = lower;
 110     }
 111   else
 112     {
 113       int i;
 114       for(i = 2; bin_width_in < interval/i; i++);
 115       binwidth = interval/i;
 116       *adjusted_min = floor((min - lower)/binwidth)*binwidth + lower;
 117     }
 118
 119   nbins = ceil((max-*adjusted_min)/binwidth);
 120   *adjusted_max = nbins*binwidth + *adjusted_min;
 121
 122   return nbins;
 123 }
 124
 125
 126 struct histogram *
 127 histogram_create (double bin_width_in, double min, double max)
 128 {
 129   struct histogram *h;
 130   struct statistic *stat;
 131   int bins;
 132   double adjusted_min, adjusted_max;
 133
 134   if (max == min)
 135     {
 136       msg (MW, _("Not creating histogram because the data contains less than 2 distinct values"));
 137       return NULL;
 138     }
 139
 140   assert (bin_width_in > 0);
 141
 142   bins = hist_find_pretty_no_of_bins(bin_width_in, min, max, &adjusted_min, &adjusted_max);
 143
 144   h = xmalloc (sizeof *h);
 145
 146   h->gsl_hist = gsl_histogram_alloc (bins);
 147
 148   gsl_histogram_set_ranges_uniform (h->gsl_hist, adjusted_min, adjusted_max);
 149
 150   stat = &h->parent;
 151   stat->accumulate = acc;
 152   stat->destroy = destroy;
 153
 154   return h;
 155 }
 156