Fix GCC 4.3 warning about uninitialized structure member.
[pspp] / 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_vars = n_variables;
58
59   for (i = 0; i < n_variables; i++)
60     {
61       v = v_variables[i];
62       assert ((dm->vars + i) != NULL);
63       (dm->vars + i)->v = v;    /* Allows us to look up the variable from
64                                    the design matrix. */
65       (dm->vars + i)->first_column = n_cols;
66       if (var_is_numeric (v))
67         {
68           (dm->vars + i)->last_column = n_cols;
69           n_cols++;
70         }
71       else if (var_is_alpha (v))
72         {
73           size_t n_categories = cat_get_n_categories (v);
74           (dm->vars + i)->last_column =
75             (dm->vars + i)->first_column + n_categories - 2;
76           n_cols += n_categories - 1;
77         }
78     }
79   dm->m = gsl_matrix_calloc (n_data, n_cols);
80   col = 0;
81
82   return dm;
83 }
84
85 void
86 design_matrix_destroy (struct design_matrix *dm)
87 {
88   free (dm->vars);
89   gsl_matrix_free (dm->m);
90   free (dm);
91 }
92
93 /*
94   Return the index of the variable for the
95   given column.
96  */
97 const struct variable *
98 design_matrix_col_to_var (const struct design_matrix *dm, size_t col)
99 {
100   size_t i;
101   struct design_matrix_var v;
102
103   for (i = 0; i < dm->n_vars; i++)
104     {
105       v = dm->vars[i];
106       if (v.first_column <= col && col <= v.last_column)
107         return v.v;
108     }
109   return NULL;
110 }
111
112 /*
113   Return the number of the first column which holds the
114   values for variable v.
115  */
116 size_t
117 design_matrix_var_to_column (const struct design_matrix * dm,
118                              const struct variable * v)
119 {
120   size_t i;
121   struct design_matrix_var tmp;
122
123   for (i = 0; i < dm->n_vars; i++)
124     {
125       tmp = dm->vars[i];
126       if (tmp.v == v)
127         {
128           return tmp.first_column;
129         }
130     }
131   return DM_COLUMN_NOT_FOUND;
132 }
133
134 /* Last column. */
135 static size_t
136 dm_var_to_last_column (const struct design_matrix *dm,
137                        const struct variable *v)
138 {
139   size_t i;
140   struct design_matrix_var tmp;
141
142   for (i = 0; i < dm->n_vars; i++)
143     {
144       tmp = dm->vars[i];
145       if (tmp.v == v)
146         {
147           return tmp.last_column;
148         }
149     }
150   return DM_COLUMN_NOT_FOUND;
151 }
152
153 /*
154   Set the appropriate value in the design matrix,
155   whether that value is from a categorical or numeric
156   variable. For a categorical variable, only the usual
157   binary encoding is allowed.
158  */
159 void
160 design_matrix_set_categorical (struct design_matrix *dm, size_t row,
161                                const struct variable *var,
162                                const union value *val)
163 {
164   size_t col;
165   size_t is_one;
166   size_t fc;
167   size_t lc;
168   double entry;
169
170   assert (var_is_alpha (var));
171   fc = design_matrix_var_to_column (dm, var);
172   lc = dm_var_to_last_column (dm, var);
173   assert (lc != DM_COLUMN_NOT_FOUND);
174   assert (fc != DM_COLUMN_NOT_FOUND);
175   is_one = fc + cat_value_find (var, val);
176   for (col = fc; col <= lc; col++)
177     {
178       entry = (col == is_one) ? 1.0 : 0.0;
179       gsl_matrix_set (dm->m, row, col, entry);
180     }
181 }
182
183 void
184 design_matrix_set_numeric (struct design_matrix *dm, size_t row,
185                            const struct variable *var, const union value *val)
186 {
187   size_t col;
188
189   assert (var_is_numeric (var));
190   col = design_matrix_var_to_column ((const struct design_matrix *) dm, var);
191   assert (col != DM_COLUMN_NOT_FOUND);
192   gsl_matrix_set (dm->m, row, col, val->f);
193 }
194
195 struct design_matrix *
196 design_matrix_clone (const struct design_matrix *dm)
197 {
198   struct design_matrix *result;
199   size_t i;
200   
201   assert (dm != NULL);
202   result = xmalloc (sizeof *result);
203   result->vars = xnmalloc (dm->n_vars, sizeof *dm->vars);
204   result->n_vars = dm->n_vars;
205   result->m = gsl_matrix_alloc (dm->m->size1, dm->m->size2);
206   
207   gsl_matrix_memcpy (result->m, dm->m);
208   for (i = 0; i < result->n_vars; i++)
209     {
210       result->vars[i] = dm->vars[i];
211     }
212   return result;
213 }
214