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 struct design_matrix *
  47 design_matrix_create (int n_variables,
  48                       const struct variable *v_variables[],
  49                       const size_t n_data)
  50 {
  51   struct design_matrix *dm;
  52   const struct variable *v;
  53   size_t i;
  54   size_t n_cols = 0;
  55   size_t col;
  56
  57   dm = xmalloc (sizeof *dm);
  58   dm->vars = xnmalloc (n_variables, sizeof *dm->vars);
  59   dm->n_vars = n_variables;
  60
  61   for (i = 0; i < n_variables; i++)
  62     {
  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   return dm;
  85 }
  86
  87 void
  88 design_matrix_destroy (struct design_matrix *dm)
  89 {
  90   free (dm->vars);
  91   gsl_matrix_free (dm->m);
  92   free (dm);
  93 }
  94
  95 /*
  96   Return the index of the variable for the
  97   given column.
  98  */
  99 const struct variable *
 100 design_matrix_col_to_var (const struct design_matrix *dm, size_t col)
 101 {
 102   size_t i;
 103   struct design_matrix_var v;
 104
 105   for (i = 0; i < dm->n_vars; i++)
 106     {
 107       v = dm->vars[i];
 108       if (v.first_column <= col && col <= v.last_column)
 109         return v.v;
 110     }
 111   return NULL;
 112 }
 113
 114 /*
 115   Return the number of the first column which holds the
 116   values for variable v.
 117  */
 118 size_t
 119 design_matrix_var_to_column (const struct design_matrix * dm,
 120                              const struct variable * v)
 121 {
 122   size_t i;
 123   struct design_matrix_var tmp;
 124
 125   for (i = 0; i < dm->n_vars; i++)
 126     {
 127       tmp = dm->vars[i];
 128       if (tmp.v == v)
 129         {
 130           return tmp.first_column;
 131         }
 132     }
 133   return DM_COLUMN_NOT_FOUND;
 134 }
 135
 136 /* Last column. */
 137 static size_t
 138 dm_var_to_last_column (const struct design_matrix *dm,
 139                        const struct variable *v)
 140 {
 141   size_t i;
 142   struct design_matrix_var tmp;
 143
 144   for (i = 0; i < dm->n_vars; i++)
 145     {
 146       tmp = dm->vars[i];
 147       if (tmp.v == v)
 148         {
 149           return tmp.last_column;
 150         }
 151     }
 152   return DM_COLUMN_NOT_FOUND;
 153 }
 154
 155 /*
 156   Set the appropriate value in the design matrix,
 157   whether that value is from a categorical or numeric
 158   variable. For a categorical variable, only the usual
 159   binary encoding is allowed.
 160  */
 161 void
 162 design_matrix_set_categorical (struct design_matrix *dm, size_t row,
 163                                const struct variable *var,
 164                                const union value *val)
 165 {
 166   size_t col;
 167   size_t is_one;
 168   size_t fc;
 169   size_t lc;
 170   double entry;
 171
 172   assert (var_is_alpha (var));
 173   fc = design_matrix_var_to_column (dm, var);
 174   lc = dm_var_to_last_column (dm, var);
 175   assert (lc != DM_COLUMN_NOT_FOUND);
 176   assert (fc != DM_COLUMN_NOT_FOUND);
 177   is_one = fc + cat_value_find (var, val);
 178   for (col = fc; col <= lc; col++)
 179     {
 180       entry = (col == is_one) ? 1.0 : 0.0;
 181       gsl_matrix_set (dm->m, row, col, entry);
 182     }
 183 }
 184
 185 void
 186 design_matrix_set_numeric (struct design_matrix *dm, size_t row,
 187                            const struct variable *var, const union value *val)
 188 {
 189   size_t col;
 190
 191   assert (var_is_numeric (var));
 192   col = design_matrix_var_to_column ((const struct design_matrix *) dm, var);
 193   assert (col != DM_COLUMN_NOT_FOUND);
 194   gsl_matrix_set (dm->m, row, col, val->f);
 195 }