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 static gsl_vector_view cr_value_to_vector (const union value *,
  57                                            struct recoded_categorical *);
  58
  59 struct recoded_categorical *
  60 cr_recoded_categorical_create (const struct variable *v)
  61 {
  62   struct recoded_categorical *rc;
  63
  64   rc = xmalloc (sizeof (*rc));
  65   rc->v = v;
  66   rc->n_categories = 0;
  67   rc->n_allocated_categories = N_INITIAL_CATEGORIES;
  68   rc->vals = xmalloc (N_INITIAL_CATEGORIES *
  69                       sizeof (*rc->vals));
  70
  71   return rc;
  72 }
  73
  74 void
  75 cr_recoded_categorical_destroy (struct recoded_categorical *r)
  76 {
  77   free (r->vals);
  78   free (r);
  79 }
  80
  81 struct recoded_categorical_array *
  82 cr_recoded_cat_ar_create (int n_variables, struct variable *v_variables[])
  83 {
  84   size_t n_categoricals = 0;
  85   size_t i;
  86   struct recoded_categorical_array *ca;
  87   struct variable *v;
  88
  89   ca = (struct recoded_categorical_array *) xmalloc (sizeof (*ca));
  90   for (i = 0; i < n_variables; i++)
  91     {
  92       v = v_variables[i];
  93       if (v->type == ALPHA)
  94         {
  95           n_categoricals++;
  96         }
  97     }
  98   ca->n_vars = n_categoricals;
  99   ca->a = xmalloc (n_categoricals * sizeof (*(ca->a)));
 100   for (i = 0; i < n_categoricals; i++)
 101     {
 102       *(ca->a + i) = cr_recoded_categorical_create (v_variables[i]);
 103     }
 104
 105   return ca;
 106 }
 107
 108 int
 109 cr_free_recoded_array (struct recoded_categorical_array *r)
 110 {
 111   int rc = 0;
 112   size_t i;
 113
 114   for (i = 0; i < r->n_vars; i++)
 115     {
 116       cr_recoded_categorical_destroy (*(r->a + i));
 117     }
 118   return rc;
 119 }
 120
 121 static size_t
 122 cr_value_find (struct recoded_categorical *rc, const union value *v)
 123 {
 124   size_t i;
 125   const union value *val;
 126
 127   for (i = 0; i < rc->n_categories; i++)
 128     {
 129       val = rc->vals + i;
 130       if (!compare_values (val, v, rc->v->width))
 131         {
 132           return i;
 133         }
 134     }
 135   return CR_VALUE_NOT_FOUND;
 136 }
 137
 138 /*
 139    Add the new value unless it is already present.
 140  */
 141 void
 142 cr_value_update (struct recoded_categorical *rc, const union value *v)
 143 {
 144   if (cr_value_find (rc, v) == CR_VALUE_NOT_FOUND)
 145     {
 146       if (rc->n_categories >= rc->n_allocated_categories)
 147         {
 148           rc->n_allocated_categories *= 2;
 149           rc->vals =  xrealloc (rc->vals, rc->n_allocated_categories
 150                                 * sizeof (*(rc->vals)));
 151         }
 152       rc->vals[rc->n_categories] = *v;
 153       rc->n_categories++;
 154     }
 155 }
 156
 157 /*
 158    Create a set of gsl_matrix's, each of whose rows correspond to
 159    values of a categorical variable. Since n categories have n-1
 160    degrees of freedom, the gsl_matrix is n-by-(n-1), with the first
 161    category encoded as the zero vector.
 162  */
 163 void
 164 cr_create_value_matrices (struct recoded_categorical_array *r)
 165 {
 166   size_t i;
 167   size_t row;
 168   size_t col;
 169   size_t n_rows;
 170   size_t n_cols;
 171
 172   for (i = 0; i < r->n_vars; i++)
 173     {
 174       n_rows = (*(r->a + i))->n_categories;
 175       n_cols = (*(r->a + i))->n_categories - 1;
 176       (*(r->a + i))->m = gsl_matrix_calloc (n_rows, n_cols);
 177       for (row = 1; row < n_rows; row++)
 178         {
 179           col = row - 1;
 180           gsl_matrix_set ((*(r->a + i))->m, row, col, 1.0);
 181         }
 182     }
 183 }
 184
 185 static size_t
 186 cr_value_to_subscript (const union value *val, struct recoded_categorical *cr)
 187 {
 188   const union value *v;
 189   size_t subscript;
 190   int different;
 191
 192   subscript = cr->n_categories - 1;
 193   while (subscript > 0)
 194     {
 195       v = cr->vals + subscript;
 196       different = compare_values (val, v, cr->v->width);
 197       if (!different)
 198         {
 199           return subscript;
 200         }
 201       subscript--;
 202     }
 203   return subscript;
 204 }
 205
 206 static union value *
 207 cr_subscript_to_value (const size_t s, struct recoded_categorical *cr)
 208 {
 209   if (s < cr->n_categories)
 210     {
 211       return (cr->vals + s);
 212     }
 213   else
 214     {
 215       return NULL;
 216     }
 217 }
 218
 219 /*
 220   Return the row of the matrix corresponding
 221   to the value v.
 222  */
 223 static gsl_vector_view
 224 cr_value_to_vector (const union value * v, struct recoded_categorical * cr)
 225 {
 226   size_t row;
 227   row = cr_value_to_subscript (v, cr);
 228   return gsl_matrix_row (cr->m, row);
 229 }
 230
 231 /*
 232   Which element of a vector is equal to the value x?
 233  */
 234 static size_t
 235 cr_which_element_eq (const gsl_vector * vec, double x)
 236 {
 237   size_t i;
 238
 239   for (i = 0; i < vec->size; i++)
 240     {
 241       if (fabs (gsl_vector_get (vec, i) - x) < GSL_DBL_EPSILON)
 242         {
 243           return i;
 244         }
 245     }
 246   return CR_VALUE_NOT_FOUND;
 247 }
 248 static int
 249 cr_is_zero_vector (const gsl_vector * vec)
 250 {
 251   size_t i;
 252
 253   for (i = 0; i < vec->size; i++)
 254     {
 255       if (gsl_vector_get (vec, i) != 0.0)
 256         {
 257           return 0;
 258         }
 259     }
 260   return 1;
 261 }
 262
 263 /*
 264   Return the value corresponding to the vector.
 265   To avoid searching the matrix, this routine takes
 266   advantage of the fact that element (i,i+1) is 1
 267   when i is between 1 and cr->n_categories - 1 and
 268   i is 0 otherwise.
 269  */
 270 union value *
 271 cr_vector_to_value (const gsl_vector * vec, struct recoded_categorical *cr)
 272 {
 273   size_t i;
 274
 275   i = cr_which_element_eq (vec, 1.0);
 276   if (i != CR_VALUE_NOT_FOUND)
 277     {
 278       return cr_subscript_to_value (i + 1, cr);
 279     }
 280   if (cr_is_zero_vector (vec))
 281     {
 282       return cr_subscript_to_value (0, cr);
 283     }
 284   return NULL;
 285 }
 286
 287 /*
 288   Given a variable, return a pointer to its recoded
 289   structure. BUSTED IN HERE.
 290  */
 291 struct recoded_categorical *
 292 cr_var_to_recoded_categorical (const struct variable *v,
 293                                struct recoded_categorical_array *ca)
 294 {
 295   struct recoded_categorical *rc;
 296   size_t i;
 297
 298   for (i = 0; i < ca->n_vars; i++)
 299     {
 300       rc = *(ca->a + i);
 301       if (rc->v->index == v->index)
 302         {
 303           return rc;
 304         }
 305     }
 306   return NULL;
 307 }
 308
 309 struct design_matrix *
 310 design_matrix_create (int n_variables,
 311                       const struct variable *v_variables[],
 312                       struct recoded_categorical_array *ca,
 313                       const size_t n_data)
 314 {
 315   struct design_matrix *dm;
 316   struct design_matrix_var *tmp;
 317   struct recoded_categorical *rc;
 318   const struct variable *v;
 319   size_t i;
 320   size_t n_cols = 0;
 321   size_t col;
 322
 323   dm = xmalloc (sizeof (*dm));
 324   dm->vars = xmalloc (n_variables * sizeof (struct variable *));
 325   dm->n_vars = n_variables;
 326
 327   for (i = 0; i < n_variables; i++)
 328     {
 329       v = v_variables[i];
 330       if (v->type == NUMERIC)
 331         {
 332           n_cols++;
 333         }
 334       else if (v->type == ALPHA)
 335         {
 336           assert (ca != NULL);
 337           rc = cr_var_to_recoded_categorical (v, ca);
 338           assert (rc != NULL);
 339           rc->first_column = n_cols;
 340           rc->last_column = rc->first_column + rc->n_categories - 2;
 341           n_cols += rc->n_categories - 1;
 342         }
 343     }
 344   dm->m = gsl_matrix_calloc (n_data, n_cols);
 345   dm->vars = xmalloc (dm->n_vars * sizeof (*(dm->vars)));
 346   assert (dm->vars != NULL);
 347   col = 0;
 348
 349   for (i = 0; i < n_variables; i++)
 350     {
 351       v = v_variables[i];
 352       (dm->vars[i]).v = v;
 353       if (v->type == NUMERIC)
 354         {
 355           tmp = &(dm->vars[col]);
 356           tmp->v = v;
 357           tmp->first_column = col;
 358           col++;
 359         }
 360       else if (v->type == ALPHA)
 361         {
 362           assert (ca != NULL);
 363           rc = cr_var_to_recoded_categorical (v, ca);
 364           assert (rc != NULL);
 365           tmp = &(dm->vars[col]);
 366           tmp->v = v;
 367           tmp->last_column = rc->last_column;
 368           col = rc->last_column + 1;
 369         }
 370     }
 371   return dm;
 372 }
 373
 374 void
 375 design_matrix_destroy (struct design_matrix *dm)
 376 {
 377   free (dm->vars);
 378   gsl_matrix_free (dm->m);
 379   free (dm);
 380 }
 381
 382 /*
 383   Return the index of the variable for the
 384   given column.
 385  */
 386 static size_t
 387 design_matrix_col_to_var_index (const struct design_matrix *dm, size_t col)
 388 {
 389   size_t i;
 390   struct design_matrix_var v;
 391
 392   for (i = 0; i < dm->n_vars; i++)
 393     {
 394       v = dm->vars[i];
 395       if (v.first_column <= col && col <= v.last_column)
 396         return (v.v)->index;
 397     }
 398   return CR_INDEX_NOT_FOUND;
 399 }
 400
 401 /*
 402   Return a pointer to the variable whose values
 403   are stored in column col.
 404  */
 405 struct variable *
 406 design_matrix_col_to_var (const struct design_matrix *dm, size_t col)
 407 {
 408   size_t index;
 409   size_t i;
 410   struct design_matrix_var dmv;
 411
 412   index = design_matrix_col_to_var_index (dm, col);
 413   for (i = 0; i < dm->n_vars; i++)
 414     {
 415       dmv = dm->vars[i];
 416       if ((dmv.v)->index == index)
 417         {
 418           return dmv.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   const gsl_vector_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 }