src/math/design-matrix.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2005 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 /*
  18   Create design matrices for procedures that need them.
  19 */
  20 #include <config.h>
  21
  22 #include "design-matrix.h"
  23
  24 #include <assert.h>
  25 #include <math.h>
  26 #include <stdlib.h>
  27 #include <string.h>
  28
  29 #include <libpspp/message.h>
  30 #include <data/variable.h>
  31 #include <data/category.h>
  32 #include <data/value.h>
  33
  34 #include <gsl/gsl_machine.h>
  35 #include <gsl/gsl_vector.h>
  36 #include <gsl/gsl_matrix.h>
  37
  38 #include "xalloc.h"
  39
  40 #define DM_COLUMN_NOT_FOUND -1
  41 #define DM_INDEX_NOT_FOUND -3
  42
  43
  44 struct design_matrix *
  45 design_matrix_create (int n_variables,
  46                       const struct variable *v_variables[],
  47                       const size_t n_data)
  48 {
  49   struct design_matrix *dm;
  50   const struct variable *v;
  51   size_t i;
  52   size_t n_cols = 0;
  53   size_t col;
  54
  55   dm = xmalloc (sizeof *dm);
  56   dm->vars = xnmalloc (n_variables, sizeof *dm->vars);
  57   dm->n_cases = xnmalloc (n_variables, sizeof (*dm->n_cases));
  58   dm->n_vars = n_variables;
  59
  60   for (i = 0; i < n_variables; i++)
  61     {
  62       dm->n_cases[i] = 0;
  63       v = v_variables[i];
  64       assert ((dm->vars + i) != NULL);
  65       (dm->vars + i)->v = v;    /* Allows us to look up the variable from
  66                                    the design matrix. */
  67       (dm->vars + i)->first_column = n_cols;
  68       if (var_is_numeric (v))
  69         {
  70           (dm->vars + i)->last_column = n_cols;
  71           n_cols++;
  72         }
  73       else if (var_is_alpha (v))
  74         {
  75           size_t n_categories = cat_get_n_categories (v);
  76           (dm->vars + i)->last_column =
  77             (dm->vars + i)->first_column + n_categories - 2;
  78           n_cols += n_categories - 1;
  79         }
  80     }
  81   dm->m = gsl_matrix_calloc (n_data, n_cols);
  82   col = 0;
  83
  84
  85   return dm;
  86 }
  87
  88 void
  89 design_matrix_destroy (struct design_matrix *dm)
  90 {
  91   free (dm->vars);
  92   gsl_matrix_free (dm->m);
  93   free (dm->n_cases);
  94   free (dm);
  95 }
  96
  97 /*
  98   Return the index of the variable for the
  99   given column.
 100  */
 101 const struct variable *
 102 design_matrix_col_to_var (const struct design_matrix *dm, size_t col)
 103 {
 104   size_t i;
 105   struct design_matrix_var v;
 106
 107   for (i = 0; i < dm->n_vars; i++)
 108     {
 109       v = dm->vars[i];
 110       if (v.first_column <= col && col <= v.last_column)
 111         return v.v;
 112     }
 113   return NULL;
 114 }
 115
 116 /*
 117   Return the number of the first column which holds the
 118   values for variable v.
 119  */
 120 size_t
 121 design_matrix_var_to_column (const struct design_matrix * dm,
 122                              const struct variable * v)
 123 {
 124   size_t i;
 125   struct design_matrix_var tmp;
 126
 127   for (i = 0; i < dm->n_vars; i++)
 128     {
 129       tmp = dm->vars[i];
 130       if (tmp.v == v)
 131         {
 132           return tmp.first_column;
 133         }
 134     }
 135   return DM_COLUMN_NOT_FOUND;
 136 }
 137
 138 /* Last column. */
 139 static size_t
 140 dm_var_to_last_column (const struct design_matrix *dm,
 141                        const struct variable *v)
 142 {
 143   size_t i;
 144   struct design_matrix_var tmp;
 145
 146   for (i = 0; i < dm->n_vars; i++)
 147     {
 148       tmp = dm->vars[i];
 149       if (tmp.v == v)
 150         {
 151           return tmp.last_column;
 152         }
 153     }
 154   return DM_COLUMN_NOT_FOUND;
 155 }
 156
 157 /*
 158   Set the appropriate value in the design matrix,
 159   whether that value is from a categorical or numeric
 160   variable. For a categorical variable, only the usual
 161   binary encoding is allowed.
 162  */
 163 void
 164 design_matrix_set_categorical (struct design_matrix *dm, size_t row,
 165                                const struct variable *var,
 166                                const union value *val)
 167 {
 168   size_t col;
 169   size_t is_one;
 170   size_t fc;
 171   size_t lc;
 172   double entry;
 173
 174   assert (var_is_alpha (var));
 175   fc = design_matrix_var_to_column (dm, var);
 176   lc = dm_var_to_last_column (dm, var);
 177   assert (lc != DM_COLUMN_NOT_FOUND);
 178   assert (fc != DM_COLUMN_NOT_FOUND);
 179   is_one = fc + cat_value_find (var, val);
 180   for (col = fc; col <= lc; col++)
 181     {
 182       entry = (col == is_one) ? 1.0 : 0.0;
 183       gsl_matrix_set (dm->m, row, col, entry);
 184     }
 185 }
 186
 187 void
 188 design_matrix_set_numeric (struct design_matrix *dm, size_t row,
 189                            const struct variable *var, const union value *val)
 190 {
 191   size_t col;
 192
 193   assert (var_is_numeric (var));
 194   col = design_matrix_var_to_column ((const struct design_matrix *) dm, var);
 195   assert (col != DM_COLUMN_NOT_FOUND);
 196   gsl_matrix_set (dm->m, row, col, val->f);
 197 }
 198
 199 struct design_matrix *
 200 design_matrix_clone (const struct design_matrix *dm)
 201 {
 202   struct design_matrix *result;
 203   size_t i;
 204
 205   assert (dm != NULL);
 206   result = xmalloc (sizeof *result);
 207   result->vars = xnmalloc (dm->n_vars, sizeof *dm->vars);
 208   result->n_vars = dm->n_vars;
 209   result->m = gsl_matrix_alloc (dm->m->size1, dm->m->size2);
 210
 211   gsl_matrix_memcpy (result->m, dm->m);
 212   for (i = 0; i < result->n_vars; i++)
 213     {
 214       result->vars[i] = dm->vars[i];
 215     }
 216   return result;
 217 }
 218
 219 /*
 220   Increment the number of cases for V.
 221  */
 222 void
 223 design_matrix_increment_case_count (struct design_matrix *dm, const struct variable *v)
 224 {
 225   size_t i;
 226   assert (dm != NULL);
 227   assert (dm->n_cases != NULL);
 228   assert (v != NULL);
 229   i = design_matrix_var_to_column (dm, v);
 230   dm->n_cases[i]++;
 231 }
 232
 233 /*
 234   Set the number of cases for V.
 235  */
 236 void
 237 design_matrix_set_case_count (struct design_matrix *dm, const struct variable *v, size_t n)
 238 {
 239   size_t i;
 240   assert (dm != NULL);
 241   assert (dm->n_cases != NULL);
 242   assert (v != NULL);
 243   i = design_matrix_var_to_column (dm, v);
 244   dm->n_cases[i] = n;
 245 }
 246
 247 /*
 248   Get the number of cases for V.
 249  */
 250 size_t
 251 design_matrix_get_case_count (const struct design_matrix *dm, const struct variable *v)
 252 {
 253   size_t i;
 254   assert (dm != NULL);
 255   assert (dm->n_cases != NULL);
 256   assert (v != NULL);
 257   i = design_matrix_var_to_column (dm, v);
 258   return dm->n_cases[i];
 259 }
 260
 261 size_t
 262 design_matrix_get_n_cols (const struct design_matrix *d)
 263 {
 264   return d->m->size2;
 265 }
 266
 267 size_t
 268 design_matrix_get_n_rows (const struct design_matrix *d)
 269 {
 270   return d->m->size1;
 271 }
 272
 273 double
 274 design_matrix_get_element (const struct design_matrix *d, size_t row, size_t col)
 275 {
 276   return (gsl_matrix_get (d->m, row, col));
 277 }
 278
 279 void
 280 design_matrix_set_element (const struct design_matrix *d, size_t row, size_t col, double x)
 281 {
 282   gsl_matrix_set (d->m, row, col, x);
 283 }
 284 /*
 285   If VAR is categorical with d categories, its first category should
 286   correspond to the origin in d-dimensional Euclidean space.
 287  */
 288 static bool
 289 is_origin (const struct variable *var, const union value *val)
 290 {
 291   if (var_is_numeric (var))
 292     {
 293       return false;
 294     }
 295   if (cat_value_find (var, val) == 0)
 296     {
 297       return true;
 298     }
 299   return false;
 300 }
 301
 302 /*
 303   Return the subscript of the column of the design matrix
 304   corresponding to VAL. If VAR is categorical with d categories, its
 305   first category should correspond to the origin in d-dimensional
 306   Euclidean space, so there is no subscript for this value.
 307  */
 308 size_t
 309 dm_get_exact_subscript (const struct design_matrix *dm, const struct variable *var,
 310                      const union value *val)
 311 {
 312   size_t result;
 313
 314   result = design_matrix_var_to_column (dm, var);
 315   if (var_is_alpha (var))
 316     {
 317       if (is_origin (var, val))
 318         {
 319           return -1u;
 320         }
 321       result += cat_value_find (var, val) - 1;
 322     }
 323   return result;
 324 }