src/output/charts/plot-hist.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright  (C) 2004, 2009 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
  18 #include <config.h>
  19
  20 #include <stdio.h>
  21 #include <math.h>
  22 #include <gsl/gsl_histogram.h>
  23 #include <gsl/gsl_randist.h>
  24 #include <assert.h>
  25
  26 #include <output/charts/plot-hist.h>
  27 #include <output/charts/plot-chart.h>
  28 #include <output/chart-provider.h>
  29
  30 #include <data/variable.h>
  31 #include <libpspp/hash.h>
  32 #include <output/chart.h>
  33 #include <math/histogram.h>
  34 #include <math/moments.h>
  35
  36 #include "gettext.h"
  37 #define _(msgid) gettext (msgid)
  38
  39 static const struct chart_class histogram_chart_class;
  40
  41 /* Write the legend of the chart */
  42 static void
  43 histogram_write_legend (plPlotter *lp, const struct chart_geometry *geom,
  44                         double n, double mean, double stddev)
  45 {
  46   double y = geom->data_bottom;
  47   pl_savestate_r (lp);
  48
  49   if (n != SYSMIS)
  50     {
  51       char *buf = xasprintf ("N = %.2f", n);
  52       pl_move_r (lp, geom->legend_left, y);
  53       pl_alabel_r (lp, 0, 'b', buf);
  54       y += geom->font_size * 1.5;
  55       free (buf);
  56     }
  57
  58   if (mean != SYSMIS)
  59     {
  60       char *buf = xasprintf ("Mean = %.1f", mean);
  61       pl_fmove_r (lp,geom->legend_left, y);
  62       pl_alabel_r (lp, 0, 'b', buf);
  63       y += geom->font_size * 1.5;
  64       free (buf);
  65     }
  66
  67   if (stddev != SYSMIS)
  68     {
  69       char *buf = xasprintf ("Std. Dev = %.2f", stddev);
  70       pl_fmove_r (lp, geom->legend_left, y);
  71       pl_alabel_r (lp, 0, 'b', buf);
  72       free (buf);
  73     }
  74
  75   pl_restorestate_r (lp);
  76 }
  77
  78 static void
  79 hist_draw_bar (plPlotter *lp, const struct chart_geometry *geom,
  80                const gsl_histogram *h, int bar)
  81 {
  82   double upper;
  83   double lower;
  84   double height;
  85
  86   const size_t bins = gsl_histogram_bins (h);
  87   const double x_pos = (geom->data_right - geom->data_left) * bar / (double) bins ;
  88   const double width = (geom->data_right - geom->data_left) / (double) bins ;
  89
  90   assert ( 0 == gsl_histogram_get_range (h, bar, &lower, &upper));
  91
  92   assert ( upper >= lower);
  93
  94   height = gsl_histogram_get (h, bar) *
  95     (geom->data_top - geom->data_bottom) / gsl_histogram_max_val (h);
  96
  97   pl_savestate_r (lp);
  98   pl_move_r (lp,geom->data_left, geom->data_bottom);
  99   pl_fillcolorname_r (lp, geom->fill_colour);
 100   pl_filltype_r (lp,1);
 101
 102
 103   pl_fboxrel_r (lp,
 104                 x_pos, 0,
 105                 x_pos + width, height);
 106
 107   pl_restorestate_r (lp);
 108
 109   draw_tick (lp, geom, TICK_ABSCISSA,
 110              x_pos + width / 2.0, "%g", (upper + lower) / 2.0);
 111 }
 112
 113 struct histogram_chart
 114   {
 115     struct chart chart;
 116     gsl_histogram *gsl_hist;
 117     char *label;
 118     double n;
 119     double mean;
 120     double stddev;
 121     bool show_normal;
 122   };
 123
 124 /* Plots a histogram of the data in HIST with the given LABEL.
 125    Labels the histogram with each of N, MEAN, and STDDEV that is
 126    not SYSMIS.  If all three are not SYSMIS and SHOW_NORMAL is
 127    true, also draws a normal curve on the histogram. */
 128 struct chart *
 129 histogram_chart_create (const struct histogram *hist, const char *label,
 130                         double n, double mean, double stddev,
 131                         bool show_normal)
 132 {
 133   struct histogram_chart *h;
 134
 135   h = xmalloc (sizeof *h);
 136   chart_init (&h->chart, &histogram_chart_class);
 137   h->gsl_hist = hist->gsl_hist ? gsl_histogram_clone (hist->gsl_hist) : NULL;
 138   h->label = xstrdup (label);
 139   h->n = n;
 140   h->mean = mean;
 141   h->stddev = stddev;
 142   h->show_normal = show_normal;
 143   return &h->chart;
 144 }
 145
 146 static void
 147 histogram_chart_draw (const struct chart *chart, plPlotter *lp)
 148 {
 149   struct histogram_chart *h = (struct histogram_chart *) chart;
 150   struct chart_geometry geom;
 151   int i;
 152   int bins;
 153
 154   chart_geometry_init (lp, &geom);
 155
 156   chart_write_title (lp, &geom, _("HISTOGRAM"));
 157
 158   chart_write_ylabel (lp, &geom, _("Frequency"));
 159   chart_write_xlabel (lp, &geom, h->label);
 160
 161   if (h->gsl_hist == NULL)
 162     {
 163       /* Probably all values are SYSMIS. */
 164       return;
 165     }
 166
 167   bins = gsl_histogram_bins (h->gsl_hist);
 168
 169   chart_write_yscale (lp, &geom, 0, gsl_histogram_max_val (h->gsl_hist), 5);
 170
 171   for (i = 0; i < bins; i++)
 172     hist_draw_bar (lp, &geom, h->gsl_hist, i);
 173
 174   histogram_write_legend (lp, &geom, h->n, h->mean, h->stddev);
 175
 176   if (h->show_normal
 177       && h->n != SYSMIS && h->mean != SYSMIS && h->stddev != SYSMIS)
 178     {
 179       /* Draw the normal curve */
 180       double d;
 181       double x_min, x_max, not_used;
 182       double abscissa_scale;
 183       double ordinate_scale;
 184       double range;
 185
 186       gsl_histogram_get_range (h->gsl_hist, 0, &x_min, &not_used);
 187       range = not_used - x_min;
 188       gsl_histogram_get_range (h->gsl_hist, bins - 1, &not_used, &x_max);
 189
 190       abscissa_scale = (geom.data_right - geom.data_left) / (x_max - x_min);
 191       ordinate_scale = (geom.data_top - geom.data_bottom) /
 192         gsl_histogram_max_val (h->gsl_hist);
 193
 194       pl_move_r (lp, geom.data_left, geom.data_bottom);
 195       for (d = geom.data_left;
 196            d <= geom.data_right;
 197            d += (geom.data_right - geom.data_left) / 100.0)
 198         {
 199           const double x = (d - geom.data_left) / abscissa_scale + x_min;
 200           const double y = h->n * range *
 201             gsl_ran_gaussian_pdf (x - h->mean, h->stddev);
 202
 203           pl_fcont_r (lp,  d,  geom.data_bottom  + y * ordinate_scale);
 204
 205         }
 206       pl_endpath_r (lp);
 207     }
 208
 209   chart_geometry_free (lp);
 210 }
 211
 212
 213 static void
 214 histogram_chart_destroy (struct chart *chart)
 215 {
 216   struct histogram_chart *h = (struct histogram_chart *) chart;
 217   if (h->gsl_hist != NULL)
 218     gsl_histogram_free (h->gsl_hist);
 219   free (h->label);
 220   free (h);
 221 }
 222
 223 static const struct chart_class histogram_chart_class =
 224   {
 225     histogram_chart_draw,
 226     histogram_chart_destroy
 227   };