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