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