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