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