f8663428096ff06e434d24f3d0976b880cdf57b0
[pspp] / src / output / charts / plot-hist-cairo.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2009, 2011, 2014, 2015 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 #include "math/chart-geometry.h"
19 #include "output/charts/plot-hist.h"
20
21 #include <float.h>
22 #include <gsl/gsl_randist.h>
23
24 #include "data/val-type.h"
25 #include "output/cairo-chart.h"
26
27 #include "gl/xvasprintf.h"
28 #include "gl/minmax.h"
29
30 #include "gettext.h"
31 #define _(msgid) gettext (msgid)
32
33 /* Write the legend of the chart */
34 static void
35 histogram_write_legend (cairo_t *cr, const struct xrchart_geometry *geom,
36                         double n, double mean, double stddev)
37 {
38   double y = geom->axis[SCALE_ORDINATE].data_min;
39   cairo_save (cr);
40
41   if (n != SYSMIS)
42     {
43       char *buf = xasprintf (_("N = %.2f"), n);
44       cairo_move_to (cr, geom->legend_left, y);
45       xrchart_label (cr, 'l', 'b', geom->font_size, buf);
46       y += geom->font_size * 1.5;
47       free (buf);
48     }
49
50   if (mean != SYSMIS)
51     {
52       char *buf = xasprintf (_("Mean = %.1f"), mean);
53       cairo_move_to (cr,geom->legend_left, y);
54       xrchart_label (cr, 'l', 'b', geom->font_size, buf);
55       y += geom->font_size * 1.5;
56       free (buf);
57     }
58
59   if (stddev != SYSMIS)
60     {
61       char *buf = xasprintf (_("Std. Dev = %.2f"), stddev);
62       cairo_move_to (cr, geom->legend_left, y);
63       xrchart_label (cr, 'l', 'b', geom->font_size, buf);
64       free (buf);
65     }
66
67   cairo_restore (cr);
68 }
69
70 static void
71 hist_draw_bar (cairo_t *cr, const struct xrchart_geometry *geom,
72                const gsl_histogram *h, int bar)
73 {
74   double upper;
75   double lower;
76   double height;
77
78   assert (0 == gsl_histogram_get_range (h, bar, &lower, &upper));
79   assert (upper >= lower);
80
81   const double x_pos =
82     (lower - geom->axis[SCALE_ABSCISSA].min) * geom->axis[SCALE_ABSCISSA].scale
83     +  geom->axis[SCALE_ABSCISSA].data_min;
84   const double width = (upper - lower) * geom->axis[SCALE_ABSCISSA].scale;
85
86   height = geom->axis[SCALE_ORDINATE].scale * gsl_histogram_get (h, bar);
87
88   cairo_rectangle (cr,
89                    x_pos,
90                    geom->axis[SCALE_ORDINATE].data_min,
91                    width, height);
92   cairo_save (cr);
93   cairo_set_source_rgb (cr,
94                         geom->fill_colour.red / 255.0,
95                         geom->fill_colour.green / 255.0,
96                         geom->fill_colour.blue / 255.0);
97   cairo_fill_preserve (cr);
98   cairo_restore (cr);
99   cairo_stroke (cr);
100 }
101
102 void
103 xrchart_draw_histogram (const struct chart *chart, cairo_t *cr,
104                         struct xrchart_geometry *geom)
105 {
106   struct histogram_chart *h = to_histogram_chart (chart);
107   int i;
108   int bins;
109
110   xrchart_write_title (cr, geom, _("HISTOGRAM"));
111
112   xrchart_write_ylabel (cr, geom, _("Frequency"));
113   xrchart_write_xlabel (cr, geom, chart_get_title (chart));
114
115   if (h->gsl_hist == NULL)
116     {
117       /* Probably all values are SYSMIS. */
118       return;
119     }
120
121   if (! xrchart_write_yscale (cr, geom, 0, gsl_histogram_max_val (h->gsl_hist)))
122     return;
123   if (! xrchart_write_xscale (cr, geom, gsl_histogram_min (h->gsl_hist),
124                               gsl_histogram_max (h->gsl_hist)))
125     return;
126
127
128   /* Draw the ticks and compute if the rendered tick text is wider than the bin */
129   bins = gsl_histogram_bins (h->gsl_hist);
130
131   for (i = 0; i < bins; i++)
132     {
133       hist_draw_bar (cr, geom, h->gsl_hist, i);
134     }
135
136   histogram_write_legend (cr, geom, h->n, h->mean, h->stddev);
137
138   if (h->show_normal
139       && h->n != SYSMIS && h->mean != SYSMIS && h->stddev != SYSMIS)
140     {
141       /* Draw the normal curve */
142       double x_min, x_max;
143       double ordinate_scale;
144       double binwidth;
145       double x;
146
147       gsl_histogram_get_range (h->gsl_hist, 0, &x_min, &x_max);
148       binwidth = x_max - x_min;
149
150       /* The integral over the histogram is binwidth * sum(bin_i), while the integral over the pdf is 1 */
151       /* Therefore the pdf has to be scaled accordingly such that the integrals are equal               */
152       ordinate_scale = binwidth * gsl_histogram_sum(h->gsl_hist);
153
154       /* Clip normal curve to the rectangle formed by the axes. */
155       cairo_save (cr);
156       cairo_rectangle (cr, geom->axis[SCALE_ABSCISSA].data_min, geom->axis[SCALE_ORDINATE].data_min,
157                        geom->axis[SCALE_ABSCISSA].data_max - geom->axis[SCALE_ABSCISSA].data_min,
158                        geom->axis[SCALE_ORDINATE].data_max - geom->axis[SCALE_ORDINATE].data_min);
159       cairo_clip (cr);
160
161       cairo_move_to (cr, geom->axis[SCALE_ABSCISSA].data_min, geom->axis[SCALE_ORDINATE].data_min);
162       for (x = geom->axis[SCALE_ABSCISSA].min;
163            x <= geom->axis[SCALE_ABSCISSA].max;
164            x += (geom->axis[SCALE_ABSCISSA].max - geom->axis[SCALE_ABSCISSA].min) / 100.0)
165         {
166           const double y = gsl_ran_gaussian_pdf (x - h->mean, h->stddev) * ordinate_scale;
167           /* Transform to drawing coordinates */
168           const double x_pos = (x - geom->axis[SCALE_ABSCISSA].min) * geom->axis[SCALE_ABSCISSA].scale + geom->axis[SCALE_ABSCISSA].data_min;
169           const double y_pos = (y - geom->axis[SCALE_ORDINATE].min) * geom->axis[SCALE_ORDINATE].scale + geom->axis[SCALE_ORDINATE].data_min;
170           cairo_line_to (cr, x_pos, y_pos);
171         }
172       cairo_stroke (cr);
173
174       cairo_restore (cr);
175     }
176 }