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_const_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 = (union value **) 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 = (union value **)
 150             xrealloc (rc->vals, rc->n_allocated_categories
 151                       * sizeof (*(rc->vals)));
 152         }
 153       *(rc->vals + rc->n_categories) = v;
 154       rc->n_categories++;
 155     }
 156 }
 157
 158 /*
 159    Create a set of gsl_matrix's, each of whose rows correspond to
 160    values of a categorical variable. Since n categories have n-1
 161    degrees of freedom, the gsl_matrix is n-by-(n-1), with the first
 162    category encoded as the zero vector.
 163  */
 164 void
 165 cr_create_value_matrices (struct recoded_categorical_array *r)
 166 {
 167   size_t i;
 168   size_t row;
 169   size_t col;
 170   size_t n_rows;
 171   size_t n_cols;
 172
 173   for (i = 0; i < r->n_vars; i++)
 174     {
 175       n_rows = (*(r->a + i))->n_categories;
 176       n_cols = (*(r->a + i))->n_categories - 1;
 177       (*(r->a + i))->m = gsl_matrix_calloc (n_rows, n_cols);
 178       for (row = 1; row < n_rows; row++)
 179         {
 180           col = row - 1;
 181           gsl_matrix_set ((*(r->a + i))->m, row, col, 1.0);
 182         }
 183     }
 184 }
 185
 186 static size_t
 187 cr_value_to_subscript (const union value *val, struct recoded_categorical *cr)
 188 {
 189   const union value *v;
 190   size_t subscript;
 191   int different;
 192
 193   subscript = cr->n_categories - 1;
 194   while (subscript > 0)
 195     {
 196       v = *(cr->vals + subscript);
 197       different = compare_values (val, v, cr->v->width);
 198       if (!different)
 199         {
 200           return subscript;
 201         }
 202       subscript--;
 203     }
 204   return subscript;
 205 }
 206
 207 static const union value *
 208 cr_subscript_to_value (const size_t s, struct recoded_categorical *cr)
 209 {
 210   if (s < cr->n_categories)
 211     {
 212       return cr->vals[s];
 213     }
 214   else
 215     {
 216       return NULL;
 217     }
 218 }
 219
 220 /*
 221   Return the row of the matrix corresponding
 222   to the value v.
 223  */
 224 static gsl_vector_const_view
 225 cr_value_to_vector (const union value * v, struct recoded_categorical * cr)
 226 {
 227   size_t row;
 228   row = cr_value_to_subscript (v, cr);
 229   return gsl_matrix_const_row (cr->m, row);
 230 }
 231
 232 /*
 233   Which element of a vector is equal to the value x?
 234  */
 235 static size_t
 236 cr_which_element_eq (const gsl_vector * vec, double x)
 237 {
 238   size_t i;
 239
 240   for (i = 0; i < vec->size; i++)
 241     {
 242       if (fabs (gsl_vector_get (vec, i) - x) < GSL_DBL_EPSILON)
 243         {
 244           return i;
 245         }
 246     }
 247   return CR_VALUE_NOT_FOUND;
 248 }
 249 static int
 250 cr_is_zero_vector (const gsl_vector * vec)
 251 {
 252   size_t i;
 253
 254   for (i = 0; i < vec->size; i++)
 255     {
 256       if (gsl_vector_get (vec, i) != 0.0)
 257         {
 258           return 0;
 259         }
 260     }
 261   return 1;
 262 }
 263
 264 /*
 265   Return the value corresponding to the vector.
 266   To avoid searching the matrix, this routine takes
 267   advantage of the fact that element (i,i+1) is 1
 268   when i is between 1 and cr->n_categories - 1 and
 269   i is 0 otherwise.
 270  */
 271 const union value *
 272 cr_vector_to_value (const gsl_vector * vec, struct recoded_categorical *cr)
 273 {
 274   size_t i;
 275
 276   i = cr_which_element_eq (vec, 1.0);
 277   if (i != CR_VALUE_NOT_FOUND)
 278     {
 279       return cr_subscript_to_value (i + 1, cr);
 280     }
 281   if (cr_is_zero_vector (vec))
 282     {
 283       return cr_subscript_to_value (0, cr);
 284     }
 285   return NULL;
 286 }
 287
 288 /*
 289   Given a variable, return a pointer to its recoded
 290   structure. BUSTED IN HERE.
 291  */
 292 struct recoded_categorical *
 293 cr_var_to_recoded_categorical (const struct variable *v,
 294                                struct recoded_categorical_array *ca)
 295 {
 296   struct recoded_categorical *rc;
 297   size_t i;
 298
 299   for (i = 0; i < ca->n_vars; i++)
 300     {
 301       rc = *(ca->a + i);
 302       if (rc->v->index == v->index)
 303         {
 304           return rc;
 305         }
 306     }
 307   return NULL;
 308 }
 309
 310 struct design_matrix *
 311 design_matrix_create (int n_variables,
 312                       const struct variable *v_variables[],
 313                       struct recoded_categorical_array *ca,
 314                       const size_t n_data)
 315 {
 316   struct design_matrix *dm;
 317   struct design_matrix_var *tmp;
 318   struct recoded_categorical *rc;
 319   const struct variable *v;
 320   size_t i;
 321   size_t n_cols = 0;
 322   size_t col;
 323
 324   dm = xmalloc (sizeof (*dm));
 325   dm->vars = xmalloc (n_variables * sizeof (struct variable *));
 326   dm->n_vars = n_variables;
 327
 328   for (i = 0; i < n_variables; i++)
 329     {
 330       v = v_variables[i];
 331       if (v->type == NUMERIC)
 332         {
 333           n_cols++;
 334         }
 335       else if (v->type == ALPHA)
 336         {
 337           assert (ca != NULL);
 338           rc = cr_var_to_recoded_categorical (v, ca);
 339           assert (rc != NULL);
 340           rc->first_column = n_cols;
 341           rc->last_column = rc->first_column + rc->n_categories - 2;
 342           n_cols += rc->n_categories - 1;
 343         }
 344     }
 345   dm->m = gsl_matrix_calloc (n_data, n_cols);
 346   dm->vars = xmalloc (dm->n_vars * sizeof (*(dm->vars)));
 347   assert (dm->vars != NULL);
 348   col = 0;
 349
 350   for (i = 0; i < n_variables; i++)
 351     {
 352       v = v_variables[i];
 353       (dm->vars[i]).v = v;
 354       if (v->type == NUMERIC)
 355         {
 356           tmp = &(dm->vars[col]);
 357           tmp->v = v;
 358           tmp->first_column = col;
 359           col++;
 360         }
 361       else if (v->type == ALPHA)
 362         {
 363           assert (ca != NULL);
 364           rc = cr_var_to_recoded_categorical (v, ca);
 365           assert (rc != NULL);
 366           tmp = &(dm->vars[col]);
 367           tmp->v = v;
 368           tmp->last_column = rc->last_column;
 369           col = rc->last_column + 1;
 370         }
 371     }
 372   return dm;
 373 }
 374
 375 void
 376 design_matrix_destroy (struct design_matrix *dm)
 377 {
 378   free (dm->vars);
 379   gsl_matrix_free (dm->m);
 380   free (dm);
 381 }
 382
 383 /*
 384   Return the index of the variable for the
 385   given column.
 386  */
 387 static const size_t
 388 design_matrix_col_to_var_index (const struct design_matrix *dm, size_t col)
 389 {
 390   size_t i;
 391   struct design_matrix_var v;
 392
 393   for (i = 0; i < dm->n_vars; i++)
 394     {
 395       v = dm->vars[i];
 396       if (v.first_column <= col && col <= v.last_column)
 397         return (v.v)->index;
 398     }
 399   return CR_INDEX_NOT_FOUND;
 400 }
 401
 402 /*
 403   Return a pointer to the variable whose values
 404   are stored in column col.
 405  */
 406 const struct variable *
 407 design_matrix_col_to_var (const struct design_matrix *dm, size_t col)
 408 {
 409   size_t index;
 410   size_t i;
 411   struct design_matrix_var dmv;
 412   const struct variable *v;
 413
 414   index = design_matrix_col_to_var_index (dm, col);
 415   for (i = 0; i < dm->n_vars; i++)
 416     {
 417       dmv = dm->vars[i];
 418       v = (dmv.v)->index;
 419       if (v->index == index)
 420         {
 421           return v;
 422         }
 423     }
 424   return NULL;
 425 }
 426
 427 static size_t
 428 cmp_dm_var_index (const struct design_matrix_var *dmv, size_t index)
 429 {
 430   if (dmv->v->index == index)
 431     return 1;
 432   return 0;
 433 }
 434
 435 /*
 436   Return the number of the first column storing the
 437   values for variable v.
 438  */
 439 size_t
 440 design_matrix_var_to_column (const struct design_matrix * dm,
 441                              const struct variable * v)
 442 {
 443   size_t i;
 444   struct design_matrix_var tmp;
 445
 446   for (i = 0; i < dm->n_vars; i++)
 447     {
 448       tmp = dm->vars[i];
 449       if (cmp_dm_var_index (&tmp, v->index))
 450         {
 451           return tmp.first_column;
 452         }
 453     }
 454   return CR_COLUMN_NOT_FOUND;
 455 }
 456
 457 /*
 458   Set the appropriate value in the design matrix,
 459   whether that value is from a categorical or numeric
 460   variable.
 461  */
 462 void
 463 design_matrix_set_categorical (struct design_matrix *dm, size_t row,
 464                                const struct variable *var,
 465                                const union value *val,
 466                                struct recoded_categorical *rc)
 467 {
 468   size_t col;
 469   double x;
 470
 471   assert (var->type == ALPHA);
 472   gsl_vector_const_view vec = cr_value_to_vector (val, rc);
 473
 474   /*
 475      Copying values here is not the 'most efficient' way,
 476      but it will work even if we change the vector encoding later.
 477    */
 478   for (col = rc->first_column; col <= rc->last_column; col++)
 479     {
 480       x = gsl_vector_get (&vec.vector, col);
 481       gsl_matrix_set (dm->m, row, col, x);
 482     }
 483 }
 484 void
 485 design_matrix_set_numeric (struct design_matrix *dm, size_t row,
 486                            const struct variable *var, const union value *val)
 487 {
 488   size_t col;
 489
 490   assert (var->type == NUMERIC);
 491   col = design_matrix_var_to_column ((const struct design_matrix *) dm, var);
 492   assert (col != CR_COLUMN_NOT_FOUND);
 493   gsl_matrix_set (dm->m, row, col, val->f);
 494 }