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