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