gui: Allow File|Open to select an encoding for system files.
[pspp] / src / math / percentiles.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2008, 2009, 2011 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/percentiles.h"
20
21 #include "data/casereader.h"
22 #include "data/val-type.h"
23 #include "data/variable.h"
24 #include "libpspp/assertion.h"
25 #include "libpspp/cast.h"
26 #include "math/order-stats.h"
27
28 #include "gl/xalloc.h"
29
30 #include "gettext.h"
31 #define _(msgid) gettext (msgid)
32 #define N_(msgid) msgid
33
34 const char *const ptile_alg_desc[] = {
35   "",
36   N_("HAverage"),
37   N_("Weighted Average"),
38   N_("Rounded"),
39   N_("Empirical"),
40   N_("Empirical with averaging")
41 };
42
43
44
45 double
46 percentile_calculate (const struct percentile *ptl, enum pc_alg alg)
47 {
48   struct percentile *mutable = CONST_CAST (struct percentile *, ptl);
49   const struct order_stats *os = &ptl->parent;
50
51   if ( ptl->g1 == SYSMIS)
52     mutable->g1 = (os->k[0].tc - os->k[0].cc) / os->k[0].c_p1;
53
54   if ( ptl->g1_star == SYSMIS)
55     mutable->g1_star = os->k[0].tc - os->k[0].cc;
56
57   if ( ptl->g2 == SYSMIS)
58     {
59       if ( os->k[1].c == 0 )
60         mutable->g2 = os->k[1].tc / os->k[1].c_p1;
61       else if ( os->k[1].c_p1 == 0 )
62         mutable->g2 = 0;
63       else
64         mutable->g2 = (os->k[1].tc - os->k[1].cc) / os->k[1].c_p1;
65     }
66
67   if ( ptl->g2_star == SYSMIS)
68     {
69       if ( os->k[1].c == 0 )
70         mutable->g2_star = os->k[1].tc;
71       else if ( os->k[1].c_p1 == 0 )
72         mutable->g2_star = 0;
73       else
74         mutable->g2_star = os->k[1].tc - os->k[1].cc;
75     }
76
77   switch (alg)
78     {
79     case PC_WAVERAGE:
80       if ( ptl->g1_star >= 1.0)
81         return os->k[0].y_p1;
82       else
83         {
84           double a = ( os->k[0].y == SYSMIS ) ? 0 : os->k[0].y;
85
86           if (os->k[0].c_p1 >= 1.0)
87             return (1 - ptl->g1_star) * a + ptl->g1_star * os->k[0].y_p1;
88           else
89             return (1 - ptl->g1) * a + ptl->g1 * os->k[0].y_p1;
90         }
91       break;
92
93     case PC_ROUND:
94       {
95         double a = ( os->k[0].y == SYSMIS ) ? 0 : os->k[0].y;
96
97         if (os->k[0].c_p1 >= 1.0)
98           return (ptl->g1_star < 0.5) ? a : os->k[0].y_p1;
99         else
100           return (ptl->g1 < 0.5) ? a : os->k[0].y_p1;
101       }
102       break;
103
104     case PC_EMPIRICAL:
105       if ( ptl->g1_star == 0 )
106         return os->k[0].y;
107       else
108         return os->k[0].y_p1;
109       break;
110
111     case PC_HAVERAGE:
112       if ( ptl->g2_star >= 1.0)
113         {
114           return os->k[1].y_p1;
115         }
116       else
117         {
118           double a = ( os->k[1].y == SYSMIS ) ? 0 : os->k[1].y;
119
120           if ( os->k[1].c_p1 >= 1.0)
121             {
122               if ( ptl->g2_star == 0)
123                 return os->k[1].y;
124
125               return (1 - ptl->g2_star) * a + ptl->g2_star * os->k[1].y_p1;
126             }
127           else
128             {
129               return (1 - ptl->g2) * a + ptl->g2 * os->k[1].y_p1;
130             }
131         }
132
133       break;
134
135     case PC_AEMPIRICAL:
136       if ( ptl->g1_star == 0 )
137         return (os->k[0].y + os->k[0].y_p1)/ 2.0;
138       else
139         return os->k[0].y_p1;
140       break;
141
142     default:
143       NOT_REACHED ();
144       break;
145     }
146
147   NOT_REACHED ();
148
149   return SYSMIS;
150 }
151
152
153 static void
154 destroy (struct statistic *stat)
155 {
156   struct percentile *ptl = UP_CAST (stat, struct percentile, parent.parent);
157   struct order_stats *os = &ptl->parent;
158   free (os->k);
159   free (ptl);
160 }
161
162
163 struct percentile *
164 percentile_create (double p, double W)
165 {
166   struct percentile *ptl = xzalloc (sizeof (*ptl));
167   struct order_stats *os = &ptl->parent;
168   struct statistic *stat = &os->parent;
169
170   assert (p >= 0);
171   assert (p <= 1.0);
172
173   ptl->ptile = p;
174   ptl->w = W;
175
176   os->n_k = 2;
177   os->k = xcalloc (2, sizeof (*os->k));
178   os->k[0].tc = W * p;
179   os->k[1].tc = (W + 1.0) * p;
180
181   ptl->g1 = ptl->g1_star = SYSMIS;
182   ptl->g2 = ptl->g2_star = SYSMIS;
183
184   os->k[1].y_p1 = os->k[1].y = SYSMIS;
185   os->k[0].y_p1 = os->k[0].y = SYSMIS;
186
187   stat->destroy = destroy;
188
189   return ptl;
190 }
191