gui: Allow File|Open to select an encoding for system files.
[pspp] / src / math / histogram.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2004, 2008, 2009, 2011, 2012 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
19 #include "math/histogram.h"
20
21 #include <gsl/gsl_histogram.h>
22 #include <math.h>
23
24 #include "data/settings.h"
25 #include "libpspp/message.h"
26 #include "libpspp/assertion.h"
27 #include "libpspp/cast.h"
28 #include "math/chart-geometry.h"
29
30 #include "gettext.h"
31 #define _(msgid) gettext (msgid)
32 #define N_(msgid) msgid
33
34
35 #include "gl/xalloc.h"
36
37 void
38 histogram_add (struct histogram *h, double y, double c)
39 {
40   struct statistic *stat = &h->parent;
41   stat->accumulate (stat, NULL, c, 0, y);
42 }
43
44 static void
45 acc (struct statistic *s, const struct ccase *cx UNUSED, double c, double cc UNUSED, double y)
46 {
47   struct histogram *hist = UP_CAST (s, struct histogram, parent);
48
49   gsl_histogram_accumulate (hist->gsl_hist, y, c);
50 }
51
52 static void
53 destroy (struct statistic *s)
54 {
55   struct histogram *h = UP_CAST (s, struct histogram, parent);
56   gsl_histogram_free (h->gsl_hist);
57   free (s);
58 }
59
60
61 static
62 double get_slack (double limit, double half_bin_width, int *n_half_bins)
63 {
64   double ipart, remainder;
65
66   assert (half_bin_width > 0);
67
68   remainder =  modf (limit / half_bin_width, &ipart);
69
70   /* In C modf and % behave in an unexpected (to me at any rate) manner
71      when presented with a negative value
72
73      For example, modf (-7.0 / 3.0) returns -2.0 R -0.3333
74   */
75
76   
77   *n_half_bins = ipart;
78
79   return remainder * half_bin_width;
80 }
81
82
83 /* This functions adjusts the upper and lower range of the histogram to make them fit BIN_WIDTH
84    MIN and MAX are the lowest and highest data to be plotted in the histogram.
85    ADJ_MIN and ADJ_MAX are locations of the adjusted values of MIN and MAX (the range will
86    always be  equal or slightly larger).
87    Returns the number of bins.
88
89    The "testing_assert" expressions in this function should be algebraically correct.
90    However, due to floating point rounding they could fail, especially when small numbers
91    are involved.  In normal use, therefore, testing_assert does nothing.
92  */
93 static int
94 adjust_bin_ranges (double bin_width, double min, double max, double *adj_min, double *adj_max)
95 {
96   const double half_bin_width = bin_width / 2.0;
97
98   /* The lower and upper limits of the histogram, in units of half
99      bin widths */
100   int lower_limit, upper_limit;
101
102   double lower_slack =  get_slack (min, half_bin_width, &lower_limit);
103   double upper_slack = -get_slack (max, half_bin_width, &upper_limit);
104
105   testing_assert (max > min);
106
107   /* If min is negative, then lower_slack may be less than zero.
108      In this case, the lower bound must be extended in the negative direction
109      so that it is less than OR EQUAL to min.
110    */
111   if (lower_slack < 0)
112     {
113       lower_limit--;
114       lower_slack += half_bin_width;
115     }
116   testing_assert (lower_limit * half_bin_width <= min);
117
118   /* However, the upper bound must be extended regardless, because histogram bins
119      span the range [lower, upper). In other words, the upper bound must be
120      greater than max.
121   */
122   upper_limit++;;
123   upper_slack += half_bin_width;
124   testing_assert (upper_limit * half_bin_width > max);
125
126   /* The range must be an EVEN number of half bin_widths */
127   if ( (upper_limit - lower_limit) % 2)
128     {
129       /* Extend the range at the end which gives the least unused space */
130       if (upper_slack > lower_slack)
131         {
132           lower_limit--;
133           lower_slack += half_bin_width;
134         }
135       else
136         {
137           upper_limit++;
138           upper_slack += half_bin_width;
139         }
140     }
141
142   /* But the range should be aligned to an ODD number of
143      half bin widths, so that the labels are aesthetically pleasing ones.
144      Otherwise we are likely to get labels such as -3 -1 1 3 instead of -2 0 2 4
145   */
146   if ( lower_limit % 2 == 0)
147     {
148       /* If there is not enough slack at either end to perform a shift,
149          then we must extend the range so that there is.  We must extend
150          by two half bin widths in order to preserve the EVEN condition
151          established above.  Also, we extend on the end with the least
152          slack, in order to keep things as balanced as possible. */
153       if ( upper_slack > lower_slack && upper_slack <= half_bin_width)
154         {
155           lower_limit -= 2;
156           lower_slack += 2 * half_bin_width;
157         }
158            
159       if (lower_slack > upper_slack && lower_slack < half_bin_width)
160         {
161           upper_limit += 2;
162           upper_slack += 2 * half_bin_width;
163         }
164
165       if (upper_slack > lower_slack)
166         {
167           testing_assert (upper_slack > half_bin_width);
168
169           /* Adjust the range to the left */
170           lower_limit --;
171           upper_limit --;
172           upper_slack -= half_bin_width;
173           lower_slack += half_bin_width;
174         }
175       else
176         {
177           testing_assert (lower_slack >= half_bin_width);
178
179           /* Adjust the range to the right */
180           lower_limit ++;
181           upper_limit ++;
182           lower_slack -= half_bin_width;
183           upper_slack += half_bin_width;
184         }
185     }
186
187   /* If there are any completely empty bins, then remove them,
188      since empty bins don't really add much information to the histogram.
189    */
190   if (upper_slack > 2 * half_bin_width)
191     {
192       upper_slack -= 2 * half_bin_width;
193       upper_limit -=2;
194     }
195
196   if (lower_slack >= 2 * half_bin_width)
197     {
198       lower_slack -= 2 * half_bin_width;
199       lower_limit +=2;
200     }
201
202   *adj_min = lower_limit * half_bin_width;
203   *adj_max = upper_limit * half_bin_width;
204
205   testing_assert (*adj_max > max);
206   testing_assert (*adj_min <= min);
207
208   return (upper_limit - lower_limit) / 2.0;
209 }
210
211
212
213 struct histogram *
214 histogram_create (double bin_width, double min, double max)
215 {
216   const int MAX_BINS = 25;
217   struct histogram *h;
218   struct statistic *stat;
219   int bins;
220   double adjusted_min, adjusted_max;
221
222   if (max == min)
223     {
224       msg (MW, _("Not creating histogram because the data contains less than 2 distinct values"));
225       return NULL;
226     }
227
228   assert (bin_width > 0);
229
230   bin_width = chart_rounded_tick (bin_width);
231   bins = adjust_bin_ranges (bin_width, min, max, &adjusted_min, &adjusted_max);
232
233   /* Force the number of bins to lie in a sensible range. */
234   if (bins > MAX_BINS) 
235     {
236       bin_width = chart_rounded_tick ((max - min) / (double) (MAX_BINS - 1));
237       bins = adjust_bin_ranges (bin_width,
238                                 min, max, &adjusted_min, &adjusted_max);
239     }
240
241   /* Can this ever happen? */
242   if (bins < 1)
243     bins = 1;
244
245   h = xmalloc (sizeof *h);
246
247   h->gsl_hist = gsl_histogram_alloc (bins);
248
249   gsl_histogram_set_ranges_uniform (h->gsl_hist, adjusted_min, adjusted_max);
250
251   stat = &h->parent;
252   stat->accumulate = acc;
253   stat->destroy = destroy;
254
255   return h;
256 }
257