Merged changes from branch interaction-review
[pspp-builds.git] / 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 /*
286   Return the subscript of the column of the design matrix
287   corresponding to VAL. If VAR is categorical with d categories, its
288   first category should correspond to the origin in d-dimensional
289   Euclidean space, so there is no subscript for this value.
290  */
291 size_t
292 dm_get_exact_subscript (const struct design_matrix *dm, const struct variable *var,
293                      const union value *val)
294 {
295   size_t result;
296
297   result = design_matrix_var_to_column (dm, var);
298   if (var_is_alpha (var))
299     {
300       if (cat_is_origin (var, val))
301         {
302           return -1u;
303         }
304       result += cat_value_find (var, val) - 1;
305     }
306   return result;
307 }