49 |
#include <gsl/gsl_matrix.h> |
#include <gsl/gsl_matrix.h> |
50 |
|
|
51 |
#define N_INITIAL_CATEGORIES 1 |
#define N_INITIAL_CATEGORIES 1 |
52 |
#define CR_COLUMN_NOT_FOUND -1 |
#define CAT_COLUMN_NOT_FOUND -1 |
53 |
#define CR_VALUE_NOT_FOUND -2 |
#define CAT_VALUE_NOT_FOUND -2 |
54 |
#define CR_INDEX_NOT_FOUND -3 |
#define CAT_INDEX_NOT_FOUND -3 |
55 |
|
|
56 |
static gsl_vector_view cr_value_to_vector (const union value *, |
void |
57 |
struct recoded_categorical *); |
cat_stored_values_create (struct variable *v) |
|
|
|
|
struct recoded_categorical * |
|
|
cr_recoded_categorical_create (const struct variable *v) |
|
58 |
{ |
{ |
59 |
struct recoded_categorical *rc; |
if (v->obs_vals == NULL) |
60 |
|
{ |
61 |
rc = xmalloc (sizeof (*rc)); |
v->obs_vals = xmalloc (sizeof (*v->obs_vals)); |
62 |
rc->v = v; |
v->obs_vals->n_categories = 0; |
63 |
rc->n_categories = 0; |
v->obs_vals->n_allocated_categories = N_INITIAL_CATEGORIES; |
64 |
rc->n_allocated_categories = N_INITIAL_CATEGORIES; |
v->obs_vals->vals = |
65 |
rc->vals = xnmalloc (N_INITIAL_CATEGORIES, sizeof *rc->vals); |
xnmalloc (N_INITIAL_CATEGORIES, sizeof *v->obs_vals->vals); |
66 |
|
} |
|
return rc; |
|
67 |
} |
} |
68 |
|
|
69 |
void |
void |
70 |
cr_recoded_categorical_destroy (struct recoded_categorical *r) |
cat_stored_values_destroy (struct variable *v) |
71 |
{ |
{ |
72 |
free (r->vals); |
assert (v != NULL); |
73 |
free (r); |
if (v->obs_vals != NULL) |
74 |
|
{ |
75 |
|
free (v->obs_vals); |
76 |
|
} |
77 |
} |
} |
78 |
|
|
79 |
|
#if 0 |
80 |
struct recoded_categorical_array * |
struct recoded_categorical_array * |
81 |
cr_recoded_cat_ar_create (int n_variables, struct variable *v_variables[]) |
cr_recoded_cat_ar_create (int n_variables, struct variable *v_variables[]) |
82 |
{ |
{ |
116 |
} |
} |
117 |
return rc; |
return rc; |
118 |
} |
} |
119 |
|
#endif |
120 |
static size_t |
static size_t |
121 |
cr_value_find (struct recoded_categorical *rc, const union value *v) |
cat_value_find (const struct variable *v, const union value *val) |
122 |
{ |
{ |
123 |
size_t i; |
size_t i; |
124 |
const union value *val; |
const union value *candidate; |
125 |
|
|
126 |
for (i = 0; i < rc->n_categories; i++) |
assert (val != NULL); |
127 |
|
assert (v != NULL); |
128 |
|
assert (v->obs_vals != NULL); |
129 |
|
for (i = 0; i < v->obs_vals->n_categories; i++) |
130 |
{ |
{ |
131 |
val = rc->vals + i; |
candidate = v->obs_vals->vals + i; |
132 |
if (!compare_values (val, v, rc->v->width)) |
assert (candidate != NULL); |
133 |
|
if (!compare_values (candidate, val, v->width)) |
134 |
{ |
{ |
135 |
return i; |
return i; |
136 |
} |
} |
137 |
} |
} |
138 |
return CR_VALUE_NOT_FOUND; |
return CAT_VALUE_NOT_FOUND; |
139 |
} |
} |
140 |
|
|
141 |
/* |
/* |
142 |
Add the new value unless it is already present. |
Add the new value unless it is already present. |
143 |
*/ |
*/ |
144 |
void |
void |
145 |
cr_value_update (struct recoded_categorical *rc, const union value *v) |
cat_value_update (struct variable *v, const union value *val) |
146 |
{ |
{ |
147 |
if (cr_value_find (rc, v) == CR_VALUE_NOT_FOUND) |
struct cat_vals *cv; |
148 |
|
|
149 |
|
if (v->type == ALPHA) |
150 |
{ |
{ |
151 |
if (rc->n_categories >= rc->n_allocated_categories) |
assert (val != NULL); |
152 |
{ |
assert (v != NULL); |
153 |
rc->n_allocated_categories *= 2; |
cv = v->obs_vals; |
154 |
rc->vals = xnrealloc (rc->vals, |
if (cat_value_find (v, val) == CAT_VALUE_NOT_FOUND) |
155 |
rc->n_allocated_categories, sizeof *rc->vals); |
{ |
156 |
|
if (cv->n_categories >= cv->n_allocated_categories) |
157 |
|
{ |
158 |
|
cv->n_allocated_categories *= 2; |
159 |
|
cv->vals = xnrealloc (cv->vals, |
160 |
|
cv->n_allocated_categories, |
161 |
|
sizeof *cv->vals); |
162 |
|
} |
163 |
|
cv->vals[cv->n_categories] = *val; |
164 |
|
cv->n_categories++; |
165 |
} |
} |
|
rc->vals[rc->n_categories] = *v; |
|
|
rc->n_categories++; |
|
166 |
} |
} |
167 |
} |
} |
168 |
|
|
169 |
/* |
/* |
170 |
Create a set of gsl_matrix's, each of whose rows correspond to |
Create a gsl_matrix, whose rows correspond to values of a |
171 |
values of a categorical variable. Since n categories have n-1 |
categorical variable. Since n categories have n-1 degrees of |
172 |
degrees of freedom, the gsl_matrix is n-by-(n-1), with the first |
freedom, the gsl_matrix is n-by-(n-1), with the first category |
173 |
category encoded as the zero vector. |
encoded as the zero vector. |
174 |
*/ |
*/ |
175 |
|
#if 0 |
176 |
void |
void |
177 |
cr_create_value_matrices (struct recoded_categorical_array *r) |
cat_create_value_matrix (struct variable *v) |
178 |
{ |
{ |
179 |
size_t i; |
size_t i; |
180 |
size_t row; |
size_t row; |
182 |
size_t n_rows; |
size_t n_rows; |
183 |
size_t n_cols; |
size_t n_cols; |
184 |
|
|
185 |
for (i = 0; i < r->n_vars; i++) |
assert (v != NULL); |
186 |
|
if (v->type == ALPHA) |
187 |
{ |
{ |
188 |
n_rows = (*(r->a + i))->n_categories; |
assert (v->rc != NULL); |
189 |
n_cols = (*(r->a + i))->n_categories - 1; |
n_rows = v->rc->n_categories; |
190 |
(*(r->a + i))->m = gsl_matrix_calloc (n_rows, n_cols); |
n_cols = v->rc->n_categories - 1; |
191 |
|
v->rc->m = gsl_matrix_calloc (n_rows, n_cols); |
192 |
for (row = 1; row < n_rows; row++) |
for (row = 1; row < n_rows; row++) |
193 |
{ |
{ |
194 |
col = row - 1; |
col = row - 1; |
195 |
gsl_matrix_set ((*(r->a + i))->m, row, col, 1.0); |
gsl_matrix_set (v->rc->m, row, col, 1.0); |
196 |
} |
} |
197 |
} |
} |
198 |
} |
} |
199 |
|
#endif |
200 |
static size_t |
/* |
201 |
cr_value_to_subscript (const union value *val, struct recoded_categorical *cr) |
Return the subscript of the binary vector corresponding |
202 |
|
to this value. |
203 |
|
*/ |
204 |
|
size_t |
205 |
|
cat_value_to_subscript (const union value *val, struct variable *v) |
206 |
{ |
{ |
207 |
const union value *v; |
const union value *val2; |
208 |
size_t subscript; |
size_t subscript; |
209 |
int different; |
int different; |
210 |
|
|
211 |
subscript = cr->n_categories - 1; |
assert (v != NULL); |
212 |
|
assert (val != NULL); |
213 |
|
assert (v->obs_vals != NULL); |
214 |
|
subscript = v->obs_vals->n_categories - 1; |
215 |
while (subscript > 0) |
while (subscript > 0) |
216 |
{ |
{ |
217 |
v = cr->vals + subscript; |
val2 = v->obs_vals->vals + subscript; |
218 |
different = compare_values (val, v, cr->v->width); |
assert (val2 != NULL); |
219 |
|
different = compare_values (val, val2, v->width); |
220 |
if (!different) |
if (!different) |
221 |
{ |
{ |
222 |
return subscript; |
return subscript; |
226 |
return subscript; |
return subscript; |
227 |
} |
} |
228 |
|
|
229 |
static union value * |
union value * |
230 |
cr_subscript_to_value (const size_t s, struct recoded_categorical *cr) |
cat_subscript_to_value (const size_t s, struct variable *v) |
231 |
{ |
{ |
232 |
if (s < cr->n_categories) |
assert (v->obs_vals != NULL); |
233 |
|
if (s < v->obs_vals->n_categories) |
234 |
{ |
{ |
235 |
return (cr->vals + s); |
return (v->obs_vals->vals + s); |
236 |
} |
} |
237 |
else |
else |
238 |
{ |
{ |
240 |
} |
} |
241 |
} |
} |
242 |
|
|
243 |
|
#if 0 |
244 |
/* |
/* |
245 |
Return the row of the matrix corresponding |
Return the row of the matrix corresponding |
246 |
to the value v. |
to the value v. |
252 |
row = cr_value_to_subscript (v, cr); |
row = cr_value_to_subscript (v, cr); |
253 |
return gsl_matrix_row (cr->m, row); |
return gsl_matrix_row (cr->m, row); |
254 |
} |
} |
255 |
|
#endif |
256 |
/* |
/* |
257 |
Which element of a vector is equal to the value x? |
Which element of a vector is equal to the value x? |
258 |
*/ |
*/ |
259 |
static size_t |
static size_t |
260 |
cr_which_element_eq (const gsl_vector * vec, double x) |
cat_which_element_eq (const gsl_vector * vec, double x) |
261 |
{ |
{ |
262 |
size_t i; |
size_t i; |
263 |
|
|
268 |
return i; |
return i; |
269 |
} |
} |
270 |
} |
} |
271 |
return CR_VALUE_NOT_FOUND; |
return CAT_VALUE_NOT_FOUND; |
272 |
} |
} |
273 |
static int |
static int |
274 |
cr_is_zero_vector (const gsl_vector * vec) |
cat_is_zero_vector (const gsl_vector * vec) |
275 |
{ |
{ |
276 |
size_t i; |
size_t i; |
277 |
|
|
293 |
i is 0 otherwise. |
i is 0 otherwise. |
294 |
*/ |
*/ |
295 |
union value * |
union value * |
296 |
cr_vector_to_value (const gsl_vector * vec, struct recoded_categorical *cr) |
cr_vector_to_value (const gsl_vector * vec, struct variable *v) |
297 |
{ |
{ |
298 |
size_t i; |
size_t i; |
299 |
|
|
300 |
i = cr_which_element_eq (vec, 1.0); |
i = cat_which_element_eq (vec, 1.0); |
301 |
if (i != CR_VALUE_NOT_FOUND) |
if (i != CAT_VALUE_NOT_FOUND) |
302 |
{ |
{ |
303 |
return cr_subscript_to_value (i + 1, cr); |
return cat_subscript_to_value (i + 1, v); |
304 |
} |
} |
305 |
if (cr_is_zero_vector (vec)) |
if (cat_is_zero_vector (vec)) |
306 |
{ |
{ |
307 |
return cr_subscript_to_value (0, cr); |
return cat_subscript_to_value (0, v); |
308 |
} |
} |
309 |
return NULL; |
return NULL; |
310 |
} |
} |
311 |
|
|
312 |
|
#if 0 |
313 |
/* |
/* |
314 |
Given a variable, return a pointer to its recoded |
Given a variable, return a pointer to its recoded |
315 |
structure. BUSTED IN HERE. |
structure. BUSTED IN HERE. |
331 |
} |
} |
332 |
return NULL; |
return NULL; |
333 |
} |
} |
334 |
|
#endif |
335 |
struct design_matrix * |
struct design_matrix * |
336 |
design_matrix_create (int n_variables, |
design_matrix_create (int n_variables, |
337 |
const struct variable *v_variables[], |
const struct variable *v_variables[], |
|
struct recoded_categorical_array *ca, |
|
338 |
const size_t n_data) |
const size_t n_data) |
339 |
{ |
{ |
340 |
struct design_matrix *dm; |
struct design_matrix *dm; |
|
struct design_matrix_var *tmp; |
|
|
struct recoded_categorical *rc; |
|
341 |
const struct variable *v; |
const struct variable *v; |
342 |
size_t i; |
size_t i; |
343 |
size_t n_cols = 0; |
size_t n_cols = 0; |
350 |
for (i = 0; i < n_variables; i++) |
for (i = 0; i < n_variables; i++) |
351 |
{ |
{ |
352 |
v = v_variables[i]; |
v = v_variables[i]; |
353 |
|
assert ((dm->vars + i) != NULL); |
354 |
|
(dm->vars + i)->v = v; /* Allows us to look up the variable from |
355 |
|
the design matrix. */ |
356 |
|
(dm->vars + i)->first_column = n_cols; |
357 |
if (v->type == NUMERIC) |
if (v->type == NUMERIC) |
358 |
{ |
{ |
359 |
n_cols++; |
n_cols++; |
360 |
|
(dm->vars + i)->last_column = n_cols; |
361 |
} |
} |
362 |
else if (v->type == ALPHA) |
else if (v->type == ALPHA) |
363 |
{ |
{ |
364 |
assert (ca != NULL); |
assert (v->obs_vals != NULL); |
365 |
rc = cr_var_to_recoded_categorical (v, ca); |
(dm->vars + i)->last_column = |
366 |
assert (rc != NULL); |
(dm->vars + i)->first_column + v->obs_vals->n_categories - 2; |
367 |
rc->first_column = n_cols; |
n_cols += v->obs_vals->n_categories - 1; |
|
rc->last_column = rc->first_column + rc->n_categories - 2; |
|
|
n_cols += rc->n_categories - 1; |
|
368 |
} |
} |
369 |
} |
} |
370 |
dm->m = gsl_matrix_calloc (n_data, n_cols); |
dm->m = gsl_matrix_calloc (n_data, n_cols); |
|
dm->vars = xnmalloc (dm->n_vars, sizeof *dm->vars); |
|
|
assert (dm->vars != NULL); |
|
371 |
col = 0; |
col = 0; |
372 |
|
|
|
for (i = 0; i < n_variables; i++) |
|
|
{ |
|
|
v = v_variables[i]; |
|
|
(dm->vars[i]).v = v; |
|
|
if (v->type == NUMERIC) |
|
|
{ |
|
|
tmp = &(dm->vars[col]); |
|
|
tmp->v = v; |
|
|
tmp->first_column = col; |
|
|
tmp->last_column = col; |
|
|
col++; |
|
|
} |
|
|
else if (v->type == ALPHA) |
|
|
{ |
|
|
assert (ca != NULL); |
|
|
rc = cr_var_to_recoded_categorical (v, ca); |
|
|
assert (rc != NULL); |
|
|
tmp = &(dm->vars[col]); |
|
|
tmp->v = v; |
|
|
tmp->last_column = rc->last_column; |
|
|
col = rc->last_column + 1; |
|
|
} |
|
|
} |
|
373 |
return dm; |
return dm; |
374 |
} |
} |
375 |
|
|
397 |
if (v.first_column <= col && col <= v.last_column) |
if (v.first_column <= col && col <= v.last_column) |
398 |
return (v.v)->index; |
return (v.v)->index; |
399 |
} |
} |
400 |
return CR_INDEX_NOT_FOUND; |
return CAT_INDEX_NOT_FOUND; |
401 |
} |
} |
402 |
|
|
403 |
/* |
/* |
432 |
} |
} |
433 |
|
|
434 |
/* |
/* |
435 |
Return the number of the first column storing the |
Return the number of the first column which holds the |
436 |
values for variable v. |
values for variable v. |
437 |
*/ |
*/ |
438 |
size_t |
size_t |
450 |
return tmp.first_column; |
return tmp.first_column; |
451 |
} |
} |
452 |
} |
} |
453 |
return CR_COLUMN_NOT_FOUND; |
return CAT_COLUMN_NOT_FOUND; |
454 |
|
} |
455 |
|
|
456 |
|
/* Last column. */ |
457 |
|
static size_t |
458 |
|
dm_var_to_last_column (const struct design_matrix *dm, |
459 |
|
const struct variable *v) |
460 |
|
{ |
461 |
|
size_t i; |
462 |
|
struct design_matrix_var tmp; |
463 |
|
|
464 |
|
for (i = 0; i < dm->n_vars; i++) |
465 |
|
{ |
466 |
|
tmp = dm->vars[i]; |
467 |
|
if (cmp_dm_var_index (&tmp, v->index)) |
468 |
|
{ |
469 |
|
return tmp.last_column; |
470 |
|
} |
471 |
|
} |
472 |
|
return CAT_COLUMN_NOT_FOUND; |
473 |
} |
} |
474 |
|
|
475 |
/* |
/* |
476 |
Set the appropriate value in the design matrix, |
Set the appropriate value in the design matrix, |
477 |
whether that value is from a categorical or numeric |
whether that value is from a categorical or numeric |
478 |
variable. |
variable. For a categorical variable, only the usual |
479 |
|
binary encoding is allowed. |
480 |
*/ |
*/ |
481 |
void |
void |
482 |
design_matrix_set_categorical (struct design_matrix *dm, size_t row, |
design_matrix_set_categorical (struct design_matrix *dm, size_t row, |
483 |
const struct variable *var, |
const struct variable *var, |
484 |
const union value *val, |
const union value *val) |
|
struct recoded_categorical *rc) |
|
485 |
{ |
{ |
486 |
size_t col; |
size_t col; |
487 |
double x; |
size_t is_one; |
488 |
|
size_t fc; |
489 |
|
size_t lc; |
490 |
|
double entry; |
491 |
|
|
492 |
assert (var->type == ALPHA); |
assert (var->type == ALPHA); |
493 |
const gsl_vector_view vec = cr_value_to_vector (val, rc); |
fc = design_matrix_var_to_column (dm, var); |
494 |
|
lc = dm_var_to_last_column (dm, var); |
495 |
/* |
assert (lc != CAT_COLUMN_NOT_FOUND); |
496 |
Copying values here is not the 'most efficient' way, |
assert (fc != CAT_COLUMN_NOT_FOUND); |
497 |
but it will work even if we change the vector encoding later. |
is_one = fc + cat_value_find (var, val); |
498 |
*/ |
for (col = fc; col <= lc; col++) |
|
for (col = rc->first_column; col <= rc->last_column; col++) |
|
499 |
{ |
{ |
500 |
x = gsl_vector_get (&vec.vector, col); |
entry = (col == is_one) ? 1.0 : 0.0; |
501 |
gsl_matrix_set (dm->m, row, col, x); |
gsl_matrix_set (dm->m, row, col, entry); |
502 |
} |
} |
503 |
} |
} |
504 |
void |
void |
509 |
|
|
510 |
assert (var->type == NUMERIC); |
assert (var->type == NUMERIC); |
511 |
col = design_matrix_var_to_column ((const struct design_matrix *) dm, var); |
col = design_matrix_var_to_column ((const struct design_matrix *) dm, var); |
512 |
assert (col != CR_COLUMN_NOT_FOUND); |
assert (col != CAT_COLUMN_NOT_FOUND); |
513 |
gsl_matrix_set (dm->m, row, col, val->f); |
gsl_matrix_set (dm->m, row, col, val->f); |
514 |
} |
} |