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 store values of a categorical
  22   variable, and to recode those values into binary vectors.
  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. For example, 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 #include <config.h>
  35 #include <stdlib.h>
  36 #include <error.h>
  37 #include "alloc.h"
  38 #include "error.h"
  39 #include "var.h"
  40 #include "cat.h"
  41 #include <string.h>
  42 #include <math.h>
  43 #include <gsl/gsl_machine.h>
  44 #include <gsl/gsl_vector.h>
  45 #include <gsl/gsl_matrix.h>
  46
  47 #define N_INITIAL_CATEGORIES 1
  48 #define CAT_COLUMN_NOT_FOUND -1
  49 #define CAT_VALUE_NOT_FOUND -2
  50 #define CAT_INDEX_NOT_FOUND -3
  51
  52 void
  53 cat_stored_values_create (struct variable *v)
  54 {
  55   if (v->obs_vals == NULL)
  56     {
  57       v->obs_vals = xmalloc (sizeof (*v->obs_vals));
  58       v->obs_vals->n_categories = 0;
  59       v->obs_vals->n_allocated_categories = N_INITIAL_CATEGORIES;
  60       v->obs_vals->vals =
  61         xnmalloc (N_INITIAL_CATEGORIES, sizeof *v->obs_vals->vals);
  62     }
  63 }
  64
  65 void
  66 cat_stored_values_destroy (struct variable *v)
  67 {
  68   assert (v != NULL);
  69   if (v->obs_vals != NULL)
  70     {
  71       free (v->obs_vals);
  72     }
  73 }
  74
  75 static size_t
  76 cat_value_find (const struct variable *v, const union value *val)
  77 {
  78   size_t i;
  79   const union value *candidate;
  80
  81   assert (val != NULL);
  82   assert (v != NULL);
  83   assert (v->obs_vals != NULL);
  84   for (i = 0; i < v->obs_vals->n_categories; i++)
  85     {
  86       candidate = v->obs_vals->vals + i;
  87       assert (candidate != NULL);
  88       if (!compare_values (candidate, val, v->width))
  89         {
  90           return i;
  91         }
  92     }
  93   return CAT_VALUE_NOT_FOUND;
  94 }
  95
  96 /*
  97    Add the new value unless it is already present.
  98  */
  99 void
 100 cat_value_update (struct variable *v, const union value *val)
 101 {
 102   struct cat_vals *cv;
 103
 104   if (v->type == ALPHA)
 105     {
 106       assert (val != NULL);
 107       assert (v != NULL);
 108       cv = v->obs_vals;
 109       if (cat_value_find (v, val) == CAT_VALUE_NOT_FOUND)
 110         {
 111           if (cv->n_categories >= cv->n_allocated_categories)
 112             {
 113               cv->n_allocated_categories *= 2;
 114               cv->vals = xnrealloc (cv->vals,
 115                                     cv->n_allocated_categories,
 116                                     sizeof *cv->vals);
 117             }
 118           cv->vals[cv->n_categories] = *val;
 119           cv->n_categories++;
 120         }
 121     }
 122 }
 123
 124 /*
 125   Return the subscript of the binary vector corresponding
 126   to this value.
 127  */
 128 size_t
 129 cat_value_to_subscript (const union value *val, struct variable *v)
 130 {
 131   const union value *val2;
 132   size_t subscript;
 133   int different;
 134
 135   assert (v != NULL);
 136   assert (val != NULL);
 137   assert (v->obs_vals != NULL);
 138   subscript = v->obs_vals->n_categories - 1;
 139   while (subscript > 0)
 140     {
 141       val2 = v->obs_vals->vals + subscript;
 142       assert (val2 != NULL);
 143       different = compare_values (val, val2, v->width);
 144       if (!different)
 145         {
 146           return subscript;
 147         }
 148       subscript--;
 149     }
 150   return subscript;
 151 }
 152
 153 union value *
 154 cat_subscript_to_value (const size_t s, struct variable *v)
 155 {
 156   assert (v->obs_vals != NULL);
 157   if (s < v->obs_vals->n_categories)
 158     {
 159       return (v->obs_vals->vals + s);
 160     }
 161   else
 162     {
 163       return NULL;
 164     }
 165 }
 166
 167 /*
 168   Which element of a vector is equal to the value x?
 169  */
 170 static size_t
 171 cat_which_element_eq (const gsl_vector * vec, double x)
 172 {
 173   size_t i;
 174
 175   for (i = 0; i < vec->size; i++)
 176     {
 177       if (fabs (gsl_vector_get (vec, i) - x) < GSL_DBL_EPSILON)
 178         {
 179           return i;
 180         }
 181     }
 182   return CAT_VALUE_NOT_FOUND;
 183 }
 184 static int
 185 cat_is_zero_vector (const gsl_vector * vec)
 186 {
 187   size_t i;
 188
 189   for (i = 0; i < vec->size; i++)
 190     {
 191       if (gsl_vector_get (vec, i) != 0.0)
 192         {
 193           return 0;
 194         }
 195     }
 196   return 1;
 197 }
 198
 199 /*
 200   Return the value of v corresponding to the vector vec.
 201  */
 202 union value *
 203 cat_vector_to_value (const gsl_vector * vec, struct variable *v)
 204 {
 205   size_t i;
 206
 207   i = cat_which_element_eq (vec, 1.0);
 208   if (i != CAT_VALUE_NOT_FOUND)
 209     {
 210       return cat_subscript_to_value (i + 1, v);
 211     }
 212   if (cat_is_zero_vector (vec))
 213     {
 214       return cat_subscript_to_value (0, v);
 215     }
 216   return NULL;
 217 }
 218
 219 struct design_matrix *
 220 design_matrix_create (int n_variables,
 221                       const struct variable *v_variables[],
 222                       const size_t n_data)
 223 {
 224   struct design_matrix *dm;
 225   const struct variable *v;
 226   size_t i;
 227   size_t n_cols = 0;
 228   size_t col;
 229
 230   dm = xmalloc (sizeof *dm);
 231   dm->vars = xnmalloc (n_variables, sizeof *dm->vars);
 232   dm->n_vars = n_variables;
 233
 234   for (i = 0; i < n_variables; i++)
 235     {
 236       v = v_variables[i];
 237       assert ((dm->vars + i) != NULL);
 238       (dm->vars + i)->v = v;    /* Allows us to look up the variable from
 239                                    the design matrix. */
 240       (dm->vars + i)->first_column = n_cols;
 241       if (v->type == NUMERIC)
 242         {
 243           n_cols++;
 244           (dm->vars + i)->last_column = n_cols;
 245         }
 246       else if (v->type == ALPHA)
 247         {
 248           assert (v->obs_vals != NULL);
 249           (dm->vars + i)->last_column =
 250             (dm->vars + i)->first_column + v->obs_vals->n_categories - 2;
 251           n_cols += v->obs_vals->n_categories - 1;
 252         }
 253     }
 254   dm->m = gsl_matrix_calloc (n_data, n_cols);
 255   col = 0;
 256
 257   return dm;
 258 }
 259
 260 void
 261 design_matrix_destroy (struct design_matrix *dm)
 262 {
 263   free (dm->vars);
 264   gsl_matrix_free (dm->m);
 265   free (dm);
 266 }
 267
 268 /*
 269   Return the index of the variable for the
 270   given column.
 271  */
 272 static size_t
 273 design_matrix_col_to_var_index (const struct design_matrix *dm, size_t col)
 274 {
 275   size_t i;
 276   struct design_matrix_var v;
 277
 278   for (i = 0; i < dm->n_vars; i++)
 279     {
 280       v = dm->vars[i];
 281       if (v.first_column <= col && col <= v.last_column)
 282         return (v.v)->index;
 283     }
 284   return CAT_INDEX_NOT_FOUND;
 285 }
 286
 287 /*
 288   Return a pointer to the variable whose values
 289   are stored in column col.
 290  */
 291 struct variable *
 292 design_matrix_col_to_var (const struct design_matrix *dm, size_t col)
 293 {
 294   size_t index;
 295   size_t i;
 296   struct design_matrix_var dmv;
 297
 298   index = design_matrix_col_to_var_index (dm, col);
 299   for (i = 0; i < dm->n_vars; i++)
 300     {
 301       dmv = dm->vars[i];
 302       if ((dmv.v)->index == index)
 303         {
 304           return (struct variable *) dmv.v;
 305         }
 306     }
 307   return NULL;
 308 }
 309
 310 static size_t
 311 cmp_dm_var_index (const struct design_matrix_var *dmv, size_t index)
 312 {
 313   if (dmv->v->index == index)
 314     return 1;
 315   return 0;
 316 }
 317
 318 /*
 319   Return the number of the first column which holds the
 320   values for variable v.
 321  */
 322 size_t
 323 design_matrix_var_to_column (const struct design_matrix * dm,
 324                              const struct variable * v)
 325 {
 326   size_t i;
 327   struct design_matrix_var tmp;
 328
 329   for (i = 0; i < dm->n_vars; i++)
 330     {
 331       tmp = dm->vars[i];
 332       if (cmp_dm_var_index (&tmp, v->index))
 333         {
 334           return tmp.first_column;
 335         }
 336     }
 337   return CAT_COLUMN_NOT_FOUND;
 338 }
 339
 340 /* Last column. */
 341 static size_t
 342 dm_var_to_last_column (const struct design_matrix *dm,
 343                        const struct variable *v)
 344 {
 345   size_t i;
 346   struct design_matrix_var tmp;
 347
 348   for (i = 0; i < dm->n_vars; i++)
 349     {
 350       tmp = dm->vars[i];
 351       if (cmp_dm_var_index (&tmp, v->index))
 352         {
 353           return tmp.last_column;
 354         }
 355     }
 356   return CAT_COLUMN_NOT_FOUND;
 357 }
 358
 359 /*
 360   Set the appropriate value in the design matrix,
 361   whether that value is from a categorical or numeric
 362   variable. For a categorical variable, only the usual
 363   binary encoding is allowed.
 364  */
 365 void
 366 design_matrix_set_categorical (struct design_matrix *dm, size_t row,
 367                                const struct variable *var,
 368                                const union value *val)
 369 {
 370   size_t col;
 371   size_t is_one;
 372   size_t fc;
 373   size_t lc;
 374   double entry;
 375
 376   assert (var->type == ALPHA);
 377   fc = design_matrix_var_to_column (dm, var);
 378   lc = dm_var_to_last_column (dm, var);
 379   assert (lc != CAT_COLUMN_NOT_FOUND);
 380   assert (fc != CAT_COLUMN_NOT_FOUND);
 381   is_one = fc + cat_value_find (var, val);
 382   for (col = fc; col <= lc; col++)
 383     {
 384       entry = (col == is_one) ? 1.0 : 0.0;
 385       gsl_matrix_set (dm->m, row, col, entry);
 386     }
 387 }
 388 void
 389 design_matrix_set_numeric (struct design_matrix *dm, size_t row,
 390                            const struct variable *var, const union value *val)
 391 {
 392   size_t col;
 393
 394   assert (var->type == NUMERIC);
 395   col = design_matrix_var_to_column ((const struct design_matrix *) dm, var);
 396   assert (col != CAT_COLUMN_NOT_FOUND);
 397   gsl_matrix_set (dm->m, row, col, val->f);
 398 }