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 (var_is_numeric (v))
 119         {
 120           (dm->vars + i)->last_column = n_cols;
 121           n_cols++;
 122         }
 123       else if (var_is_alpha (v))
 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 struct variable *
 150 design_matrix_col_to_var (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 (struct variable *) v.v;
 160     }
 161   return NULL;
 162 }
 163
 164 /*
 165   Return the number of the first column which holds the
 166   values for variable v.
 167  */
 168 size_t
 169 design_matrix_var_to_column (const struct design_matrix * dm,
 170                              const struct variable * v)
 171 {
 172   size_t i;
 173   struct design_matrix_var tmp;
 174
 175   for (i = 0; i < dm->n_vars; i++)
 176     {
 177       tmp = dm->vars[i];
 178       if (tmp.v == v)
 179         {
 180           return tmp.first_column;
 181         }
 182     }
 183   return DM_COLUMN_NOT_FOUND;
 184 }
 185
 186 /* Last column. */
 187 static size_t
 188 dm_var_to_last_column (const struct design_matrix *dm,
 189                        const struct variable *v)
 190 {
 191   size_t i;
 192   struct design_matrix_var tmp;
 193
 194   for (i = 0; i < dm->n_vars; i++)
 195     {
 196       tmp = dm->vars[i];
 197       if (tmp.v == v)
 198         {
 199           return tmp.last_column;
 200         }
 201     }
 202   return DM_COLUMN_NOT_FOUND;
 203 }
 204
 205 /*
 206   Set the appropriate value in the design matrix,
 207   whether that value is from a categorical or numeric
 208   variable. For a categorical variable, only the usual
 209   binary encoding is allowed.
 210  */
 211 void
 212 design_matrix_set_categorical (struct design_matrix *dm, size_t row,
 213                                const struct variable *var,
 214                                const union value *val)
 215 {
 216   size_t col;
 217   size_t is_one;
 218   size_t fc;
 219   size_t lc;
 220   double entry;
 221
 222   assert (var_is_alpha (var));
 223   fc = design_matrix_var_to_column (dm, var);
 224   lc = dm_var_to_last_column (dm, var);
 225   assert (lc != DM_COLUMN_NOT_FOUND);
 226   assert (fc != DM_COLUMN_NOT_FOUND);
 227   is_one = fc + cat_value_find (var, val);
 228   for (col = fc; col <= lc; col++)
 229     {
 230       entry = (col == is_one) ? 1.0 : 0.0;
 231       gsl_matrix_set (dm->m, row, col, entry);
 232     }
 233 }
 234 void
 235 design_matrix_set_numeric (struct design_matrix *dm, size_t row,
 236                            const struct variable *var, const union value *val)
 237 {
 238   size_t col;
 239
 240   assert (var_is_numeric (var));
 241   col = design_matrix_var_to_column ((const struct design_matrix *) dm, var);
 242   assert (col != DM_COLUMN_NOT_FOUND);
 243   gsl_matrix_set (dm->m, row, col, val->f);
 244 }