categoricals: Mark dump_interaction() and categoricals_dump() UNUSED.
[pspp] / src / math / categoricals.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2009, 2010, 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/categoricals.h"
20 #include "math/interaction.h"
21
22 #include <stdio.h>
23
24 #include "data/case.h"
25 #include "data/value.h"
26 #include "data/variable.h"
27 #include "libpspp/array.h"
28 #include "libpspp/hmap.h"
29 #include "libpspp/pool.h"
30 #include "libpspp/str.h"
31 #include "libpspp/hash-functions.h"
32
33 #include "gl/xalloc.h"
34
35 #define CATEGORICALS_DEBUG 0
36
37 #define EFFECTS_CODING 1
38
39 struct value_node
40 {
41   struct hmap_node node;      /* Node in hash map. */
42
43   union value val;            /* The value */
44
45   int index;                  /* A zero based unique index for this value */
46 };
47
48 struct interaction_value
49 {
50   struct hmap_node node;      /* Node in hash map */
51
52   struct ccase *ccase;        /* A case (probably the first in the dataset) which matches this value */
53
54   double cc;        /* Total of the weights of cases matching this interaction */
55
56   void *user_data;            /* A pointer to data which the caller can store stuff */
57 };
58
59 static struct value_node *
60 lookup_value (const struct hmap *map, const union value *val, unsigned int hash, int width)
61 {
62   struct value_node *vn = NULL;
63   HMAP_FOR_EACH_WITH_HASH (vn, struct value_node, node, hash, map)
64     {
65       if (value_equal (&vn->val, val, width))
66         break;
67     }
68   
69   return vn;
70 }
71
72 struct variable_node
73 {
74   struct hmap_node node;      /* Node in hash map. */
75   const struct variable *var; /* The variable */
76
77   struct hmap valmap;         /* A map of value nodes */
78   int n_vals;                 /* Number of values for this variable */
79 };
80
81 static void UNUSED
82 dump_interaction (const struct interaction *iact)
83 {
84   if (CATEGORICALS_DEBUG)
85     {
86       struct string str = DS_EMPTY_INITIALIZER;
87       interaction_to_string (iact, &str);
88       printf ("Interaction: %s\n", ds_cstr (&str));
89       ds_destroy (&str);
90     }
91 }
92
93 static struct variable_node *
94 lookup_variable (const struct hmap *map, const struct variable *var, unsigned int hash)
95 {
96   struct variable_node *vn = NULL;
97   HMAP_FOR_EACH_WITH_HASH (vn, struct variable_node, node, hash, map)
98     {
99       if (vn->var == var)
100         break;
101       
102       fprintf (stderr, "Warning: Hash table collision\n");
103     }
104   
105   return vn;
106 }
107
108
109 struct interact_params
110 {
111   /* A map indexed by a interaction_value */
112   struct hmap ivmap;
113
114   const struct interaction *iact;
115
116   int base_subscript_short;
117   int base_subscript_long;
118
119   /* The number of distinct values of this interaction */
120   int n_cats;
121
122   /* An array of integers df_n * df_{n-1} * df_{n-2} ...
123      These are the products of the degrees of freedom for the current 
124      variable and all preceeding variables */
125   int *df_prod; 
126
127   double *enc_sum;
128
129   /* A map of interaction_values indexed by subscript */
130   struct interaction_value **reverse_interaction_value_map;
131
132   double cc;
133 };
134
135
136 /* Comparison function to sort the reverse_value_map in ascending order */
137 static int
138 compare_interaction_value_3way (const void *vn1_, const void *vn2_, const void *aux)
139 {
140   const struct interaction_value *const *vn1p = vn1_;
141   const struct interaction_value *const *vn2p = vn2_;
142
143   const struct interact_params *iap = aux;
144
145   return interaction_case_cmp_3way (iap->iact, (*vn1p)->ccase, (*vn2p)->ccase);
146 }
147
148 struct categoricals
149 {
150   /* The weight variable */
151   const struct variable *wv;
152
153   /* An array of interact_params */
154   struct interact_params *iap;
155
156   /* Map whose members are the union of the variables which comprise IAP */
157   struct hmap varmap;
158
159   /* The size of IAP. (ie, the number of interactions involved.) */
160   size_t n_iap;
161
162   /* The number of categorical variables which contain entries.
163      In the absence of missing values, this will be equal to N_IAP */
164   size_t n_vars;
165
166   size_t df_sum;
167
168   /* A map to enable the lookup of variables indexed by subscript.
169      This map considers only the N - 1 of the N variables.
170   */
171   int *reverse_variable_map_short;
172
173   /* Like the above, but uses all N variables */
174   int *reverse_variable_map_long;
175
176   size_t n_cats_total;
177
178   struct pool *pool;
179
180   /* Missing values to be excluded */
181   enum mv_class exclude;
182
183   const void *aux1;
184   void *aux2;
185
186   const struct payload *payload;
187 };
188
189 static void UNUSED
190 categoricals_dump (const struct categoricals *cat)
191 {
192   if (CATEGORICALS_DEBUG)
193     {
194       int i;
195
196       printf ("Reverse Variable Map (short):\n");
197       for (i = 0; i < cat->df_sum; ++i)
198         {
199           printf (" %d", cat->reverse_variable_map_short[i]);
200         }
201       printf ("\n");
202
203       printf ("Reverse Variable Map (long):\n");
204       for (i = 0; i < cat->n_cats_total; ++i)
205         {
206           printf (" %d", cat->reverse_variable_map_long[i]);
207         }
208       printf ("\n");
209
210
211       printf ("Number of interactions %d\n", cat->n_iap);
212       for (i = 0 ; i < cat->n_iap; ++i)
213         {
214           int v;
215           struct string str;
216           const struct interact_params *iap = &cat->iap[i];
217           const struct interaction *iact = iap->iact;
218
219           ds_init_empty (&str);
220           interaction_to_string (iact, &str);
221
222           printf ("\nInteraction: %s (n: %d); ", ds_cstr (&str), iap->n_cats);
223           ds_destroy (&str);
224           printf ("Base subscript: %d\n", iap->base_subscript_short);
225
226           printf ("\t(");
227           for (v = 0; v < hmap_count (&iap->ivmap); ++v)
228             {
229               int vv;
230               const struct interaction_value *iv = iap->reverse_interaction_value_map[v];
231           
232               if (v > 0)  printf ("   ");
233               printf ("{");
234               for (vv = 0; vv < iact->n_vars; ++vv)
235                 {
236                   const struct variable *var = iact->vars[vv];
237                   const union value *val = case_data (iv->ccase, var);
238               
239                   printf ("%g", val->f);
240                   if (vv < iact->n_vars - 1)
241                     printf (", ");
242                 }
243               printf ("}");
244             }
245           printf (")\n");
246         }
247     }
248 }
249
250 void
251 categoricals_destroy (struct categoricals *cat)
252 {
253   struct variable_node *vn = NULL;
254   int i;
255   if (NULL == cat)
256     return;
257   for (i = 0; i < cat->n_iap; ++i)
258     {
259       struct interaction_value *iv = NULL;
260       /* Interate over each interaction value, and unref any cases that we reffed */
261       HMAP_FOR_EACH (iv, struct interaction_value, node, &cat->iap[i].ivmap)
262         {
263 #if 0
264           if (cat->payload)
265             cat->payload->destroy (cat->aux1, iv->user_data);
266 #endif
267           case_unref (iv->ccase);
268         }
269
270       free (cat->iap[i].enc_sum);
271       free (cat->iap[i].df_prod);
272       hmap_destroy (&cat->iap[i].ivmap);
273     }
274
275   /* Interate over each variable and delete its value map */
276   HMAP_FOR_EACH (vn, struct variable_node, node, &cat->varmap)
277     {
278       hmap_destroy (&vn->valmap);
279     }
280
281   hmap_destroy (&cat->varmap);
282
283   pool_destroy (cat->pool);
284
285   free (cat);
286 }
287
288
289
290 static struct interaction_value *
291 lookup_case (const struct hmap *map, const struct interaction *iact, const struct ccase *c)
292 {
293   struct interaction_value *iv = NULL;
294   size_t hash = interaction_case_hash (iact, c, 0);
295
296   HMAP_FOR_EACH_WITH_HASH (iv, struct interaction_value, node, hash, map)
297     {
298       if (interaction_case_equal (iact, c, iv->ccase))
299         break;
300
301       fprintf (stderr, "Warning: Hash table collision\n");
302     }
303
304   return iv;
305 }
306
307
308 struct categoricals *
309 categoricals_create (struct interaction *const*inter, size_t n_inter,
310                      const struct variable *wv, enum mv_class exclude)
311 {
312   size_t i;
313   struct categoricals *cat = xmalloc (sizeof *cat);
314   
315   cat->n_iap = n_inter;
316   cat->wv = wv;
317   cat->n_cats_total = 0;
318   cat->n_vars = 0;
319   cat->reverse_variable_map_short = NULL;
320   cat->reverse_variable_map_long = NULL;
321   cat->pool = pool_create ();
322   cat->exclude = exclude;
323   cat->payload = NULL;
324   cat->aux2 = NULL;
325
326   cat->iap = pool_calloc (cat->pool, cat->n_iap, sizeof *cat->iap);
327
328   hmap_init (&cat->varmap);
329   for (i = 0 ; i < cat->n_iap; ++i)
330     {
331       int v;
332       hmap_init (&cat->iap[i].ivmap);
333       cat->iap[i].iact = inter[i];
334       cat->iap[i].cc = 0.0;
335       for (v = 0; v < inter[i]->n_vars; ++v)
336         {
337           const struct variable *var = inter[i]->vars[v];
338           unsigned int hash = hash_pointer (var, 0);
339           struct variable_node *vn = lookup_variable (&cat->varmap, var, hash);
340           if (vn == NULL)
341             {
342               vn = pool_malloc (cat->pool, sizeof *vn);
343               vn->var = var;
344               vn->n_vals = 0;
345               hmap_init (&vn->valmap);
346
347               hmap_insert (&cat->varmap, &vn->node,  hash);
348             }
349         }
350     }
351
352   return cat;
353 }
354
355
356
357 void
358 categoricals_update (struct categoricals *cat, const struct ccase *c)
359 {
360   int i;
361   struct variable_node *vn = NULL;
362   const double weight = cat->wv ? case_data (c, cat->wv)->f : 1.0;
363
364   assert (NULL == cat->reverse_variable_map_short);
365   assert (NULL == cat->reverse_variable_map_long);
366
367   /* Interate over each variable, and add the value of that variable
368      to the appropriate map, if it's not already present. */
369   HMAP_FOR_EACH (vn, struct variable_node, node, &cat->varmap)
370     {
371       const int width = var_get_width (vn->var);
372       const union value *val = case_data (c, vn->var);
373       unsigned int hash = value_hash (val, width, 0);
374
375       struct value_node *valn = lookup_value (&vn->valmap, val, hash, width);
376       if (valn == NULL)
377         {
378           valn = pool_malloc (cat->pool, sizeof *valn);
379           valn->index = vn->n_vals++;
380           value_init (&valn->val, width);
381           value_copy (&valn->val, val, width);
382           hmap_insert (&vn->valmap, &valn->node, hash);
383         }
384     }     
385   
386   for (i = 0 ; i < cat->n_iap; ++i)
387     {
388       const struct interaction *iact = cat->iap[i].iact;
389
390       size_t hash;
391       struct interaction_value *node;
392
393       if ( interaction_case_is_missing (iact, c, cat->exclude))
394         continue;
395
396       hash = interaction_case_hash (iact, c, 0);
397       node = lookup_case (&cat->iap[i].ivmap, iact, c);
398
399       if ( NULL == node)
400         {
401           node = pool_malloc (cat->pool, sizeof *node);
402
403           node->ccase = case_ref (c);
404           node->cc = weight;
405
406           hmap_insert (&cat->iap[i].ivmap, &node->node, hash);
407
408           if (cat->payload) 
409             {
410               node->user_data = cat->payload->create (cat->aux1, cat->aux2);
411             }
412         }
413       else
414         {
415           node->cc += weight;
416         }
417       cat->iap[i].cc += weight;
418
419       if (cat->payload)
420         {
421           double weight = cat->wv ? case_data (c, cat->wv)->f : 1.0;
422           cat->payload->update (cat->aux1, cat->aux2, node->user_data, c, weight);
423         }
424
425     }
426 }
427
428 /* Return the number of categories (distinct values) for interction N */
429 size_t
430 categoricals_n_count (const struct categoricals *cat, size_t n)
431 {
432   return hmap_count (&cat->iap[n].ivmap);
433 }
434
435
436 size_t
437 categoricals_df (const struct categoricals *cat, size_t n)
438 {
439   const struct interact_params *iap = &cat->iap[n];
440   return iap->df_prod[iap->iact->n_vars - 1];
441 }
442
443
444 /* Return the total number of categories */
445 size_t
446 categoricals_n_total (const struct categoricals *cat)
447 {
448   if (!categoricals_is_complete (cat))
449     return 0;
450
451   return cat->n_cats_total;
452 }
453
454 size_t
455 categoricals_df_total (const struct categoricals *cat)
456 {
457   return cat->df_sum;
458 }
459
460 bool
461 categoricals_is_complete (const struct categoricals *cat)
462 {
463   return (NULL != cat->reverse_variable_map_short);
464 }
465
466
467 /* This function must be called *before* any call to categoricals_get_*_by subscript and
468  *after* all calls to categoricals_update */
469 bool
470 categoricals_done (const struct categoricals *cat_)
471 {
472   /* Implementation Note: Whilst this function is O(n) in cat->n_cats_total, in most
473      uses it will be more efficient that using a tree based structure, since it
474      is called only once, and means that subsequent lookups will be O(1).
475
476      1 call of O(n) + 10^9 calls of O(1) is better than 10^9 calls of O(log n).
477   */
478   struct categoricals *cat = CONST_CAST (struct categoricals *, cat_);
479   int v;
480   int i;
481   int idx_short = 0;
482   int idx_long = 0;
483   cat->df_sum = 0;
484   cat->n_cats_total = 0;
485
486   /* Calculate the degrees of freedom, and the number of categories */
487   for (i = 0 ; i < cat->n_iap; ++i)
488     {
489       int df = 1;
490       const struct interaction *iact = cat->iap[i].iact;
491
492       cat->iap[i].df_prod = iact->n_vars ? xcalloc (iact->n_vars, sizeof (int)) : NULL;
493
494       cat->iap[i].n_cats = 1;
495       
496       for (v = 0 ; v < iact->n_vars; ++v)
497         {
498           const struct variable *var = iact->vars[v];
499
500           struct variable_node *vn = lookup_variable (&cat->varmap, var, hash_pointer (var, 0));
501
502           if  (hmap_count (&vn->valmap) == 0)
503             return false;
504
505           cat->iap[i].df_prod[v] = df * (hmap_count (&vn->valmap) - 1);
506           df = cat->iap[i].df_prod[v];
507
508           cat->iap[i].n_cats *= hmap_count (&vn->valmap);
509         }
510
511       assert (v == iact->n_vars);
512       if (v > 0)
513         cat->df_sum += cat->iap[i].df_prod [v - 1];
514
515       cat->n_cats_total += cat->iap[i].n_cats;
516     }
517
518
519   cat->reverse_variable_map_short = pool_calloc (cat->pool,
520                                                  cat->df_sum,
521                                                  sizeof *cat->reverse_variable_map_short);
522
523   cat->reverse_variable_map_long = pool_calloc (cat->pool,
524                                                 cat->n_cats_total,
525                                                 sizeof *cat->reverse_variable_map_long);
526
527   for (i = 0 ; i < cat->n_iap; ++i)
528     {
529       struct interaction_value *ivn = NULL;
530       int x = 0;
531       int ii;
532       struct interact_params *iap = &cat->iap[i];
533
534       iap->base_subscript_short = idx_short;
535       iap->base_subscript_long = idx_long;
536
537       iap->reverse_interaction_value_map = pool_calloc (cat->pool, iap->n_cats,
538                                                         sizeof *iap->reverse_interaction_value_map);
539
540       HMAP_FOR_EACH (ivn, struct interaction_value, node, &iap->ivmap)
541         {
542           iap->reverse_interaction_value_map[x++] = ivn;
543         }
544
545       assert (x <= iap->n_cats);
546
547       /* For some purposes (eg CONTRASTS in ONEWAY) the values need to be sorted */
548       sort (iap->reverse_interaction_value_map, x, sizeof (*iap->reverse_interaction_value_map),
549             compare_interaction_value_3way, iap);
550
551       /* Fill the remaining values with null */
552       for (ii = x ; ii < iap->n_cats; ++ii)
553         iap->reverse_interaction_value_map[ii] = NULL;
554
555       /* Populate the reverse variable maps. */
556       if (iap->df_prod)
557         {
558           for (ii = 0; ii < iap->df_prod [iap->iact->n_vars - 1]; ++ii)
559             cat->reverse_variable_map_short[idx_short++] = i;
560         }
561
562       for (ii = 0; ii < iap->n_cats; ++ii)
563         cat->reverse_variable_map_long[idx_long++] = i;
564     }
565
566   assert (cat->n_vars <= cat->n_iap);
567
568   //  categoricals_dump (cat);
569
570   /* Tally up the sums for all the encodings */
571   for (i = 0 ; i < cat->n_iap; ++i)
572     {
573       int x, y;
574       struct interact_params *iap = &cat->iap[i];
575       const struct interaction *iact = iap->iact;
576
577       const int df = iap->df_prod ? iap->df_prod [iact->n_vars - 1] : 0;
578
579       iap->enc_sum = xcalloc (df, sizeof (*(iap->enc_sum)));
580
581       for (y = 0; y < hmap_count (&iap->ivmap); ++y)
582         {
583           struct interaction_value *iv = iap->reverse_interaction_value_map[y];
584           for (x = iap->base_subscript_short; x < iap->base_subscript_short + df ;++x)
585             {
586               const double bin = categoricals_get_code_for_case (cat, x, iv->ccase); \
587               iap->enc_sum [x - iap->base_subscript_short] += bin * iv->cc;
588             }
589         }
590     }
591
592   return true;
593 }
594
595
596 static int
597 reverse_variable_lookup_short (const struct categoricals *cat, int subscript)
598 {
599   assert (cat->reverse_variable_map_short);
600   assert (subscript >= 0);
601   assert (subscript < cat->df_sum);
602
603   return cat->reverse_variable_map_short[subscript];
604 }
605
606 static int
607 reverse_variable_lookup_long (const struct categoricals *cat, int subscript)
608 {
609   assert (cat->reverse_variable_map_long);
610   assert (subscript >= 0);
611   assert (subscript < cat->n_cats_total);
612
613   return cat->reverse_variable_map_long[subscript];
614 }
615
616
617 /* Return the interaction corresponding to SUBSCRIPT */
618 const struct interaction *
619 categoricals_get_interaction_by_subscript (const struct categoricals *cat, int subscript)
620 {
621   int index = reverse_variable_lookup_short (cat, subscript);
622
623   return cat->iap[index].iact;
624 }
625
626 double
627 categoricals_get_weight_by_subscript (const struct categoricals *cat, int subscript)
628 {
629   int vindex = reverse_variable_lookup_short (cat, subscript);
630   const struct interact_params *vp = &cat->iap[vindex];
631
632   return vp->cc;
633 }
634
635 double
636 categoricals_get_sum_by_subscript (const struct categoricals *cat, int subscript)
637 {
638   int vindex = reverse_variable_lookup_short (cat, subscript);
639   const struct interact_params *vp = &cat->iap[vindex];
640
641   return   vp->enc_sum[subscript - vp->base_subscript_short];
642 }
643
644 /* Returns unity if the value in case C at SUBSCRIPT is equal to the category
645    for that subscript */
646 double
647 categoricals_get_code_for_case (const struct categoricals *cat, int subscript,
648                                 const struct ccase *c)
649 {
650   const struct interaction *iact = categoricals_get_interaction_by_subscript (cat, subscript);
651
652   const int i = reverse_variable_lookup_short (cat, subscript);
653
654   const int base_index = cat->iap[i].base_subscript_short;
655
656   int v;
657   double result = 1.0;
658
659   const struct interact_params *iap = &cat->iap[i];
660
661   double dfp = 1.0;
662   for (v = 0; v < iact->n_vars; ++v)
663     {
664       const struct variable *var = iact->vars[v];
665
666       const union value *val = case_data (c, var);
667       const int width = var_get_width (var);
668       const struct variable_node *vn = lookup_variable (&cat->varmap, var, hash_pointer (var, 0));
669
670       const unsigned int hash = value_hash (val, width, 0);
671       const struct value_node *valn = lookup_value (&vn->valmap, val, hash, width);
672
673       double bin = 1.0;
674
675       const double df = iap->df_prod[v] / dfp;
676
677       /* Translate the subscript into an index for the individual variable */
678       const int index = ((subscript - base_index) % iap->df_prod[v] ) / dfp;
679       dfp = iap->df_prod [v];
680
681 #if EFFECTS_CODING
682       if ( valn->index == df )
683         bin = -1.0;
684       else 
685 #endif
686         if ( valn->index  != index )
687           bin = 0;
688     
689       result *= bin;
690     }
691
692   return result;
693 }
694
695
696 size_t
697 categoricals_get_n_variables (const struct categoricals *cat)
698 {
699   printf ("%s\n", __FUNCTION__);
700   return cat->n_vars;
701 }
702
703
704 /* Return a case containing the set of values corresponding to 
705    the Nth Category of the IACTth interaction */
706 const struct ccase *
707 categoricals_get_case_by_category_real (const struct categoricals *cat, int iact, int n)
708 {
709   const struct interact_params *vp = &cat->iap[iact];
710   const struct interaction_value *vn = vp->reverse_interaction_value_map [n];
711
712   return vn->ccase;
713 }
714
715 /* Return a the user data corresponding to the Nth Category of the IACTth interaction. */
716 void *
717 categoricals_get_user_data_by_category_real (const struct categoricals *cat, int iact, int n)
718 {
719   const struct interact_params *vp = &cat->iap[iact];
720   const struct interaction_value *iv = vp->reverse_interaction_value_map [n];
721
722   return iv->user_data;
723 }
724
725
726
727 /* Return a case containing the set of values corresponding to SUBSCRIPT */
728 const struct ccase *
729 categoricals_get_case_by_category (const struct categoricals *cat, int subscript)
730 {
731   int vindex = reverse_variable_lookup_long (cat, subscript);
732   const struct interact_params *vp = &cat->iap[vindex];
733   const struct interaction_value *vn = vp->reverse_interaction_value_map [subscript - vp->base_subscript_long];
734
735   return vn->ccase;
736 }
737
738 void *
739 categoricals_get_user_data_by_category (const struct categoricals *cat, int subscript)
740 {
741   int vindex = reverse_variable_lookup_long (cat, subscript);
742   const struct interact_params *vp = &cat->iap[vindex];
743
744   const struct interaction_value *iv = vp->reverse_interaction_value_map [subscript - vp->base_subscript_long];
745   return iv->user_data;
746 }
747
748
749 \f
750
751 void
752 categoricals_set_payload (struct categoricals *cat, const struct payload *p,
753                           const void *aux1, void *aux2)
754 {
755   cat->payload = p;
756   cat->aux1 = aux1;
757   cat->aux2 = aux2;
758 }