85fb3700e6985b41b7878eddeb37b046035a001b
[pspp-builds.git] / src / plot-hist.c
1 /* PSPP - computes sample statistics.
2    Copyright (C) 2004 Free Software Foundation, Inc.
3    Written by John Darrington <john@darrington.wattle.id.au>
4
5    This program is free software; you can redistribute it and/or
6    modify it under the terms of the GNU General Public License as
7    published by the Free Software Foundation; either version 2 of the
8    License, or (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful, but
11    WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13    General Public License for more details.
14
15    You should have received a copy of the GNU General Public License
16    along with this program; if not, write to the Free Software
17    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
18    02110-1301, USA. */
19
20 /* If you add/modify any public symbols in this file, don't forget to
21    change the stubs in dummy-chart.c */
22
23 #include <config.h>
24
25 #include <stdio.h>
26 #include <plot.h>
27 #include <math.h>
28 #include <gsl/gsl_histogram.h>
29 #include <gsl/gsl_randist.h>
30 #include <assert.h>
31
32 #include "hash.h"
33
34 #include "var.h"
35 #include "chart.h"
36
37 /* Write the legend of the chart */
38 void
39 histogram_write_legend(struct chart *ch, const struct normal_curve *norm)
40 {
41   char buf[100];
42   if ( !ch )
43           return ;
44
45   pl_savestate_r(ch->lp);
46
47   sprintf(buf,"N = %.2f",norm->N);
48   pl_move_r(ch->lp, ch->legend_left, ch->data_bottom);
49   pl_alabel_r(ch->lp,0,'b',buf);
50
51   sprintf(buf,"Mean = %.1f",norm->mean);
52   pl_fmove_r(ch->lp,ch->legend_left,ch->data_bottom + ch->font_size * 1.5);
53   pl_alabel_r(ch->lp,0,'b',buf);
54
55   sprintf(buf,"Std. Dev = %.2f",norm->stddev);
56   pl_fmove_r(ch->lp,ch->legend_left,ch->data_bottom + ch->font_size * 1.5 * 2);
57   pl_alabel_r(ch->lp,0,'b',buf);
58
59   pl_restorestate_r(ch->lp);    
60 }
61
62 static void hist_draw_bar(struct chart *ch, const gsl_histogram *hist, int bar);
63
64
65 static void
66 hist_draw_bar(struct chart *ch, const gsl_histogram *hist, int bar)
67 {
68   if ( !ch ) 
69     return ;
70
71
72   {
73     double upper;
74     double lower;
75     double height;
76
77     const size_t bins = gsl_histogram_bins(hist);
78     const double x_pos = (ch->data_right - ch->data_left) * bar / (double) bins ;
79     const double width = (ch->data_right - ch->data_left) / (double) bins ;
80
81
82     assert ( 0 == gsl_histogram_get_range(hist, bar, &lower, &upper));
83
84     assert( upper >= lower);
85
86     height = gsl_histogram_get(hist, bar) * 
87       (ch->data_top - ch->data_bottom) / gsl_histogram_max_val(hist);
88
89     pl_savestate_r(ch->lp);
90     pl_move_r(ch->lp,ch->data_left, ch->data_bottom);
91     pl_fillcolorname_r(ch->lp, ch->fill_colour); 
92     pl_filltype_r(ch->lp,1);
93
94
95     pl_fboxrel_r(ch->lp,
96                  x_pos, 0,
97                  x_pos + width, height);
98
99     pl_restorestate_r(ch->lp);
100
101     {
102       char buf[5];
103       snprintf(buf,5,"%g",(upper + lower) / 2.0);
104       draw_tick(ch, TICK_ABSCISSA,
105                 x_pos + width / 2.0, buf);
106     }
107   }
108 }
109
110
111
112
113 void
114 histogram_plot(const gsl_histogram *hist,
115                const char *factorname,
116                const struct normal_curve *norm, short show_normal)
117 {
118   int i;
119   int bins;
120   
121   struct chart *ch;
122
123   bins = gsl_histogram_bins(hist);
124
125   ch = chart_create();
126   chart_write_title(ch, _("HISTOGRAM"));
127
128   chart_write_ylabel(ch, _("Frequency"));
129   chart_write_xlabel(ch, factorname);
130
131   chart_write_yscale(ch, 0, gsl_histogram_max_val(hist), 5);
132
133   for ( i = 0 ; i < bins ; ++i ) 
134       hist_draw_bar(ch, hist, i);
135
136   histogram_write_legend(ch, norm);
137
138   if ( show_normal  )
139   {
140     /* Draw the normal curve */    
141
142     double d ;
143     double x_min, x_max, not_used ;
144     double abscissa_scale ;
145     double ordinate_scale ;
146     double range ;
147
148     gsl_histogram_get_range(hist, 0, &x_min, &not_used);
149     range = not_used - x_min;
150     gsl_histogram_get_range(hist, bins - 1, &not_used, &x_max);
151     assert(range == x_max - not_used);
152
153     abscissa_scale = (ch->data_right - ch->data_left) / (x_max - x_min);
154     ordinate_scale = (ch->data_top - ch->data_bottom) / 
155       gsl_histogram_max_val(hist) ;
156
157     pl_move_r(ch->lp, ch->data_left, ch->data_bottom);    
158     for( d = ch->data_left; 
159          d <= ch->data_right ; 
160          d += (ch->data_right - ch->data_left) / 100.0)
161       {    
162         const double x = (d - ch->data_left) / abscissa_scale + x_min ; 
163         const double y = norm->N * range * 
164           gsl_ran_gaussian_pdf(x - norm->mean, norm->stddev);
165
166         pl_fcont_r(ch->lp,  d,  ch->data_bottom  + y * ordinate_scale);
167
168       }
169     pl_endpath_r(ch->lp);
170
171   }
172   chart_submit(ch);
173 }
174