src/math/design-matrix.c

   1 /* PSPP - Creates design-matrices.
   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   Create design matrices for procedures that need them.
  22 */
  23 #include <config.h>
  24 #include <stdlib.h>
  25 #include <message.h>
  26 #include "alloc.h"
  27 #include "message.h"
  28 #include "variable.h"
  29 #include "category.h"
  30 #include "design-matrix.h"
  31 #include <string.h>
  32 #include <math.h>
  33 #include <gsl/gsl_machine.h>
  34 #include <gsl/gsl_vector.h>
  35 #include <gsl/gsl_matrix.h>
  36
  37 #define DM_COLUMN_NOT_FOUND -1
  38 #define DM_INDEX_NOT_FOUND -3
  39
  40 /*
  41   Which element of a vector is equal to the value x?
  42  */
  43 static size_t
  44 cat_which_element_eq (const gsl_vector * vec, double x)
  45 {
  46   size_t i;
  47
  48   for (i = 0; i < vec->size; i++)
  49     {
  50       if (fabs (gsl_vector_get (vec, i) - x) < GSL_DBL_EPSILON)
  51         {
  52           return i;
  53         }
  54     }
  55   return CAT_VALUE_NOT_FOUND;
  56 }
  57 static int
  58 cat_is_zero_vector (const gsl_vector * vec)
  59 {
  60   size_t i;
  61
  62   for (i = 0; i < vec->size; i++)
  63     {
  64       if (gsl_vector_get (vec, i) != 0.0)
  65         {
  66           return 0;
  67         }
  68     }
  69   return 1;
  70 }
  71
  72 /*
  73   Return the value of v corresponding to the vector vec.
  74  */
  75 union value *
  76 cat_vector_to_value (const gsl_vector * vec, struct variable *v)
  77 {
  78   size_t i;
  79
  80   i = cat_which_element_eq (vec, 1.0);
  81   if (i != CAT_VALUE_NOT_FOUND)
  82     {
  83       return cat_subscript_to_value (i + 1, v);
  84     }
  85   if (cat_is_zero_vector (vec))
  86     {
  87       return cat_subscript_to_value (0, v);
  88     }
  89   return NULL;
  90 }
  91
  92 struct design_matrix *
  93 design_matrix_create (int n_variables,
  94                       const struct variable *v_variables[],
  95                       const size_t n_data)
  96 {
  97   struct design_matrix *dm;
  98   const struct variable *v;
  99   size_t i;
 100   size_t n_cols = 0;
 101   size_t col;
 102
 103   dm = xmalloc (sizeof *dm);
 104   dm->vars = xnmalloc (n_variables, sizeof *dm->vars);
 105   dm->n_vars = n_variables;
 106
 107   for (i = 0; i < n_variables; i++)
 108     {
 109       v = v_variables[i];
 110       assert ((dm->vars + i) != NULL);
 111       (dm->vars + i)->v = v;    /* Allows us to look up the variable from
 112                                    the design matrix. */
 113       (dm->vars + i)->first_column = n_cols;
 114       if (v->type == NUMERIC)
 115         {
 116           (dm->vars + i)->last_column = n_cols;
 117           n_cols++;
 118         }
 119       else if (v->type == ALPHA)
 120         {
 121           assert (v->obs_vals != NULL);
 122           (dm->vars + i)->last_column =
 123             (dm->vars + i)->first_column + v->obs_vals->n_categories - 2;
 124           n_cols += v->obs_vals->n_categories - 1;
 125         }
 126     }
 127   dm->m = gsl_matrix_calloc (n_data, n_cols);
 128   col = 0;
 129
 130   return dm;
 131 }
 132
 133 void
 134 design_matrix_destroy (struct design_matrix *dm)
 135 {
 136   free (dm->vars);
 137   gsl_matrix_free (dm->m);
 138   free (dm);
 139 }
 140
 141 /*
 142   Return the index of the variable for the
 143   given column.
 144  */
 145 static size_t
 146 design_matrix_col_to_var_index (const struct design_matrix *dm, size_t col)
 147 {
 148   size_t i;
 149   struct design_matrix_var v;
 150
 151   for (i = 0; i < dm->n_vars; i++)
 152     {
 153       v = dm->vars[i];
 154       if (v.first_column <= col && col <= v.last_column)
 155         return (v.v)->index;
 156     }
 157   return DM_INDEX_NOT_FOUND;
 158 }
 159
 160 /*
 161   Return a pointer to the variable whose values
 162   are stored in column col.
 163  */
 164 struct variable *
 165 design_matrix_col_to_var (const struct design_matrix *dm, size_t col)
 166 {
 167   size_t index;
 168   size_t i;
 169   struct design_matrix_var dmv;
 170
 171   index = design_matrix_col_to_var_index (dm, col);
 172   for (i = 0; i < dm->n_vars; i++)
 173     {
 174       dmv = dm->vars[i];
 175       if ((dmv.v)->index == index)
 176         {
 177           return (struct variable *) dmv.v;
 178         }
 179     }
 180   return NULL;
 181 }
 182
 183 static size_t
 184 cmp_dm_var_index (const struct design_matrix_var *dmv, size_t index)
 185 {
 186   if (dmv->v->index == index)
 187     return 1;
 188   return 0;
 189 }
 190
 191 /*
 192   Return the number of the first column which holds the
 193   values for variable v.
 194  */
 195 size_t
 196 design_matrix_var_to_column (const struct design_matrix * dm,
 197                              const struct variable * v)
 198 {
 199   size_t i;
 200   struct design_matrix_var tmp;
 201
 202   for (i = 0; i < dm->n_vars; i++)
 203     {
 204       tmp = dm->vars[i];
 205       if (cmp_dm_var_index (&tmp, v->index))
 206         {
 207           return tmp.first_column;
 208         }
 209     }
 210   return DM_COLUMN_NOT_FOUND;
 211 }
 212
 213 /* Last column. */
 214 static size_t
 215 dm_var_to_last_column (const struct design_matrix *dm,
 216                        const struct variable *v)
 217 {
 218   size_t i;
 219   struct design_matrix_var tmp;
 220
 221   for (i = 0; i < dm->n_vars; i++)
 222     {
 223       tmp = dm->vars[i];
 224       if (cmp_dm_var_index (&tmp, v->index))
 225         {
 226           return tmp.last_column;
 227         }
 228     }
 229   return DM_COLUMN_NOT_FOUND;
 230 }
 231
 232 /*
 233   Set the appropriate value in the design matrix,
 234   whether that value is from a categorical or numeric
 235   variable. For a categorical variable, only the usual
 236   binary encoding is allowed.
 237  */
 238 void
 239 design_matrix_set_categorical (struct design_matrix *dm, size_t row,
 240                                const struct variable *var,
 241                                const union value *val)
 242 {
 243   size_t col;
 244   size_t is_one;
 245   size_t fc;
 246   size_t lc;
 247   double entry;
 248
 249   assert (var->type == ALPHA);
 250   fc = design_matrix_var_to_column (dm, var);
 251   lc = dm_var_to_last_column (dm, var);
 252   assert (lc != DM_COLUMN_NOT_FOUND);
 253   assert (fc != DM_COLUMN_NOT_FOUND);
 254   is_one = fc + cat_value_find (var, val);
 255   for (col = fc; col <= lc; col++)
 256     {
 257       entry = (col == is_one) ? 1.0 : 0.0;
 258       gsl_matrix_set (dm->m, row, col, entry);
 259     }
 260 }
 261 void
 262 design_matrix_set_numeric (struct design_matrix *dm, size_t row,
 263                            const struct variable *var, const union value *val)
 264 {
 265   size_t col;
 266
 267   assert (var->type == NUMERIC);
 268   col = design_matrix_var_to_column ((const struct design_matrix *) dm, var);
 269   assert (col != DM_COLUMN_NOT_FOUND);
 270   gsl_matrix_set (dm->m, row, col, val->f);
 271 }