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