Fixed the categoricals such that now both GLM and ONEWAY work
[pspp] / src / math / categoricals.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2009, 2010, 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/categoricals.h"
20
21 #include <stdio.h>
22
23 #include "data/case.h"
24 #include "data/value.h"
25 #include "data/variable.h"
26 #include "libpspp/array.h"
27 #include "libpspp/hmap.h"
28 #include "libpspp/pool.h"
29 #include "libpspp/str.h"
30
31 #include "gl/xalloc.h"
32
33 struct value_node
34 {
35   struct hmap_node node;      /* Node in hash map. */
36   union value value;          /* The value being labeled. */
37   double cc;                  /* The total of the weights of cases with this value */
38
39   void *user_data;            /* A pointer to data which the caller can store stuff */
40
41   int subscript;              /* A zero based integer, unique within the variable.
42                                  Can be used as an index into an array */
43 };
44
45 struct var_params
46 {
47   /* A map indexed by a union values */
48   struct hmap map;
49
50   const struct variable *var;
51
52   int base_subscript_short;
53   int base_subscript_long;
54
55   /* The number of distinct values of this variable */
56   int n_cats;
57
58   /* A map of values indexed by subscript */
59   const struct value_node **reverse_value_map;
60
61   /* Total of the weights of this variable */
62   double cc; 
63 };
64
65
66 /* Comparison function to sort the reverse_value_map in ascending order */
67 static int
68 compare_value_node (const void *vn1_, const void *vn2_, const void *aux)
69 {
70   const struct value_node * const *vn1 = vn1_;
71   const struct value_node * const *vn2 = vn2_;
72   const struct var_params *vp = aux;
73
74   return value_compare_3way (&(*vn1)->value, &(*vn2)->value, var_get_width (vp->var));
75 }
76
77
78 struct categoricals
79 {
80   /* The weight variable */
81   const struct variable *wv;
82
83
84   /* An array of var_params */
85   struct var_params *vp;
86
87   /* The size of VP. (ie, the number of variables involved.) */
88   size_t n_vp;
89
90   /* The number of categorical variables which contain entries.
91      In the absence of missing values, this will be equal to N_VP */
92   size_t n_vars;
93
94   /* A map to enable the lookup of variables indexed by subscript.
95      This map considers only the N - 1 of the N variables.
96    */
97   int *reverse_variable_map_short;
98
99   /* Like the above, but uses all N variables */
100   int *reverse_variable_map_long;
101
102   size_t n_cats_total;
103
104   struct pool *pool;
105
106   /* Missing values to be excluded */
107   enum mv_class exclude;
108
109   /* Function to be called on each update */
110   update_func *update;
111
112   /* Function specified by the caller to create user_data */
113   user_data_create_func *user_data_create;
114
115   /* Auxilliary data to be passed to update and user_data_create_func*/
116   void *aux1;
117   void *aux2;
118 };
119
120
121 void
122 categoricals_destroy ( struct categoricals *cat)
123 {
124   int i;
125   if (cat != NULL)
126     {
127       for (i = 0 ; i < cat->n_vp; ++i)
128         hmap_destroy (&cat->vp[i].map);
129       
130       pool_destroy (cat->pool);
131       free (cat);
132     }
133 }
134
135
136 #if 0
137 void
138 categoricals_dump (const struct categoricals *cat)
139 {
140   int v;
141
142   for (v = 0 ; v < cat->n_vp; ++v)
143     {
144       const struct var_params *vp = &cat->vp[v];
145       const struct hmap *m = &vp->map;
146       struct hmap_node *node ;
147       int x;
148      
149       printf ("\n%s (%d)  CC=%g n_cats=%d:\n", 
150               var_get_name (vp->var), vp->base_subscript_long, vp->cc, vp->n_cats);
151
152       printf ("Reverse map\n");
153       for (x = 0 ; x < vp->n_cats; ++x)
154         {
155           struct string s;
156           const struct value_node *vn = vp->reverse_value_map[x];
157           ds_init_empty (&s);
158           var_append_value_name (vp->var, &vn->value, &s);
159           printf ("Value for %d is %s\n", x, ds_cstr(&s));
160           ds_destroy (&s);
161         }
162
163       printf ("\nForward map\n");
164       for (node = hmap_first (m); node; node = hmap_next (m, node))
165         {
166           struct string s;
167           const struct value_node *vn = HMAP_DATA (node, struct value_node, node);
168           ds_init_empty (&s);
169           var_append_value_name (vp->var, &vn->value, &s);
170           printf ("Value: %s; Index %d; CC %g\n",
171                   ds_cstr (&s),
172                   vn->subscript, vn->cc);
173           ds_destroy (&s);
174         }
175     }
176
177   assert (cat->n_vars <= cat->n_vp);
178
179   printf ("\n");
180   printf ("Number of categorical variables: %d\n", cat->n_vp);
181   printf ("Number of non-empty categorical variables: %d\n", cat->n_vars);
182   printf ("Total number of categories: %d\n", cat->n_cats_total);
183
184   printf ("\nReverse variable map (short):\n");
185   for (v = 0 ; v < cat->n_cats_total - cat->n_vars; ++v)
186     printf ("%d ", cat->reverse_variable_map_short[v]);
187
188   printf ("\nReverse variable map (long):\n");
189   for (v = 0 ; v < cat->n_cats_total; ++v)
190     printf ("%d ", cat->reverse_variable_map_long[v]);
191
192   printf ("\n");
193 }
194 #endif
195
196
197
198 static struct value_node *
199 lookup_value (const struct hmap *map, const struct variable *var, const union value *val)
200 {
201   struct value_node *foo;
202   unsigned int width = var_get_width (var);
203   size_t hash = value_hash (val, width, 0);
204
205   HMAP_FOR_EACH_WITH_HASH (foo, struct value_node, node, hash, map)
206     {
207       if (value_equal (val, &foo->value, width))
208         break;
209
210       fprintf (stderr, "Warning: Hash table collision\n");
211     }
212
213   return foo;
214 }
215
216
217 struct categoricals *
218 categoricals_create (const struct variable *const *v, size_t n_vars,
219                      const struct variable *wv, enum mv_class exclude,
220                      user_data_create_func *udf,
221                      update_func *update, void *aux1, void *aux2
222                      )
223 {
224   size_t i;
225   struct categoricals *cat = xmalloc (sizeof *cat);
226   
227   cat->n_vp = n_vars;
228   cat->wv = wv;
229   cat->n_cats_total = 0;
230   cat->n_vars = 0;
231   cat->reverse_variable_map_short = NULL;
232   cat->reverse_variable_map_long = NULL;
233   cat->pool = pool_create ();
234   cat->exclude = exclude;
235   cat->update = update;
236   cat->user_data_create = udf;
237
238   cat->aux1 = aux1;
239   cat->aux2 = aux2;
240
241
242   cat->vp = pool_calloc (cat->pool, cat->n_vp, sizeof *cat->vp);
243
244   for (i = 0 ; i < cat->n_vp; ++i)
245     {
246       hmap_init (&cat->vp[i].map);
247       cat->vp[i].var = v[i];
248     }
249
250   return cat;
251 }
252
253
254
255 void
256 categoricals_update (struct categoricals *cat, const struct ccase *c)
257 {
258   size_t i;
259   
260   const double weight = cat->wv ? case_data (c, cat->wv)->f : 1.0;
261
262   assert (NULL == cat->reverse_variable_map_short);
263   assert (NULL == cat->reverse_variable_map_long);
264
265   for (i = 0 ; i < cat->n_vp; ++i)
266     {
267       const struct variable *var = cat->vp[i].var;
268       unsigned int width = var_get_width (var);
269       const union value *val = case_data (c, var);
270       size_t hash;
271       struct value_node *node ;
272
273       if ( var_is_value_missing (var, val, cat->exclude))
274         continue;
275
276       hash = value_hash (val, width, 0);
277       node = lookup_value (&cat->vp[i].map, var, val);
278
279       if ( NULL == node)
280         {
281           node = pool_malloc (cat->pool, sizeof *node);
282
283           value_init (&node->value, width);
284           value_copy (&node->value, val, width);
285           node->cc = 0.0;
286
287           hmap_insert (&cat->vp[i].map, &node->node,  hash);
288           cat->n_cats_total++;
289           
290           if ( 0 == cat->vp[i].n_cats)
291             cat->n_vars++;
292
293           node->subscript = cat->vp[i].n_cats++ ;
294
295           if (cat->user_data_create)
296             node->user_data = cat->user_data_create (cat->aux1, cat->aux2);
297         }
298
299       node->cc += weight;
300       cat->vp[i].cc += weight;
301
302       if (cat->update)
303         cat->update (node->user_data, cat->exclude, cat->wv, var, c, cat->aux1, cat->aux2);
304     }
305 }
306
307 /* Return the number of categories (distinct values) for variable N */
308 size_t
309 categoricals_n_count (const struct categoricals *cat, size_t n)
310 {
311   return hmap_count (&cat->vp[n].map);
312 }
313
314
315 /* Return the total number of categories */
316 size_t
317 categoricals_total (const struct categoricals *cat)
318 {
319   return cat->n_cats_total;
320 }
321
322
323 /* This function must be called *before* any call to categoricals_get_*_by subscript and
324  *after* all calls to categoricals_update */
325 void
326 categoricals_done (const struct categoricals *cat_)
327 {
328   /* Implementation Note: Whilst this function is O(n) in cat->n_cats_total, in most
329      uses it will be more efficient that using a tree based structure, since it
330      is called only once, and means that subsequent lookups will be O(1).
331
332      1 call of O(n) + 10^9 calls of O(1) is better than 10^9 calls of O(log n).
333   */
334   struct categoricals *cat = CONST_CAST (struct categoricals *, cat_);
335   int v;
336   int idx_short = 0;
337   int idx_long = 0;
338   cat->reverse_variable_map_short = pool_calloc (cat->pool,
339                                                  cat->n_cats_total - cat->n_vars,
340                                                  sizeof *cat->reverse_variable_map_short);
341
342   cat->reverse_variable_map_long = pool_calloc (cat->pool,
343                                                 cat->n_cats_total,
344                                                 sizeof *cat->reverse_variable_map_long);
345   
346   for (v = 0 ; v < cat->n_vp; ++v)
347     {
348       int i;
349       struct var_params *vp = &cat->vp[v];
350       int n_cats_total = categoricals_n_count (cat, v);
351       struct hmap_node *node ;
352
353       vp->reverse_value_map = pool_calloc (cat->pool, n_cats_total, sizeof *vp->reverse_value_map);
354
355       vp->base_subscript_short = idx_short;
356       vp->base_subscript_long = idx_long;
357
358       for (node = hmap_first (&vp->map); node; node = hmap_next (&vp->map, node))
359         {
360           const struct value_node *vn = HMAP_DATA (node, struct value_node, node);
361           vp->reverse_value_map[vn->subscript] = vn;
362         }
363
364       /* For some purposes (eg CONTRASTS in ONEWAY) the values need to be sorted */
365       sort (vp->reverse_value_map, vp->n_cats, sizeof (const struct value_node *),
366             compare_value_node, vp);
367
368       /* Populate the reverse variable maps. */
369       for (i = 0; i < vp->n_cats - 1; ++i)
370         cat->reverse_variable_map_short[idx_short++] = v;
371
372       for (i = 0; i < vp->n_cats; ++i)
373         cat->reverse_variable_map_long[idx_long++] = v;
374     }
375
376   assert (cat->n_vars <= cat->n_vp);
377 }
378
379
380 static int
381 reverse_variable_lookup_short (const struct categoricals *cat, int subscript)
382 {
383   assert (cat->reverse_variable_map_short);
384   assert (subscript >= 0);
385   assert (subscript < cat->n_cats_total - cat->n_vars);
386
387   return cat->reverse_variable_map_short[subscript];
388 }
389
390 static int
391 reverse_variable_lookup_long (const struct categoricals *cat, int subscript)
392 {
393   assert (cat->reverse_variable_map_long);
394   assert (subscript >= 0);
395   assert (subscript < cat->n_cats_total);
396
397   return cat->reverse_variable_map_long[subscript];
398 }
399
400
401
402 /* Return the categorical variable corresponding to SUBSCRIPT */
403 const struct variable *
404 categoricals_get_variable_by_subscript (const struct categoricals *cat, int subscript)
405 {
406   int index = reverse_variable_lookup_short (cat, subscript);
407
408   return cat->vp[index].var;
409 }
410
411 /* Return the value corresponding to SUBSCRIPT */
412 static const union value *
413 categoricals_get_value_by_subscript (const struct categoricals *cat, int subscript)
414 {
415   int vindex = reverse_variable_lookup_short (cat, subscript);
416   const struct var_params *vp = &cat->vp[vindex];
417   const struct value_node *vn = vp->reverse_value_map [subscript - vp->base_subscript_short];
418
419   return &vn->value;
420 }
421
422
423 double
424 categoricals_get_weight_by_subscript (const struct categoricals *cat, int subscript)
425 {
426   int vindex = reverse_variable_lookup_short (cat, subscript);
427   const struct var_params *vp = &cat->vp[vindex];
428
429   return vp->cc;
430 }
431
432 double
433 categoricals_get_sum_by_subscript (const struct categoricals *cat, int subscript)
434 {
435   int vindex = reverse_variable_lookup_short (cat, subscript);
436   const struct var_params *vp = &cat->vp[vindex];
437
438   const struct value_node *vn = vp->reverse_value_map [subscript - vp->base_subscript_short];
439   return vn->cc;
440 }
441
442
443 /* Returns unity if the value in case C at SUBSCRIPT is equal to the category
444    for that subscript */
445 double
446 categoricals_get_binary_by_subscript (const struct categoricals *cat, int subscript,
447                                       const struct ccase *c)
448 {
449   const struct variable *var = categoricals_get_variable_by_subscript (cat, subscript);
450   int width = var_get_width (var);
451
452   const union value *val = case_data (c, var);
453
454   return value_equal (val, categoricals_get_value_by_subscript (cat, subscript), width);
455 }
456
457
458 size_t
459 categoricals_get_n_variables (const struct categoricals *cat)
460 {
461   return cat->n_vars;
462 }
463
464
465
466 /* Return the value corresponding to SUBSCRIPT */
467 const union value *
468 categoricals_get_value_by_category (const struct categoricals *cat, int subscript)
469 {
470   int vindex = reverse_variable_lookup_long (cat, subscript);
471   const struct var_params *vp = &cat->vp[vindex];
472   const struct value_node *vn = vp->reverse_value_map [subscript - vp->base_subscript_long];
473
474   return &vn->value;
475 }
476
477
478 void *
479 categoricals_get_user_data_by_category (const struct categoricals *cat, int subscript)
480 {
481   int vindex = reverse_variable_lookup_long (cat, subscript);
482   const struct var_params *vp = &cat->vp[vindex];
483
484   const struct value_node *vn = vp->reverse_value_map [subscript - vp->base_subscript_long];
485   return vn->user_data;
486 }