New module to perform decimal floating point arithmetic for charts.
[pspp] / src / output / charts / plot-hist-cairo.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2009, 2011, 2014 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/decimal.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
29 #include "gettext.h"
30 #define _(msgid) gettext (msgid)
31
32 /* Write the legend of the chart */
33 static void
34 histogram_write_legend (cairo_t *cr, const struct xrchart_geometry *geom,
35                         double n, double mean, double stddev)
36 {
37   double y = geom->axis[SCALE_ORDINATE].data_min;
38   cairo_save (cr);
39
40   if (n != SYSMIS)
41     {
42       char *buf = xasprintf (_("N = %.2f"), n);
43       cairo_move_to (cr, geom->legend_left, y);
44       xrchart_label (cr, 'l', 'b', geom->font_size, buf);
45       y += geom->font_size * 1.5;
46       free (buf);
47     }
48
49   if (mean != SYSMIS)
50     {
51       char *buf = xasprintf (_("Mean = %.1f"), mean);
52       cairo_move_to (cr,geom->legend_left, y);
53       xrchart_label (cr, 'l', 'b', geom->font_size, buf);
54       y += geom->font_size * 1.5;
55       free (buf);
56     }
57
58   if (stddev != SYSMIS)
59     {
60       char *buf = xasprintf (_("Std. Dev = %.2f"), stddev);
61       cairo_move_to (cr, geom->legend_left, y);
62       xrchart_label (cr, 'l', 'b', geom->font_size, buf);
63       free (buf);
64     }
65
66   cairo_restore (cr);
67 }
68
69 static void
70 hist_draw_bar (cairo_t *cr, const struct xrchart_geometry *geom,
71                const gsl_histogram *h, int bar, bool label)
72 {
73   double upper;
74   double lower;
75   double height;
76
77   const size_t bins = gsl_histogram_bins (h);
78
79   const double x_pos =
80     (geom->axis[SCALE_ABSCISSA].data_max - geom->axis[SCALE_ABSCISSA].data_min) *
81     bar / (double) bins ;
82
83   const double width =
84     (geom->axis[SCALE_ABSCISSA].data_max - geom->axis[SCALE_ABSCISSA].data_min) / (double) bins ;
85
86   assert ( 0 == gsl_histogram_get_range (h, bar, &lower, &upper));
87
88   assert ( upper >= lower);
89
90   height = geom->axis[SCALE_ORDINATE].scale * gsl_histogram_get (h, bar);
91
92   cairo_rectangle (cr,
93                    geom->axis[SCALE_ABSCISSA].data_min + x_pos,
94                    geom->axis[SCALE_ORDINATE].data_min,
95                    width, height);
96   cairo_save (cr);
97   cairo_set_source_rgb (cr,
98                         geom->fill_colour.red / 255.0,
99                         geom->fill_colour.green / 255.0,
100                         geom->fill_colour.blue / 255.0);
101   cairo_fill_preserve (cr);
102   cairo_restore (cr);
103   cairo_stroke (cr);
104
105   if (label)
106     {
107       struct decimal decupper;
108       struct decimal declower;
109       struct decimal middle;
110       decimal_from_double (&declower, lower);
111       decimal_from_double (&decupper, upper);
112       middle = declower;
113       decimal_add (&middle, &decupper);
114       decimal_int_divide (&middle, 2);
115       char *str = decimal_to_string (&middle);
116       draw_tick (cr, geom, SCALE_ABSCISSA, bins > 10,
117                  x_pos + width / 2.0, "%s", str);
118       free (str);
119     }
120 }
121
122 void
123 xrchart_draw_histogram (const struct chart_item *chart_item, cairo_t *cr,
124                         struct xrchart_geometry *geom)
125 {
126   struct histogram_chart *h = to_histogram_chart (chart_item);
127   int i;
128   int bins;
129
130   xrchart_write_title (cr, geom, _("HISTOGRAM"));
131
132   xrchart_write_ylabel (cr, geom, _("Frequency"));
133   xrchart_write_xlabel (cr, geom, chart_item_get_title (chart_item));
134
135   if (h->gsl_hist == NULL)
136     {
137       /* Probably all values are SYSMIS. */
138       return;
139     }
140
141   bins = gsl_histogram_bins (h->gsl_hist);
142
143   xrchart_write_yscale (cr, geom, 0, gsl_histogram_max_val (h->gsl_hist));
144
145   for (i = 0; i < bins; i++)
146     {
147       hist_draw_bar (cr, geom, h->gsl_hist, i, true);
148     }
149
150   histogram_write_legend (cr, geom, h->n, h->mean, h->stddev);
151
152   if (h->show_normal
153       && h->n != SYSMIS && h->mean != SYSMIS && h->stddev != SYSMIS)
154     {
155       /* Draw the normal curve */
156       double d;
157       double x_min, x_max, not_used;
158       double abscissa_scale;
159       double ordinate_scale;
160       double range;
161
162       gsl_histogram_get_range (h->gsl_hist, 0, &x_min, &not_used);
163       range = not_used - x_min;
164       gsl_histogram_get_range (h->gsl_hist, bins - 1, &not_used, &x_max);
165
166       abscissa_scale = (geom->axis[SCALE_ABSCISSA].data_max - geom->axis[SCALE_ABSCISSA].data_min) / (x_max - x_min);
167       ordinate_scale = (geom->axis[SCALE_ORDINATE].data_max - geom->axis[SCALE_ORDINATE].data_min) /
168         gsl_histogram_max_val (h->gsl_hist);
169
170       cairo_move_to (cr, geom->axis[SCALE_ABSCISSA].data_min, geom->axis[SCALE_ORDINATE].data_min);
171       for (d = geom->axis[SCALE_ABSCISSA].data_min;
172            d <= geom->axis[SCALE_ABSCISSA].data_max;
173            d += (geom->axis[SCALE_ABSCISSA].data_max - geom->axis[SCALE_ABSCISSA].data_min) / 100.0)
174         {
175           const double x = (d - geom->axis[SCALE_ABSCISSA].data_min) / abscissa_scale + x_min;
176           const double y = h->n * range *
177             gsl_ran_gaussian_pdf (x - h->mean, h->stddev);
178
179           cairo_line_to (cr, d, geom->axis[SCALE_ORDINATE].data_min  + y * ordinate_scale);
180
181         }
182       cairo_stroke (cr);
183     }
184 }