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_fillcolor_r (lp,
 100                   geom->fill_colour.red * 257,
 101                   geom->fill_colour.green * 257,
 102                   geom->fill_colour.blue * 257);
 103   pl_filltype_r (lp,1);
 104
 105
 106   pl_fboxrel_r (lp,
 107                 x_pos, 0,
 108                 x_pos + width, height);
 109
 110   pl_restorestate_r (lp);
 111
 112   draw_tick (lp, geom, TICK_ABSCISSA,
 113              x_pos + width / 2.0, "%g", (upper + lower) / 2.0);
 114 }
 115
 116 struct histogram_chart
 117   {
 118     struct chart chart;
 119     gsl_histogram *gsl_hist;
 120     char *label;
 121     double n;
 122     double mean;
 123     double stddev;
 124     bool show_normal;
 125   };
 126
 127 /* Plots a histogram of the data in HIST with the given LABEL.
 128    Labels the histogram with each of N, MEAN, and STDDEV that is
 129    not SYSMIS.  If all three are not SYSMIS and SHOW_NORMAL is
 130    true, also draws a normal curve on the histogram. */
 131 struct chart *
 132 histogram_chart_create (const struct histogram *hist, const char *label,
 133                         double n, double mean, double stddev,
 134                         bool show_normal)
 135 {
 136   struct histogram_chart *h;
 137
 138   h = xmalloc (sizeof *h);
 139   chart_init (&h->chart, &histogram_chart_class);
 140   h->gsl_hist = hist->gsl_hist ? gsl_histogram_clone (hist->gsl_hist) : NULL;
 141   h->label = xstrdup (label);
 142   h->n = n;
 143   h->mean = mean;
 144   h->stddev = stddev;
 145   h->show_normal = show_normal;
 146   return &h->chart;
 147 }
 148
 149 static void
 150 histogram_chart_draw (const struct chart *chart, plPlotter *lp)
 151 {
 152   struct histogram_chart *h = (struct histogram_chart *) chart;
 153   struct chart_geometry geom;
 154   int i;
 155   int bins;
 156
 157   chart_geometry_init (lp, &geom);
 158
 159   chart_write_title (lp, &geom, _("HISTOGRAM"));
 160
 161   chart_write_ylabel (lp, &geom, _("Frequency"));
 162   chart_write_xlabel (lp, &geom, h->label);
 163
 164   if (h->gsl_hist == NULL)
 165     {
 166       /* Probably all values are SYSMIS. */
 167       return;
 168     }
 169
 170   bins = gsl_histogram_bins (h->gsl_hist);
 171
 172   chart_write_yscale (lp, &geom, 0, gsl_histogram_max_val (h->gsl_hist), 5);
 173
 174   for (i = 0; i < bins; i++)
 175     hist_draw_bar (lp, &geom, h->gsl_hist, i);
 176
 177   histogram_write_legend (lp, &geom, h->n, h->mean, h->stddev);
 178
 179   if (h->show_normal
 180       && h->n != SYSMIS && h->mean != SYSMIS && h->stddev != SYSMIS)
 181     {
 182       /* Draw the normal curve */
 183       double d;
 184       double x_min, x_max, not_used;
 185       double abscissa_scale;
 186       double ordinate_scale;
 187       double range;
 188
 189       gsl_histogram_get_range (h->gsl_hist, 0, &x_min, &not_used);
 190       range = not_used - x_min;
 191       gsl_histogram_get_range (h->gsl_hist, bins - 1, &not_used, &x_max);
 192
 193       abscissa_scale = (geom.data_right - geom.data_left) / (x_max - x_min);
 194       ordinate_scale = (geom.data_top - geom.data_bottom) /
 195         gsl_histogram_max_val (h->gsl_hist);
 196
 197       pl_move_r (lp, geom.data_left, geom.data_bottom);
 198       for (d = geom.data_left;
 199            d <= geom.data_right;
 200            d += (geom.data_right - geom.data_left) / 100.0)
 201         {
 202           const double x = (d - geom.data_left) / abscissa_scale + x_min;
 203           const double y = h->n * range *
 204             gsl_ran_gaussian_pdf (x - h->mean, h->stddev);
 205
 206           pl_fcont_r (lp,  d,  geom.data_bottom  + y * ordinate_scale);
 207
 208         }
 209       pl_endpath_r (lp);
 210     }
 211
 212   chart_geometry_free (lp);
 213 }
 214
 215
 216 static void
 217 histogram_chart_destroy (struct chart *chart)
 218 {
 219   struct histogram_chart *h = (struct histogram_chart *) chart;
 220   if (h->gsl_hist != NULL)
 221     gsl_histogram_free (h->gsl_hist);
 222   free (h->label);
 223   free (h);
 224 }
 225
 226 static const struct chart_class histogram_chart_class =
 227   {
 228     histogram_chart_draw,
 229     histogram_chart_destroy
 230   };