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           tmp->last_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. BUG IN HERE
 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 (struct variable *) 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 }