1 /* PSPP - Creates design-matrices.
2 Copyright (C) 2005 Free Software Foundation, Inc.
3 Written by Jason H Stover <jason@sakla.net>.
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.
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.
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
21 Create design matrices for procedures that need them.
25 #include <libpspp/alloc.h>
26 #include <libpspp/message.h>
27 #include <data/variable.h>
28 #include <data/category.h>
29 #include "design-matrix.h"
32 #include <gsl/gsl_machine.h>
33 #include <gsl/gsl_vector.h>
34 #include <gsl/gsl_matrix.h>
36 #define DM_COLUMN_NOT_FOUND -1
37 #define DM_INDEX_NOT_FOUND -3
40 Which element of a vector is equal to the value x?
43 cat_which_element_eq (const gsl_vector * vec, double x)
47 for (i = 0; i < vec->size; i++)
49 if (fabs (gsl_vector_get (vec, i) - x) < GSL_DBL_EPSILON)
54 return CAT_VALUE_NOT_FOUND;
57 cat_is_zero_vector (const gsl_vector * vec)
61 for (i = 0; i < vec->size; i++)
63 if (gsl_vector_get (vec, i) != 0.0)
72 Return the value of v corresponding to the vector vec.
75 cat_vector_to_value (const gsl_vector * vec, struct variable *v)
79 i = cat_which_element_eq (vec, 1.0);
80 if (i != CAT_VALUE_NOT_FOUND)
82 return cat_subscript_to_value (i + 1, v);
84 if (cat_is_zero_vector (vec))
86 return cat_subscript_to_value (0, v);
91 struct design_matrix *
92 design_matrix_create (int n_variables,
93 const struct variable *v_variables[],
96 struct design_matrix *dm;
97 const struct variable *v;
102 dm = xmalloc (sizeof *dm);
103 dm->vars = xnmalloc (n_variables, sizeof *dm->vars);
104 dm->n_vars = n_variables;
106 for (i = 0; i < n_variables; i++)
109 assert ((dm->vars + i) != NULL);
110 (dm->vars + i)->v = v; /* Allows us to look up the variable from
111 the design matrix. */
112 (dm->vars + i)->first_column = n_cols;
113 if (v->type == NUMERIC)
115 (dm->vars + i)->last_column = n_cols;
118 else if (v->type == ALPHA)
120 assert (v->obs_vals != NULL);
121 (dm->vars + i)->last_column =
122 (dm->vars + i)->first_column + v->obs_vals->n_categories - 2;
123 n_cols += v->obs_vals->n_categories - 1;
126 dm->m = gsl_matrix_calloc (n_data, n_cols);
133 design_matrix_destroy (struct design_matrix *dm)
136 gsl_matrix_free (dm->m);
141 Return the index of the variable for the
145 design_matrix_col_to_var_index (const struct design_matrix *dm, size_t col)
148 struct design_matrix_var v;
150 for (i = 0; i < dm->n_vars; i++)
153 if (v.first_column <= col && col <= v.last_column)
156 return DM_INDEX_NOT_FOUND;
160 Return a pointer to the variable whose values
161 are stored in column col.
164 design_matrix_col_to_var (const struct design_matrix *dm, size_t col)
168 struct design_matrix_var dmv;
170 index = design_matrix_col_to_var_index (dm, col);
171 for (i = 0; i < dm->n_vars; i++)
174 if ((dmv.v)->index == index)
176 return (struct variable *) dmv.v;
183 cmp_dm_var_index (const struct design_matrix_var *dmv, size_t index)
185 if (dmv->v->index == index)
191 Return the number of the first column which holds the
192 values for variable v.
195 design_matrix_var_to_column (const struct design_matrix * dm,
196 const struct variable * v)
199 struct design_matrix_var tmp;
201 for (i = 0; i < dm->n_vars; i++)
204 if (cmp_dm_var_index (&tmp, v->index))
206 return tmp.first_column;
209 return DM_COLUMN_NOT_FOUND;
214 dm_var_to_last_column (const struct design_matrix *dm,
215 const struct variable *v)
218 struct design_matrix_var tmp;
220 for (i = 0; i < dm->n_vars; i++)
223 if (cmp_dm_var_index (&tmp, v->index))
225 return tmp.last_column;
228 return DM_COLUMN_NOT_FOUND;
232 Set the appropriate value in the design matrix,
233 whether that value is from a categorical or numeric
234 variable. For a categorical variable, only the usual
235 binary encoding is allowed.
238 design_matrix_set_categorical (struct design_matrix *dm, size_t row,
239 const struct variable *var,
240 const union value *val)
248 assert (var->type == ALPHA);
249 fc = design_matrix_var_to_column (dm, var);
250 lc = dm_var_to_last_column (dm, var);
251 assert (lc != DM_COLUMN_NOT_FOUND);
252 assert (fc != DM_COLUMN_NOT_FOUND);
253 is_one = fc + cat_value_find (var, val);
254 for (col = fc; col <= lc; col++)
256 entry = (col == is_one) ? 1.0 : 0.0;
257 gsl_matrix_set (dm->m, row, col, entry);
261 design_matrix_set_numeric (struct design_matrix *dm, size_t row,
262 const struct variable *var, const union value *val)
266 assert (var->type == NUMERIC);
267 col = design_matrix_var_to_column ((const struct design_matrix *) dm, var);
268 assert (col != DM_COLUMN_NOT_FOUND);
269 gsl_matrix_set (dm->m, row, col, val->f);