src/cat.c

   1 /* PSPP - linear regression.
   2    Copyright (C) 2005 Free Software Foundation, Inc.
   3    Written by Jason H Stover <jason@sakla.net>.
   4
   5    This program is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU General Public License as
   7    published by the Free Software Foundation; either version 2 of the
   8    License, or (at your option) any later version.
   9
  10    This program is distributed in the hope that it will be useful, but
  11    WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    General Public License for more details.
  14
  15    You should have received a copy of the GNU General Public License
  16    along with this program; if not, write to the Free Software
  17    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  18    02110-1301, USA. */
  19
  20 /*
  21   Functions and data structures to recode categorical variables into
  22   vectors and sub-rows of matrices.
  23
  24   For some statistical models, it is necessary to change each value
  25   of a categorical variable to a vector with binary entries. These
  26   vectors are then stored as sub-rows within a matrix during
  27   model-fitting. E.g., we need functions and data strucutres to map a
  28   value, say 'a', of a variable named 'cat_var', to a vector, say (0
  29   1 0 0 0), and vice versa.  We also need to be able to map the
  30   vector back to the value 'a', and if the vector is a sub-row of a
  31   matrix, we need to know which sub-row corresponds to the variable
  32   'cat_var'.
  33
  34   The data structures defined here will be placed in the variable
  35   structure in the future. When that happens, the useful code
  36   in this file will be that which refers to design matrices.
  37 */
  38 #include <config.h>
  39 #include <stdlib.h>
  40 #include <error.h>
  41 #include "alloc.h"
  42 #include "error.h"
  43 #include "var.h"
  44 #include "cat.h"
  45 #include <string.h>
  46 #include <math.h>
  47 #include <gsl/gsl_machine.h>
  48 #include <gsl/gsl_vector.h>
  49 #include <gsl/gsl_matrix.h>
  50
  51 #define N_INITIAL_CATEGORIES 1
  52 #define CR_COLUMN_NOT_FOUND -1
  53 #define CR_VALUE_NOT_FOUND -2
  54 #define CR_INDEX_NOT_FOUND -3
  55
  56 struct recoded_categorical *
  57 cr_recoded_categorical_create (const struct variable *v)
  58 {
  59   struct recoded_categorical *rc;
  60
  61   rc = xmalloc (sizeof (*rc));
  62   rc->v = v;
  63   rc->n_categories = 0;
  64   rc->n_allocated_categories = N_INITIAL_CATEGORIES;
  65   rc->vals = (union value **) xmalloc (N_INITIAL_CATEGORIES *
  66                                        sizeof (*rc->vals));
  67
  68   return rc;
  69 }
  70
  71 void
  72 cr_recoded_categorical_destroy (struct recoded_categorical *r)
  73 {
  74   free (r->vals);
  75   free (r);
  76 }
  77
  78 struct recoded_categorical_array *
  79 cr_recoded_cat_ar_create (int n_variables, struct variable *v_variables[])
  80 {
  81   size_t n_categoricals = 0;
  82   size_t i;
  83   struct recoded_categorical_array *ca;
  84   struct variable *v;
  85
  86   ca = (struct recoded_categorical_array *) xmalloc (sizeof (*ca));
  87   for (i = 0; i < n_variables; i++)
  88     {
  89       v = v_variables[i];
  90       if (v->type == ALPHA)
  91         {
  92           n_categoricals++;
  93         }
  94     }
  95   ca->n_vars = n_categoricals;
  96   ca->a = xmalloc (n_categoricals * sizeof (*(ca->a)));
  97   for (i = 0; i < n_categoricals; i++)
  98     {
  99       *(ca->a + i) = cr_recoded_categorical_create (v_variables[i]);
 100     }
 101
 102   return ca;
 103 }
 104
 105 int
 106 cr_free_recoded_array (struct recoded_categorical_array *r)
 107 {
 108   int rc = 0;
 109   size_t i;
 110
 111   for (i = 0; i < r->n_vars; i++)
 112     {
 113       cr_recoded_categorical_destroy (*(r->a + i));
 114     }
 115   return rc;
 116 }
 117
 118 static size_t
 119 cr_value_find (struct recoded_categorical *rc, const union value *v)
 120 {
 121   size_t i;
 122   const union value *val;
 123
 124   for (i = 0; i < rc->n_categories; i++)
 125     {
 126       val = *(rc->vals + i);
 127       if (!compare_values (val, v, rc->v->width))
 128         {
 129           return i;
 130         }
 131     }
 132   return CR_VALUE_NOT_FOUND;
 133 }
 134
 135 /*
 136    Add the new value unless it is already present.
 137  */
 138 void
 139 cr_value_update (struct recoded_categorical *rc, const union value *v)
 140 {
 141   if (cr_value_find (rc, v) == CR_VALUE_NOT_FOUND)
 142     {
 143       if (rc->n_categories >= rc->n_allocated_categories)
 144         {
 145           rc->n_allocated_categories *= 2;
 146           rc->vals = (union value **)
 147             xrealloc (rc->vals, rc->n_allocated_categories
 148                       * sizeof (*(rc->vals)));
 149         }
 150       *(rc->vals + rc->n_categories) = v;
 151       rc->n_categories++;
 152     }
 153 }
 154
 155 /*
 156    Create a set of gsl_matrix's, each of whose rows correspond to
 157    values of a categorical variable. Since n categories have n-1
 158    degrees of freedom, the gsl_matrix is n-by-(n-1), with the first
 159    category encoded as the zero vector.
 160  */
 161 void
 162 cr_create_value_matrices (struct recoded_categorical_array *r)
 163 {
 164   size_t i;
 165   size_t row;
 166   size_t col;
 167   size_t n_rows;
 168   size_t n_cols;
 169
 170   for (i = 0; i < r->n_vars; i++)
 171     {
 172       n_rows = (*(r->a + i))->n_categories;
 173       n_cols = (*(r->a + i))->n_categories - 1;
 174       (*(r->a + i))->m = gsl_matrix_calloc (n_rows, n_cols);
 175       for (row = 1; row < n_rows; row++)
 176         {
 177           col = row - 1;
 178           gsl_matrix_set ((*(r->a + i))->m, row, col, 1.0);
 179         }
 180     }
 181 }
 182
 183 static size_t
 184 cr_value_to_subscript (const union value *val, struct recoded_categorical *cr)
 185 {
 186   const union value *v;
 187   size_t subscript;
 188   int different;
 189
 190   subscript = cr->n_categories - 1;
 191   while (subscript > 0)
 192     {
 193       v = *(cr->vals + subscript);
 194       different = compare_values (val, v, cr->v->width);
 195       if (!different)
 196         {
 197           return subscript;
 198         }
 199       subscript--;
 200     }
 201   return subscript;
 202 }
 203
 204 static const union value *
 205 cr_subscript_to_value (const size_t s, struct recoded_categorical *cr)
 206 {
 207   if (s < cr->n_categories)
 208     {
 209       return cr->vals[s];
 210     }
 211   else
 212     {
 213       return NULL;
 214     }
 215 }
 216
 217 /*
 218   Return the row of the matrix corresponding
 219   to the value v.
 220  */
 221 gsl_vector_const_view
 222 cr_value_to_vector (const union value * v, struct recoded_categorical * cr)
 223 {
 224   size_t row;
 225   row = cr_value_to_subscript (v, cr);
 226   return gsl_matrix_const_row (cr->m, row);
 227 }
 228
 229 /*
 230   Which element of a vector is equal to the value x?
 231  */
 232 static size_t
 233 cr_which_element_eq (const gsl_vector * vec, double x)
 234 {
 235   size_t i;
 236
 237   for (i = 0; i < vec->size; i++)
 238     {
 239       if (fabs (gsl_vector_get (vec, i) - x) < GSL_DBL_EPSILON)
 240         {
 241           return i;
 242         }
 243     }
 244   return CR_VALUE_NOT_FOUND;
 245 }
 246 static int
 247 cr_is_zero_vector (const gsl_vector * vec)
 248 {
 249   size_t i;
 250
 251   for (i = 0; i < vec->size; i++)
 252     {
 253       if (gsl_vector_get (vec, i) != 0.0)
 254         {
 255           return 0;
 256         }
 257     }
 258   return 1;
 259 }
 260
 261 /*
 262   Return the value corresponding to the vector.
 263   To avoid searching the matrix, this routine takes
 264   advantage of the fact that element (i,i+1) is 1
 265   when i is between 1 and cr->n_categories - 1 and
 266   i is 0 otherwise.
 267  */
 268 const union value *
 269 cr_vector_to_value (const gsl_vector * vec, struct recoded_categorical *cr)
 270 {
 271   size_t i;
 272
 273   i = cr_which_element_eq (vec, 1.0);
 274   if (i != CR_VALUE_NOT_FOUND)
 275     {
 276       return cr_subscript_to_value (i + 1, cr);
 277     }
 278   if (cr_is_zero_vector (vec))
 279     {
 280       return cr_subscript_to_value (0, cr);
 281     }
 282   return NULL;
 283 }
 284
 285 /*
 286   Given a variable, return a pointer to its recoded
 287   structure. BUSTED IN HERE.
 288  */
 289 struct recoded_categorical *
 290 cr_var_to_recoded_categorical (const struct variable *v,
 291                                struct recoded_categorical_array *ca)
 292 {
 293   struct recoded_categorical *rc;
 294   size_t i;
 295
 296   for (i = 0; i < ca->n_vars; i++)
 297     {
 298       rc = *(ca->a + i);
 299       if (rc->v->index == v->index)
 300         {
 301           return rc;
 302         }
 303     }
 304   return NULL;
 305 }
 306
 307 struct design_matrix *
 308 design_matrix_create (int n_variables,
 309                       const struct variable *v_variables[],
 310                       struct recoded_categorical_array *ca,
 311                       const size_t n_data)
 312 {
 313   struct design_matrix *dm;
 314   struct design_matrix_var *tmp;
 315   struct recoded_categorical *rc;
 316   const struct variable *v;
 317   size_t i;
 318   size_t n_cols = 0;
 319   size_t col;
 320
 321   dm = xmalloc (sizeof (*dm));
 322   dm->vars = xmalloc (n_variables * sizeof (struct variable *));
 323   dm->n_vars = n_variables;
 324
 325   for (i = 0; i < n_variables; i++)
 326     {
 327       v = v_variables[i];
 328       if (v->type == NUMERIC)
 329         {
 330           n_cols++;
 331         }
 332       else if (v->type == ALPHA)
 333         {
 334           assert (ca != NULL);
 335           rc = cr_var_to_recoded_categorical (v, ca);
 336           assert (rc != NULL);
 337           rc->first_column = n_cols;
 338           rc->last_column = rc->first_column + rc->n_categories - 2;
 339           n_cols += rc->n_categories - 1;
 340         }
 341     }
 342   dm->m = gsl_matrix_calloc (n_data, n_cols);
 343   dm->vars = xmalloc (dm->n_vars * sizeof (*(dm->vars)));
 344   assert (dm->vars != NULL);
 345   col = 0;
 346
 347   for (i = 0; i < n_variables; i++)
 348     {
 349       v = v_variables[i];
 350       (dm->vars[i]).v = v;
 351       if (v->type == NUMERIC)
 352         {
 353           tmp = &(dm->vars[col]);
 354           tmp->v = v;
 355           tmp->first_column = col;
 356           col++;
 357         }
 358       else if (v->type == ALPHA)
 359         {
 360           assert (ca != NULL);
 361           rc = cr_var_to_recoded_categorical (v, ca);
 362           assert (rc != NULL);
 363           tmp = &(dm->vars[col]);
 364           tmp->v = v;
 365           tmp->last_column = rc->last_column;
 366           col = rc->last_column + 1;
 367         }
 368     }
 369   return dm;
 370 }
 371
 372 void
 373 design_matrix_destroy (struct design_matrix *dm)
 374 {
 375   free (dm->vars);
 376   gsl_matrix_free (dm->m);
 377   free (dm);
 378 }
 379
 380 /*
 381   Return the index of the variable for the
 382   given column.
 383  */
 384 static const size_t
 385 design_matrix_col_to_var_index (const struct design_matrix *dm, size_t col)
 386 {
 387   size_t i;
 388   struct design_matrix_var v;
 389
 390   for (i = 0; i < dm->n_vars; i++)
 391     {
 392       v = dm->vars[i];
 393       if (v.first_column <= col && col <= v.last_column)
 394         return (v.v)->index;
 395     }
 396   return CR_INDEX_NOT_FOUND;
 397 }
 398
 399 /*
 400   Return a pointer to the variable whose values
 401   are stored in column col.
 402  */
 403 const struct variable *
 404 design_matrix_col_to_var (const struct design_matrix *dm, size_t col)
 405 {
 406   size_t index;
 407   size_t i;
 408   struct design_matrix_var dmv;
 409   const struct variable *v;
 410
 411   index = design_matrix_col_to_var_index (dm, col);
 412   for (i = 0; i < dm->n_vars; i++)
 413     {
 414       dmv = dm->vars[i];
 415       v = (dmv.v)->index;
 416       if (v->index == index)
 417         {
 418           return v;
 419         }
 420     }
 421   return NULL;
 422 }
 423
 424 static size_t
 425 cmp_dm_var_index (const struct design_matrix_var *dmv, size_t index)
 426 {
 427   if (dmv->v->index == index)
 428     return 1;
 429   return 0;
 430 }
 431
 432 /*
 433   Return the number of the first column storing the
 434   values for variable v.
 435  */
 436 size_t
 437 design_matrix_var_to_column (const struct design_matrix * dm,
 438                              const struct variable * v)
 439 {
 440   size_t i;
 441   struct design_matrix_var tmp;
 442
 443   for (i = 0; i < dm->n_vars; i++)
 444     {
 445       tmp = dm->vars[i];
 446       if (cmp_dm_var_index (&tmp, v->index))
 447         {
 448           return tmp.first_column;
 449         }
 450     }
 451   return CR_COLUMN_NOT_FOUND;
 452 }
 453
 454 /*
 455   Set the appropriate value in the design matrix,
 456   whether that value is from a categorical or numeric
 457   variable.
 458  */
 459 void
 460 design_matrix_set_categorical (struct design_matrix *dm, size_t row,
 461                                const struct variable *var,
 462                                const union value *val,
 463                                struct recoded_categorical *rc)
 464 {
 465   size_t col;
 466   double x;
 467
 468   assert (var->type == ALPHA);
 469   gsl_vector_const_view vec = cr_value_to_vector (val, rc);
 470
 471   /*
 472      Copying values here is not the 'most efficient' way,
 473      but it will work even if we change the vector encoding later.
 474    */
 475   for (col = rc->first_column; col <= rc->last_column; col++)
 476     {
 477       x = gsl_vector_get (&vec.vector, col);
 478       gsl_matrix_set (dm->m, row, col, x);
 479     }
 480 }
 481 void
 482 design_matrix_set_numeric (struct design_matrix *dm, size_t row,
 483                            const struct variable *var, const union value *val)
 484 {
 485   size_t col;
 486
 487   assert (var->type == NUMERIC);
 488   col = design_matrix_var_to_column ((const struct design_matrix *) dm, var);
 489   assert (col != CR_COLUMN_NOT_FOUND);
 490   gsl_matrix_set (dm->m, row, col, val->f);
 491 }