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