1 /*
2 * gretl -- Gnu Regression, Econometrics and Time-series Library
3 * Copyright (C) 2001 Allin Cottrell and Riccardo "Jack" Lucchetti
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU 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, see <http://www.gnu.org/licenses/>.
17 *
18 */
19
20 #include "libgretl.h"
21 #include "usermat.h"
22 #include "libset.h"
23 #include "gretl_cmatrix.h"
24 #include "matrix_extra.h"
25 #include "swap_bytes.h"
26
27 #ifdef WIN32
28 # include "gretl_win32.h"
29 #endif
30
31 #include <errno.h>
32
33 /**
34 * SECTION:matrix_extra
35 * @short_description: more matrix operations
36 * @title: Matrix extra
37 * @include: gretl/libgretl.h, gretl/matrix_extra.h
38 *
39 * Functions that convert between gretl_matrix and other
40 * datatypes; also printing of gretl_matrices and some
41 * functions relating to masking and cutting rows and/or
42 * columns of matrices.
43 */
44
45 /**
46 * gretl_vector_from_array:
47 * @x: pointer to array of elements.
48 * @n: number of elements.
49 * @mod: modifier flag: either %GRETL_MOD_NONE, or %GRETL_MOD_SQUARE
50 * to use the squares of the elements of @x.
51 *
52 * Returns: pointer to a newly allocated gretl_vector containing
53 * the elements of x (or their squares), or NULL on failure.
54 * Missing values in @x are skipped.
55 */
56
57 gretl_vector *
gretl_vector_from_array(const double * x,int n,GretlMatrixMod mod)58 gretl_vector_from_array (const double *x, int n, GretlMatrixMod mod)
59 {
60 gretl_matrix *v;
61 double xi;
62
63 v = gretl_column_vector_alloc(n);
64
65 if (v != NULL) {
66 int i = 0, j = 0;
67
68 while (j < n) {
69 xi = x[i++];
70 if (!na(xi)) {
71 if (mod == GRETL_MOD_SQUARE) {
72 v->val[j] = xi * xi;
73 } else {
74 v->val[j] = xi;
75 }
76 j++;
77 }
78 }
79 }
80
81 return v;
82 }
83
84 /**
85 * gretl_vector_from_series:
86 * @x: series from data array.
87 * @t1: starting observation.
88 * @t2: ending observation.
89 *
90 * Returns: a newly allocated gretl_vector containing the values
91 * of the given data series for the given range, or NULL on failure.
92 * Any missing values in the input array are preserved as NaNs in
93 * the output.
94 */
95
gretl_vector_from_series(const double * x,int t1,int t2)96 gretl_vector *gretl_vector_from_series (const double *x,
97 int t1, int t2)
98 {
99 gretl_matrix *v = NULL;
100 int n = t2 - t1 + 1;
101
102 if (n > 0) {
103 size_t sz = n * sizeof *x;
104
105 v = gretl_column_vector_alloc(n);
106 if (v != NULL) {
107 memcpy(v->val, x + t1, sz);
108 gretl_matrix_set_t1(v, t1);
109 gretl_matrix_set_t2(v, t2);
110 }
111 }
112
113 return v;
114 }
115
116 /**
117 * gretl_matrix_from_2d_array:
118 * @X: two-dimensional array of doubles.
119 * @rows: number of rows in target matrix.
120 * @cols: number of columns in target matrix.
121 *
122 * Returns: allocated gretl_matrix, the elements of which are set to
123 * the values in @X, or NULL on allocation failure.
124 */
125
gretl_matrix_from_2d_array(const double ** X,int rows,int cols)126 gretl_matrix *gretl_matrix_from_2d_array (const double **X,
127 int rows, int cols)
128 {
129 gretl_matrix *m;
130 int i, j, k = 0;
131
132 m = gretl_matrix_alloc(rows, cols);
133
134 if (m != NULL) {
135 for (j=0; j<cols; j++) {
136 for (i=0; i<rows; i++) {
137 m->val[k++] = X[j][i];
138 }
139 }
140 }
141
142 return m;
143 }
144
145 /**
146 * gretl_matrix_from_scalar:
147 * @x: scalar to be "promoted".
148 *
149 * Returns: allocated 1 x 1 gretl_matrix, the single element
150 * of which is set to @x, or NULL on allocation failure.
151 */
152
gretl_matrix_from_scalar(double x)153 gretl_matrix *gretl_matrix_from_scalar (double x)
154 {
155 gretl_matrix *m = gretl_matrix_alloc(1, 1);
156
157 m->val[0] = x;
158
159 return m;
160 }
161
count_selection(const char * s,int n)162 static int count_selection (const char *s, int n)
163 {
164 int i, c = 0;
165
166 for (i=0; i<n; i++) {
167 if (s[i] != 0) c++;
168 }
169
170 return c;
171 }
172
173 /**
174 * gretl_vcv_matrix_from_model:
175 * @pmod: pointer to model
176 * @select: char array indicating which rows and colums to select
177 * (or NULL for the full matrix).
178 * @err: location to receive error code.
179 *
180 * Produces all or part of the covariance matrix for @pmod
181 * in the form of a gretl_matrix. Storage is allocated, to be freed
182 * by the caller. If @select is not NULL, it should be an array
183 * with non-zero elements in positions corresponding to the
184 * desired rows (and columns), and zero elements otherwise.
185 *
186 * Returns: the covariance matrix, or NULL on error.
187 */
188
189 gretl_matrix *
gretl_vcv_matrix_from_model(MODEL * pmod,const char * select,int * err)190 gretl_vcv_matrix_from_model (MODEL *pmod, const char *select, int *err)
191 {
192 gretl_matrix *vcv;
193 int i, j, idx, nc;
194 int ii, jj;
195 int k = pmod->ncoeff;
196
197 /* first ensure the model _has_ a vcv */
198 *err = makevcv(pmod, pmod->sigma);
199 if (*err) {
200 return NULL;
201 }
202
203 if (select == NULL) {
204 nc = k;
205 } else {
206 nc = count_selection(select, k);
207 }
208
209 if (nc == 0) {
210 *err = E_DATA;
211 return NULL;
212 }
213
214 vcv = gretl_matrix_alloc(nc, nc);
215 if (vcv == NULL) {
216 *err = E_ALLOC;
217 return NULL;
218 }
219
220 ii = 0;
221 for (i=0; i<k; i++) {
222 if (select != NULL && !select[i]) {
223 continue;
224 }
225 jj = 0;
226 for (j=0; j<=i; j++) {
227 if (select != NULL && !select[j]) {
228 continue;
229 }
230 idx = ijton(i, j, pmod->ncoeff);
231 gretl_matrix_set(vcv, ii, jj, pmod->vcv[idx]);
232 if (jj != ii) {
233 gretl_matrix_set(vcv, jj, ii, pmod->vcv[idx]);
234 }
235 jj++;
236 }
237 ii++;
238 }
239
240 return vcv;
241 }
242
243 /**
244 * gretl_coeff_vector_from_model:
245 * @pmod: pointer to model
246 * @select: char array indicating which rows to select
247 * (or NULL for the full vector).
248 * @err: location to receive error code.
249 *
250 * Produces all or part of the coefficient vector for @pmod
251 * in the form of a gretl column vector. Storage is allocated, to be freed
252 * by the caller. If @select is non-NULL, it should be an array
253 * with non-zero elements in positions corresponding to the
254 * desired rows and zero elements otherwise.
255 *
256 * Returns: the coefficient vector, or NULL on error.
257 */
258
259 gretl_vector *
gretl_coeff_vector_from_model(const MODEL * pmod,const char * select,int * err)260 gretl_coeff_vector_from_model (const MODEL *pmod, const char *select, int *err)
261 {
262 gretl_vector *b;
263 int i, j, nc;
264 int k = pmod->ncoeff;
265
266 if (select == NULL) {
267 nc = k;
268 } else {
269 nc = count_selection(select, k);
270 }
271
272 if (nc == 0) {
273 *err = E_DATA;
274 return NULL;
275 }
276
277 b = gretl_column_vector_alloc(nc);
278 if (b == NULL) {
279 *err = E_ALLOC;
280 return NULL;
281 }
282
283 j = 0;
284 for (i=0; i<k; i++) {
285 if (select != NULL && !select[i]) {
286 continue;
287 }
288 b->val[j++] = pmod->coeff[i];
289 }
290
291 return b;
292 }
293
294 /**
295 * gretl_covariance_matrix_from_varlist:
296 * @list: list of variables by ID number.
297 * @dset: dataset struct.
298 * @means: pointer to pick up vector of means, or NULL to discard.
299 * @errp: pointer to receive non-zero error code in case of
300 * failure, or NULL.
301 *
302 * Returns: the variance-covariance matrix of the listed variables
303 * (over the currently defined data sample), or NULL in case of
304 * failure.
305 */
306
307 gretl_matrix *
gretl_covariance_matrix_from_varlist(const int * list,const DATASET * dset,gretl_matrix ** means,int * errp)308 gretl_covariance_matrix_from_varlist (const int *list,
309 const DATASET *dset,
310 gretl_matrix **means,
311 int *errp)
312 {
313 gretl_matrix *V;
314 gretl_vector *xbar;
315 int k = list[0];
316 int t, i, j, vi, vj;
317 double vv, x, y;
318 int n, err = 0;
319
320 V = gretl_matrix_alloc(k, k);
321 xbar = gretl_vector_alloc(k);
322
323 if (V == NULL || xbar == NULL) {
324 err = E_ALLOC;
325 goto bailout;
326 }
327
328 for (i=0; i<k && !err; i++) {
329 xbar->val[i] = gretl_mean(dset->t1, dset->t2, dset->Z[list[i+1]]);
330 if (na(xbar->val[i])) {
331 err = E_DATA;
332 }
333 }
334
335 for (i=0; i<k && !err; i++) {
336 vi = list[i+1];
337 for (j=i; j<k; j++) {
338 vj = list[j+1];
339 vv = 0.0;
340 n = 0;
341 for (t=dset->t1; t<=dset->t2; t++) {
342 x = dset->Z[vi][t];
343 y = dset->Z[vj][t];
344 if (na(x) || na(y)) {
345 continue;
346 }
347 vv += (x - xbar->val[i]) * (y - xbar->val[j]);
348 n++;
349 }
350 if (n < 2) {
351 err = E_DATA;
352 vv = NADBL;
353 } else {
354 vv /= (n - 1); /* or plain n? */
355 }
356 gretl_matrix_set(V, i, j, vv);
357 gretl_matrix_set(V, j, i, vv);
358 }
359 }
360
361 bailout:
362
363 if (means != NULL && !err) {
364 *means = xbar;
365 } else {
366 gretl_vector_free(xbar);
367 }
368
369 if (err) {
370 gretl_matrix_free(V);
371 V = NULL;
372 if (errp != NULL) {
373 *errp = err;
374 }
375 }
376
377 return V;
378 }
379
380 /**
381 * gretl_matrix_row_to_array:
382 * @m: source matrix.
383 * @i: the row from which values should be copied.
384 * @x: array of doubles large enough to hold a row from @m.
385 *
386 * Copies the values from row @i of matrix @m into the array
387 * @x, which should already be allocated to the correct size.
388 *
389 * Returns: 0 on sucess, non-zero if the row is out of bounds.
390 */
391
gretl_matrix_row_to_array(const gretl_matrix * m,int i,double * x)392 int gretl_matrix_row_to_array (const gretl_matrix *m, int i, double *x)
393 {
394 int j, err = 0;
395
396 if (i < 0 || i >= gretl_matrix_rows(m)) {
397 err = E_DATA;
398 } else {
399 for (j=0; j<m->cols; j++) {
400 x[j] = gretl_matrix_get(m, i, j);
401 }
402 }
403
404 return err;
405 }
406
407 /**
408 * gretl_matrix_get_columns:
409 * @m: source matrix.
410 * @err: location to receive error code.
411 *
412 * Constructs a two-dimensional array in which each sub-array
413 * is a pointer to a column of @m. The content of these arrays
414 * belongs to @m and must not be freed; only the returned
415 * pointer itself should be freed.
416 *
417 * Returns: the constructed array on success or NULL on failure.
418 */
419
gretl_matrix_get_columns(const gretl_matrix * m,int * err)420 double **gretl_matrix_get_columns (const gretl_matrix *m, int *err)
421 {
422 double **X = NULL;
423
424 if (gretl_is_null_matrix(m)) {
425 *err = E_DATA;
426 } else {
427 double *val = m->val;
428 int j;
429
430 X = doubles_array_new(m->cols, 0);
431
432 if (X == NULL) {
433 *err = E_ALLOC;
434 } else {
435 for (j=0; j<m->cols; j++) {
436 X[j] = val;
437 val += m->rows;
438 }
439 }
440 }
441
442 return X;
443 }
444
get_mask_count(const char * mask,int n)445 static int get_mask_count (const char *mask, int n)
446 {
447 int i, k = 0;
448
449 for (i=0; i<n; i++) {
450 if (mask[i]) k++;
451 }
452
453 return k;
454 }
455
add_dataset_colnames(gretl_matrix * M,const int * list,const DATASET * dset)456 static void add_dataset_colnames (gretl_matrix *M,
457 const int *list,
458 const DATASET *dset)
459 {
460 int i, vi, nv;
461 char **S = NULL;
462 int err = 0;
463
464 if (list != NULL) {
465 nv = list[0];
466 } else {
467 /* all series except "const" */
468 nv = dset->v - 1;
469 }
470
471 /* we won't treat errors here as fatal */
472
473 if (nv != M->cols) {
474 /* "can't happen" */
475 return;
476 }
477
478 S = strings_array_new(nv);
479 if (S == NULL) {
480 return;
481 }
482
483 for (i=1; i<=nv; i++) {
484 vi = list == NULL ? i : list[i];
485 S[i-1] = gretl_strdup(dset->varname[vi]);
486 if (S[i-1] == NULL) {
487 strings_array_free(S, nv);
488 return;
489 }
490 }
491
492 err = gretl_matrix_set_colnames(M, S);
493 if (err) {
494 strings_array_free(S, nv);
495 }
496 }
497
498 static gretl_matrix *
real_gretl_matrix_data_subset(const int * list,const DATASET * dset,int t1,int t2,const char * mask,int op,int * err)499 real_gretl_matrix_data_subset (const int *list,
500 const DATASET *dset,
501 int t1, int t2,
502 const char *mask,
503 int op, int *err)
504 {
505 gretl_matrix *M = NULL;
506 int T, Tmax = t2 - t1 + 1;
507 int do_obs_info = 1;
508 int k, mt1 = 0, mt2 = 0;
509 int skip;
510 int j, vj, s, t;
511
512 if (list != NULL) {
513 k = list[0];
514 } else {
515 k = dset->v - 1;
516 }
517
518 if (t1 < 0 && t2 < 0) {
519 /* use full dataset */
520 t1 = 0;
521 t2 = dset->n - 1;
522 T = Tmax = dset->n;
523 do_obs_info = 0;
524 }
525
526 if (k <= 0 || Tmax <= 0 || t2 >= dset->n) {
527 *err = E_DATA;
528 return NULL;
529 }
530
531 *err = 0;
532 T = Tmax;
533
534 if (mask != NULL) {
535 T -= get_mask_count(mask, Tmax);
536 } else if (op == M_MISSING_TRIM) {
537 *err = list_adjust_sample(list, &t1, &t2, dset, NULL);
538 if (*err) {
539 return NULL;
540 } else {
541 Tmax = T = t2 - t1 + 1;
542 }
543 } else if (op == M_MISSING_SKIP || op == M_MISSING_ERROR) {
544 for (t=t1; t<=t2 && !*err; t++) {
545 for (j=1; j<=k; j++) {
546 vj = list == NULL ? j : list[j];
547 if (na(dset->Z[vj][t])) {
548 if (op == M_MISSING_SKIP) {
549 T--;
550 } else {
551 *err = E_MISSDATA;
552 return NULL;
553 }
554 break;
555 }
556 }
557 }
558 }
559
560 if (T < 0) {
561 /* "can't happen" */
562 *err = E_DATA;
563 } else {
564 M = gretl_matrix_alloc(T, k);
565 if (M == NULL) {
566 *err = E_ALLOC;
567 }
568 }
569
570 if (*err || T == 0) {
571 return M;
572 }
573
574 if (do_obs_info) {
575 mt1 = t1; mt2 = t2;
576 }
577
578 s = 0;
579 for (t=t1; t<=t2; t++) {
580 skip = 0;
581 if (mask != NULL) {
582 skip = mask[t - t1];
583 } else if (op == M_MISSING_SKIP) {
584 for (j=1; j<=k; j++) {
585 vj = list == NULL ? j : list[j];
586 if (na(dset->Z[vj][t])) {
587 skip = 1;
588 break;
589 }
590 }
591 }
592 if (!skip) {
593 for (j=1; j<=k; j++) {
594 vj = list == NULL ? j : list[j];
595 gretl_matrix_set(M, s, j-1, dset->Z[vj][t]);
596 }
597 if (s == 0) {
598 mt1 = t;
599 } else if (s == T - 1) {
600 mt2 = t;
601 break;
602 }
603 s++;
604 }
605 }
606
607 if (!*err && (mask != NULL || op == M_MISSING_OK)) {
608 int n = k * T;
609
610 for (j=0; j<n && !*err; j++) {
611 if (na(M->val[j])) {
612 if (op != M_MISSING_OK) {
613 *err = E_MISSDATA;
614 }
615 }
616 }
617 }
618
619 if (*err) {
620 gretl_matrix_free(M);
621 M = NULL;
622 } else {
623 if (do_obs_info && T == mt2 - mt1 + 1) {
624 gretl_matrix_set_t1(M, mt1);
625 gretl_matrix_set_t2(M, mt2);
626 }
627 if (dset->varname != NULL) {
628 add_dataset_colnames(M, list, dset);
629 }
630 }
631
632 return M;
633 }
634
635 /**
636 * gretl_matrix_data_subset:
637 * @list: list of series to process, or %NULL to include
638 * all series except "const".
639 * @dset: pointer to dataset struct.
640 * @t1: starting observation.
641 * @t2: ending observation.
642 * @missop: how to handle missing observations.
643 * @err: location to receive error code.
644 *
645 * Creates a matrix holding the subset of variables from
646 * @dset specified by @list, over the sample range @t1 to @t2,
647 * inclusive. Series are in columns. The @missop flag can
648 * be %M_MISSING_OK to indicate that it's OK to include
649 * missing values in the matrix, %M_MISSING_ERROR (it's an
650 * error if any missing values are found), or %M_MISSING_SKIP
651 * (observations with any missing values are omitted from the
652 * matrix).
653 *
654 * Returns: allocated matrix or NULL on failure.
655 */
656
gretl_matrix_data_subset(const int * list,const DATASET * dset,int t1,int t2,int missop,int * err)657 gretl_matrix *gretl_matrix_data_subset (const int *list,
658 const DATASET *dset,
659 int t1, int t2,
660 int missop,
661 int *err)
662 {
663 return real_gretl_matrix_data_subset(list, dset, t1, t2, NULL,
664 missop, err);
665 }
666
667 /**
668 * gretl_matrix_data_subset_masked:
669 * @list: list of variable to process.
670 * @dset: dataset struct.
671 * @t1: starting observation.
672 * @t2: ending observation.
673 * @mask: missing observations mask, or NULL.
674 * @err: location to receive error code.
675 *
676 * Creates a gretl matrix holding the subset of variables from
677 * @dset specified by @list, over the sample range @t1 to @t2,
678 * inclusive. Variables are in columns. @mask should be an
679 * array of char of length (@t2 - @t1 + 1) with 1s in the positions
680 * of observations to exclude from the subset and zeros elsewhere.
681 * This apparatus can be used to exclude missing observations.
682 *
683 * Returns: allocated matrix or NULL on failure.
684 */
685
686 gretl_matrix *
gretl_matrix_data_subset_masked(const int * list,const DATASET * dset,int t1,int t2,const char * mask,int * err)687 gretl_matrix_data_subset_masked (const int *list,
688 const DATASET *dset,
689 int t1, int t2,
690 const char *mask,
691 int *err)
692 {
693 return real_gretl_matrix_data_subset(list, dset, t1, t2, mask,
694 M_MISSING_ERROR, err);
695 }
696
697 /* count the number of selected rows in the current
698 sample range */
699
mmask_row_count(const gretl_matrix * mask,const DATASET * dset)700 static int mmask_row_count (const gretl_matrix *mask,
701 const DATASET *dset)
702 {
703 int t, m = 0;
704
705 for (t=dset->t1; t<=dset->t2; t++) {
706 if (mask->val[t] != 0) {
707 m++;
708 }
709 }
710
711 return m;
712 }
713
714 /**
715 * gretl_matrix_data_subset_special:
716 * @list: list of variables to process.
717 * @dset: dataset struct.
718 * @mmask: matrix holding desired observations mask.
719 * @err: location to receive error code.
720 *
721 * Creates a gretl matrix holding the subset of variables from
722 * @dset specified by @list, using the observations (data rows)
723 * selected by non-zero elements in @mmask. This is designed
724 * to support the gretl "libset" variable %matrix_mask.
725 *
726 * The length of the vector @mmask must equal the number of
727 * observations in the dataset, the member %n of @dset.
728 *
729 * Returns: allocated matrix or NULL on failure.
730 */
731
732 gretl_matrix *
gretl_matrix_data_subset_special(const int * list,const DATASET * dset,const gretl_matrix * mmask,int * err)733 gretl_matrix_data_subset_special (const int *list,
734 const DATASET *dset,
735 const gretl_matrix *mmask,
736 int *err)
737 {
738 int n = gretl_vector_get_length(mmask);
739 gretl_matrix *X = NULL;
740
741 if (list == NULL || n != dset->n) {
742 *err = E_DATA;
743 } else if (list[0] == 0) {
744 X = gretl_null_matrix_new();
745 } else {
746 int T = mmask_row_count(mmask, dset);
747 int k = list[0];
748
749 if (T == 0) {
750 X = gretl_null_matrix_new();
751 } else {
752 X = gretl_matrix_alloc(T, k);
753 }
754
755 if (X != NULL && T > 0) {
756 const double *xi;
757 int i, s, t;
758
759 for (i=0; i<k; i++) {
760 xi = dset->Z[list[i+1]];
761 s = 0;
762 for (t=dset->t1; t<=dset->t2; t++) {
763 if (mmask->val[t] != 0) {
764 if (s == 0) {
765 gretl_matrix_set_t1(X, t);
766 } else if (s == T - 1) {
767 gretl_matrix_set_t2(X, t);
768 }
769 gretl_matrix_set(X, s++, i, xi[t]);
770 }
771 }
772 }
773 }
774 }
775
776 if (X == NULL && *err == 0) {
777 *err = E_ALLOC;
778 }
779
780 return X;
781 }
782
check_matrix_varnames(DATASET * dset)783 static void check_matrix_varnames (DATASET *dset)
784 {
785 int i, j, err = 0;
786
787 for (i=1; i<dset->v; i++) {
788 for (j=1; j<i; j++) {
789 if (!strcmp(dset->varname[i], dset->varname[j])) {
790 err = 1;
791 break;
792 }
793 }
794 }
795
796 if (err) {
797 for (i=1; i<dset->v; i++) {
798 sprintf(dset->varname[i], "v%d", i);
799 }
800 }
801 }
802
803 /**
804 * gretl_dataset_from_matrix:
805 * @m: source matrix.
806 * @list: list of columns (1-based) to include, or NULL.
807 * @opt: may include OPT_B to attempt "borrowing" of data;
808 * OPT_N to use plain numbers as variable names; OPT_R to
809 * use row-names as observation markers, if present; OPT_S
810 * when called from the "store" command.
811 * @err: location to receive error code.
812 *
813 * Creates a gretl dataset from matrix @m, either using the
814 * columns specified in @list or using all columns if @list
815 * is NULL.
816 *
817 * Returns: pointer to new dataset information struct on success,
818 * or NULL on failure.
819 */
820
gretl_dataset_from_matrix(const gretl_matrix * m,const int * list,gretlopt opt,int * err)821 DATASET *gretl_dataset_from_matrix (const gretl_matrix *m,
822 const int *list,
823 gretlopt opt,
824 int *err)
825 {
826 DATASET *dset = NULL;
827 const char **cnames = NULL;
828 const char **rnames = NULL;
829 int use_rownames;
830 int i, t, col, nv, T;
831
832 if (gretl_is_null_matrix(m)) {
833 *err = E_DATA;
834 return NULL;
835 }
836
837 use_rownames = opt & (OPT_R | OPT_S);
838
839 T = m->rows;
840 nv = m->cols;
841
842 if (list != NULL) {
843 for (i=1; i<=list[0]; i++) {
844 col = list[i];
845 if (col < 1 || col > nv) {
846 gretl_errmsg_sprintf("Variable number %d is out of bounds", col);
847 *err = E_DATA;
848 break;
849 }
850 }
851 nv = list[0];
852 }
853
854 if (!*err) {
855 gretlopt dsopt = OPT_NONE;
856
857 if (opt & OPT_B) {
858 dsopt |= OPT_B;
859 }
860 if (use_rownames) {
861 rnames = gretl_matrix_get_rownames(m);
862 if (rnames != NULL) {
863 dsopt |= OPT_M;
864 }
865 }
866 dset = create_auxiliary_dataset(nv + 1, T, dsopt);
867 if (dset == NULL) {
868 *err = E_ALLOC;
869 }
870 }
871
872 if (*err) {
873 return NULL;
874 }
875
876 cnames = gretl_matrix_get_colnames(m);
877
878 for (i=1; i<=nv; i++) {
879 double *src;
880
881 col = (list != NULL)? list[i] - 1 : i - 1;
882 src = m->val + T * col;
883 if (opt & OPT_B) {
884 /* "borrowing" */
885 dset->Z[i] = src;
886 } else {
887 /* copying */
888 memcpy(dset->Z[i], src, T * sizeof *src);
889 }
890 if (cnames != NULL) {
891 if (gretl_normalize_varname(dset->varname[i], cnames[col], 0, i)) {
892 series_record_display_name(dset, i, cnames[col]);
893 }
894 } else if (opt & OPT_N) {
895 sprintf(dset->varname[i], "%d", col + 1);
896 } else if (opt & OPT_R) {
897 sprintf(dset->varname[i], "v%d", i);
898 } else {
899 sprintf(dset->varname[i], "col%d", col + 1);
900 }
901 }
902
903 if (rnames != NULL) {
904 for (t=0; t<T; t++) {
905 dset->S[t][0] = '\0';
906 strncat(dset->S[t], rnames[t], OBSLEN-1);
907 }
908 }
909
910 if (cnames != NULL && (opt & OPT_S)) {
911 /* we must have valid varnames for "store" */
912 check_matrix_varnames(dset);
913 }
914
915 return dset;
916 }
917
write_matrix_as_dataset(const char * fname,gretlopt opt,PRN * prn)918 int write_matrix_as_dataset (const char *fname,
919 gretlopt opt,
920 PRN *prn)
921 {
922 const char *mname;
923 gretl_matrix *m;
924 DATASET *mdset;
925 int err = 0;
926
927 mname = get_optval_string(STORE, OPT_A);
928 m = get_matrix_by_name(mname);
929 if (m == NULL) {
930 return E_DATA;
931 }
932
933 mdset = gretl_dataset_from_matrix(m, NULL, OPT_S, &err);
934 if (!err) {
935 opt &= ~OPT_A;
936 err = write_data(fname, NULL, mdset, opt, prn);
937 }
938
939 destroy_dataset(mdset);
940
941 return err;
942 }
943
944 /**
945 * gretl_plotfit_matrices:
946 * @yvar: the y variable.
947 * @xvar: the x variable.
948 * @fit: type of fit sought.
949 * @t1: starting observation.
950 * @t2: ending observation.
951 * @py: location to receive y vector.
952 * @pX: location to receive X matrix.
953 *
954 * Creates a vector y and matrix X based on the input @yvar,
955 * @xvar and @fit, using the given sample range. An observation
956 * is skipped if either @yvar or @xvar is missing at that
957 * observation.
958 *
959 * Returns: 0 on success, non-zero code on error.
960 */
961
gretl_plotfit_matrices(const double * yvar,const double * xvar,FitType fit,int t1,int t2,gretl_matrix ** py,gretl_matrix ** pX)962 int gretl_plotfit_matrices (const double *yvar, const double *xvar,
963 FitType fit, int t1, int t2,
964 gretl_matrix **py, gretl_matrix **pX)
965 {
966 gretl_matrix *y = NULL;
967 gretl_matrix *X = NULL;
968 char *mask = NULL;
969 double xt;
970 int T = t2 - t1 + 1;
971 int n = 0;
972 int i, j, k, s, t;
973 int err = 0;
974
975 if (T <= 0) {
976 return E_DATA;
977 }
978
979 if (fit == PLOT_FIT_LOGLIN && !gretl_ispositive(t1, t2, yvar, 1)) {
980 gretl_errmsg_set(_("Non-positive values encountered"));
981 return E_DATA;
982 } else if (fit == PLOT_FIT_LINLOG && !gretl_ispositive(t1, t2, xvar, 1)) {
983 gretl_errmsg_set(_("Non-positive values encountered"));
984 return E_DATA;
985 }
986
987 mask = calloc(T, 1);
988 if (mask == NULL) {
989 return E_ALLOC;
990 }
991
992 for (s=0; s<T; s++) {
993 t = s + t1;
994 if (na(yvar[t]) || (xvar != NULL && na(xvar[t]))) {
995 mask[s] = 1;
996 } else {
997 n++;
998 }
999 }
1000
1001 if (n == 0) {
1002 free(mask);
1003 return E_MISSDATA;
1004 }
1005
1006 if (fit == PLOT_FIT_CUBIC) {
1007 k = 4;
1008 } else if (fit == PLOT_FIT_QUADRATIC) {
1009 k = 3;
1010 } else if (fit == PLOT_FIT_LOESS) {
1011 k = 1;
1012 } else {
1013 k = 2;
1014 }
1015
1016 y = gretl_column_vector_alloc(n);
1017 X = gretl_matrix_alloc(n, k);
1018 if (y == NULL || X == NULL) {
1019 err = E_ALLOC;
1020 goto bailout;
1021 }
1022
1023 i = 0;
1024 for (s=0; s<T; s++) {
1025 t = s + t1;
1026 if (!mask[s]) {
1027 j = 0;
1028 if (fit == PLOT_FIT_LOGLIN) {
1029 y->val[i] = log(yvar[t]);
1030 } else {
1031 y->val[i] = yvar[t];
1032 }
1033 if (fit != PLOT_FIT_LOESS) {
1034 gretl_matrix_set(X, i, j++, 1.0);
1035 }
1036 xt = (xvar != NULL)? xvar[t] : s;
1037 if (fit == PLOT_FIT_INVERSE) {
1038 gretl_matrix_set(X, i, j++, 1.0 / xt);
1039 } else if (fit == PLOT_FIT_LINLOG) {
1040 gretl_matrix_set(X, i, j++, log(xt));
1041 } else {
1042 gretl_matrix_set(X, i, j++, xt);
1043 }
1044 if (fit == PLOT_FIT_QUADRATIC || fit == PLOT_FIT_CUBIC) {
1045 gretl_matrix_set(X, i, j++, xt * xt);
1046 }
1047 if (fit == PLOT_FIT_CUBIC) {
1048 gretl_matrix_set(X, i, j, xt * xt * xt);
1049 }
1050 i++;
1051 }
1052 }
1053
1054 bailout:
1055
1056 free(mask);
1057
1058 if (err) {
1059 gretl_matrix_free(y);
1060 gretl_matrix_free(X);
1061 } else {
1062 *py = y;
1063 *pX = X;
1064 }
1065
1066 return err;
1067 }
1068
skip_matrix_comment(FILE * fp,int * rows,int * cols,int * is_complex,int * err)1069 static int skip_matrix_comment (FILE *fp, int *rows, int *cols,
1070 int *is_complex, int *err)
1071 {
1072 int c, ret = 0;
1073
1074 c = fgetc(fp);
1075
1076 if (c == '#') {
1077 char *p, test[16] = {0};
1078 int n, i = 0;
1079
1080 ret = 1;
1081 while (c != '\n' && c != EOF && c != -1) {
1082 c = fgetc(fp);
1083 if (i < 15) {
1084 test[i] = c;
1085 }
1086 i++;
1087 }
1088 if ((p = strstr(test, "rows:")) != NULL) {
1089 if (sscanf(p + 5, "%d", &n) == 1) {
1090 *rows = n;
1091 }
1092 } else if ((p = strstr(test, "columns:")) != NULL) {
1093 if (sscanf(p + 8, "%d", &n) == 1) {
1094 *cols = n;
1095 }
1096 } else if ((p = strstr(test, "complex:")) != NULL) {
1097 if (sscanf(p + 8, "%d", &n) == 1) {
1098 *is_complex = (n == 1);
1099 }
1100 }
1101 } else {
1102 ungetc(c, fp);
1103 }
1104
1105 if (c == EOF || c == -1) {
1106 fprintf(stderr, "reached premature end of file\n");
1107 *err = E_DATA;
1108 }
1109
1110 return ret;
1111 }
1112
1113 #if G_BYTE_ORDER == G_BIG_ENDIAN
1114
matrix_swap_endianness(gretl_matrix * A)1115 static void matrix_swap_endianness (gretl_matrix *A)
1116 {
1117 int i, n = A->rows * A->cols;
1118
1119 for (i=0; i<n; i++) {
1120 reverse_double(A->val[i]);
1121 }
1122 }
1123
1124 #endif
1125
read_binary_matrix_file(FILE * fp,int * err)1126 static gretl_matrix *read_binary_matrix_file (FILE *fp, int *err)
1127 {
1128 gretl_matrix *A = NULL;
1129 char header[20] = {0};
1130 int is_complex = 0;
1131 gint32 dim[2];
1132
1133 if (fread(header, 1, 19, fp) < 19) {
1134 *err = E_DATA;
1135 } else if (fread(dim, sizeof *dim, 2, fp) < 2) {
1136 *err = E_DATA;
1137 } else if (dim[0] <= 0 || dim[1] <= 0) {
1138 *err = E_DATA;
1139 }
1140
1141 if (!*err) {
1142 if (!strcmp(header, "gretl_binary_matrix")) {
1143 ; /* OK */
1144 } else if (!strcmp(header, "gretl_binar_cmatrix")) {
1145 is_complex = 1;
1146 } else {
1147 *err = E_DATA;
1148 }
1149 }
1150
1151 if (!*err) {
1152 #if G_BYTE_ORDER == G_BIG_ENDIAN
1153 reverse_int(dim[0]);
1154 reverse_int(dim[1]);
1155 #endif
1156 A = gretl_matrix_alloc(dim[0], dim[1]);
1157 if (A == NULL) {
1158 *err = E_ALLOC;
1159 }
1160 }
1161
1162 if (!*err) {
1163 size_t n = dim[0] * dim[1];
1164
1165 if (fread(A->val, sizeof *A->val, n, fp) < n) {
1166 *err = E_DATA;
1167 gretl_matrix_free(A);
1168 A = NULL;
1169 } else if (is_complex) {
1170 gretl_matrix_set_complex_full(A, 1);
1171 }
1172 }
1173
1174 #if G_BYTE_ORDER == G_BIG_ENDIAN
1175 if (A != NULL) {
1176 matrix_swap_endianness(A);
1177 }
1178 #endif
1179
1180 fclose(fp);
1181
1182 return A;
1183 }
1184
1185 #ifndef WIN32
1186
1187 /* In reading matrices, accept "NA" and "." (Stata) for NaN? */
1188
unix_scan_NA(FILE * fp,int * err)1189 static double unix_scan_NA (FILE *fp, int *err)
1190 {
1191 char test[3] = {0};
1192
1193 if (fscanf(fp, "%2s", test) < 1) {
1194 *err = E_DATA;
1195 return 0;
1196 } else if (!strcmp(test, "NA") || !strcmp(test, ".")) {
1197 return NADBL;
1198 } else {
1199 gretl_errmsg_sprintf(_("got invalid field '%s'"), test);
1200 *err = E_DATA;
1201 return 0;
1202 }
1203 }
1204
1205 #endif
1206
1207 /**
1208 * gretl_matrix_read_from_file:
1209 * @fname: name of file.
1210 * @import: non-zero means we're importing via dotdir.
1211 * @err: location to receive error code.
1212 *
1213 * Reads a matrix from a file by the name @fname; if it's
1214 * a text file the column separator must be space or tab. It is
1215 * assumed that the dimensions of the matrix (number of rows and
1216 * columns) are found on the first line of the file, so no
1217 * heuristics are necessary. In case of error @err is filled
1218 * appropriately.
1219 *
1220 * Returns: The matrix as read from file, or NULL.
1221 */
1222
gretl_matrix_read_from_file(const char * fname,int import,int * err)1223 gretl_matrix *gretl_matrix_read_from_file (const char *fname,
1224 int import, int *err)
1225 {
1226 char fullname[FILENAME_MAX];
1227 gchar *unzname = NULL;
1228 int r = -1, c = -1, n = 0;
1229 int gz, bin = 0;
1230 gretl_matrix *A = NULL;
1231 int is_complex = 0;
1232 FILE *fp = NULL;
1233
1234 gz = has_suffix(fname, ".gz");
1235
1236 if (!gz) {
1237 bin = has_suffix(fname, ".bin");
1238 }
1239
1240 if (import) {
1241 gretl_build_path(fullname, gretl_dotdir(), fname, NULL);
1242 } else {
1243 strcpy(fullname, fname);
1244 gretl_maybe_prepend_dir(fullname);
1245 }
1246
1247 if (gz) {
1248 char tmp[FILENAME_MAX];
1249
1250 gretl_build_path(tmp, gretl_dotdir(), "_unztmp_.mat", NULL);
1251 *err = gretl_gunzip(fullname, tmp);
1252 if (!*err) {
1253 fp = gretl_fopen(tmp, "rb");
1254 unzname = g_strdup(tmp);
1255 }
1256 } else {
1257 fp = gretl_fopen(fullname, "rb");
1258 }
1259
1260 if (fp == NULL) {
1261 *err = E_FOPEN;
1262 g_free(unzname);
1263 return NULL;
1264 }
1265
1266 if (bin) {
1267 return read_binary_matrix_file(fp, err);
1268 }
1269
1270 /* 'Eat' any leading comment lines starting with '#',
1271 but while we go, scan for Matlab-style rows/columns
1272 specification.
1273 */
1274 while (!*err && skip_matrix_comment(fp, &r, &c, &is_complex, err)) {
1275 ;
1276 }
1277
1278 if (!*err) {
1279 if (r >= 0 && c >= 0) {
1280 /* we got dimensions from the preamble */
1281 n = 2;
1282 } else {
1283 n = fscanf(fp, "%d %ds\n", &r, &c);
1284 }
1285 if (n < 2 || r < 0 || c < 0) {
1286 fprintf(stderr, "error reading rows, cols (r=%d, c=%d)\n",
1287 r, c);
1288 *err = E_DATA;
1289 } else {
1290 A = gretl_matrix_alloc(r, c);
1291 if (A == NULL) {
1292 *err = E_ALLOC;
1293 }
1294 }
1295 }
1296
1297 if (!*err) {
1298 long pos = 0;
1299 double x;
1300 int i, j;
1301
1302 gretl_push_c_numeric_locale();
1303
1304 for (i=0; i<r && !*err; i++) {
1305 for (j=0; j<c && !*err; j++) {
1306 pos = ftell(fp);
1307 n = fscanf(fp, "%lf", &x);
1308 if (n != 1) {
1309 fseek(fp, pos, SEEK_SET);
1310 #ifdef WIN32
1311 x = win32_fscan_nonfinite(fp, err);
1312 #else
1313 x = unix_scan_NA(fp, err);
1314 #endif
1315 }
1316 if (*err) {
1317 fprintf(stderr, "error reading row %d, column %d\n", i+1, j+1);
1318 } else {
1319 gretl_matrix_set(A, i, j, x);
1320 }
1321 }
1322 }
1323
1324 gretl_pop_c_numeric_locale();
1325 }
1326
1327 fclose(fp);
1328
1329 if (unzname != NULL) {
1330 gretl_remove(unzname);
1331 g_free(unzname);
1332 }
1333
1334 if (!*err && is_complex) {
1335 gretl_matrix_set_complex_full(A, 1);
1336 }
1337
1338 if (*err && A != NULL) {
1339 gretl_matrix_free(A);
1340 A = NULL;
1341 }
1342
1343 return A;
1344 }
1345
1346 #ifdef WIN32
1347
win32_xna_out(double x,char pad,PRN * prn)1348 static void win32_xna_out (double x, char pad, PRN *prn)
1349 {
1350 if (isnan(x) || na(x)) {
1351 pputs(prn, "nan");
1352 } else if (x < 0) {
1353 pputs(prn, "-inf");
1354 } else {
1355 pputs(prn, "inf");
1356 }
1357 pputc(prn, pad);
1358 }
1359
1360 #endif
1361
matrix_to_csv(const gretl_matrix * A,const char * fname)1362 static int matrix_to_csv (const gretl_matrix *A, const char *fname)
1363 {
1364 const char *na_str;
1365 char delim;
1366 FILE *fp;
1367 double aij;
1368 int i, j;
1369
1370 fp = gretl_fopen(fname, "wb");
1371 if (fp == NULL) {
1372 return E_FOPEN;
1373 }
1374
1375 na_str = get_csv_na_write_string();
1376 delim = get_data_export_delimiter();
1377 if (delim == '\t') {
1378 /* we can't allow this! */
1379 delim = ',';
1380 }
1381
1382 gretl_push_c_numeric_locale();
1383
1384 for (i=0; i<A->rows; i++) {
1385 for (j=0; j<A->cols; j++) {
1386 aij = gretl_matrix_get(A, i, j);
1387 if (na(aij)) {
1388 fputs(na_str, fp);
1389 } else {
1390 fprintf(fp, "%.17g", aij);
1391 }
1392 if (j < A->cols - 1) {
1393 fputc(delim, fp);
1394 } else {
1395 fputc('\n', fp);
1396 }
1397 }
1398 }
1399
1400 gretl_pop_c_numeric_locale();
1401 fclose(fp);
1402
1403 return 0;
1404 }
1405
1406 /**
1407 * gretl_matrix_write_to_file:
1408 * @A: matrix to write.
1409 * @fname: name of file.
1410 * @export: non-zero means we're exporting via dotdir.
1411 *
1412 * Writes the matrix @A to a file by the name @fname; if it's a text
1413 * file, the column separator is tab. The number of rows and columns
1414 * are written on the first line of the file.
1415 *
1416 * Returns: 0 on successful completion, non-zero code on error.
1417 */
1418
gretl_matrix_write_to_file(gretl_matrix * A,const char * fname,int export)1419 int gretl_matrix_write_to_file (gretl_matrix *A, const char *fname,
1420 int export)
1421 {
1422 char targ[FILENAME_MAX];
1423 int r, c, i, j;
1424 int is_complex = 0;
1425 PRN *prn = NULL;
1426 FILE *fp = NULL;
1427 int csv = 0;
1428 int gz, bin = 0;
1429 int err = 0;
1430
1431 gz = has_suffix(fname, ".gz");
1432
1433 if (!gz) {
1434 bin = has_suffix(fname, ".bin");
1435 }
1436
1437 if (!gz && !bin) {
1438 csv = has_suffix(fname, ".csv");
1439 if (csv && A->is_complex) {
1440 return E_CMPLX;
1441 }
1442 }
1443
1444 if (export) {
1445 gretl_build_path(targ, gretl_dotdir(), fname, NULL);
1446 } else {
1447 fname = gretl_maybe_switch_dir(fname);
1448 strcpy(targ, fname);
1449 }
1450
1451 if (csv) {
1452 return matrix_to_csv(A, targ);
1453 }
1454
1455 if (bin) {
1456 fp = gretl_fopen(targ, "wb");
1457 } else if (gz) {
1458 prn = gretl_gzip_print_new(targ, -1, &err);
1459 } else {
1460 prn = gretl_print_new_with_filename(targ, &err);
1461 }
1462
1463 if (fp == NULL && prn == NULL) {
1464 return E_FOPEN;
1465 }
1466
1467 if (A->is_complex) {
1468 /* reset to "pretend real" for writing */
1469 is_complex = 1;
1470 gretl_matrix_set_complex_full(A, 0);
1471 }
1472
1473 r = A->rows;
1474 c = A->cols;
1475
1476 if (bin) {
1477 const char *header = is_complex ? "gretl_binar_cmatrix" :
1478 "gretl_binary_matrix";
1479 gint32 dim[2] = {r, c};
1480 size_t n = r * c;
1481
1482 #if G_BYTE_ORDER == G_BIG_ENDIAN
1483 double x;
1484 int k;
1485
1486 fwrite(header, 1, strlen(header), fp);
1487 for (i=0; i<2; i++) {
1488 k = dim[i];
1489 reverse_int(k);
1490 fwrite(&k, sizeof k, 1, fp);
1491 }
1492 for (i=0; i<n; i++) {
1493 x = A->val[i];
1494 reverse_double(x);
1495 fwrite(&x, sizeof x, 1, fp);
1496 }
1497 #else
1498 fwrite(header, 1, strlen(header), fp);
1499 fwrite(dim, sizeof *dim, 2, fp);
1500 fwrite(A->val, sizeof *A->val, n, fp);
1501 #endif
1502 fclose(fp);
1503 } else {
1504 /* !bin: textual representation */
1505 int format_g = libset_get_bool(MWRITE_G);
1506 char pad, d = '\t';
1507 double x;
1508
1509 if (is_complex) {
1510 pprintf(prn, "# rows: %d\n", r);
1511 pprintf(prn, "# columns: %d\n", c);
1512 pprintf(prn, "# complex: %d\n", 1);
1513 } else {
1514 pprintf(prn, "%d%c%d\n", r, d, c);
1515 }
1516
1517 gretl_push_c_numeric_locale();
1518
1519 for (i=0; i<r; i++) {
1520 for (j=0; j<c; j++) {
1521 pad = (j == c-1)? '\n' : d;
1522 x = gretl_matrix_get(A, i, j);
1523 #ifdef WIN32
1524 if (na(x)) {
1525 win32_xna_out(x, pad, prn);
1526 continue;
1527 }
1528 #endif
1529 if (format_g) {
1530 pprintf(prn, "%g", x);
1531 } else {
1532 pprintf(prn, "%26.18E", x);
1533 }
1534 pputc(prn, pad);
1535 }
1536 }
1537
1538 gretl_pop_c_numeric_locale();
1539 gretl_print_destroy(prn);
1540 }
1541
1542 if (is_complex) {
1543 /* reset A's original status */
1544 gretl_matrix_set_complex_full(A, 1);
1545 }
1546
1547 return err;
1548 }
1549
make_numstr(char * s,double x)1550 static void make_numstr (char *s, double x)
1551 {
1552 if (x == -0.0) {
1553 x = 0.0;
1554 }
1555
1556 if (isnan(x)) {
1557 strcpy(s, "nan");
1558 return;
1559 }
1560
1561 sprintf(s, "%#.5g", x);
1562
1563 /* remove surplus, or add deficient, zeros */
1564
1565 if (strstr(s, ".00000")) {
1566 s[strlen(s) - 1] = 0;
1567 } else {
1568 char *p = s;
1569 int n = 0;
1570
1571 while (*p) {
1572 if (isdigit(*p)) {
1573 n++;
1574 } else if (isalpha(*p)) {
1575 n = 0;
1576 break;
1577 }
1578 p++;
1579 }
1580 if (n > 0 && n < 5) {
1581 strcat(s, "0");
1582 }
1583 }
1584 }
1585
max_numchars(const gretl_matrix * m)1586 static int max_numchars (const gretl_matrix *m)
1587 {
1588 char s[24];
1589 int i, n = m->rows * m->cols;
1590 int c, cmax = 0;
1591
1592 for (i=0; i<n && cmax<6; i++) {
1593 sprintf(s, "%g", m->val[i]);
1594 c = strlen(s);
1595 if (c > cmax) {
1596 cmax = c;
1597 }
1598 }
1599
1600 return cmax;
1601 }
1602
max_label_length(const char ** names,int n)1603 static int max_label_length (const char **names, int n)
1604 {
1605 int i, len, maxlen = 0;
1606
1607 for (i=0; i<n; i++) {
1608 if (names[i] != NULL) {
1609 len = g_utf8_strlen(names[i], -1);
1610 if (len > maxlen) {
1611 maxlen = len;
1612 }
1613 }
1614 }
1615
1616 return maxlen;
1617 }
1618
1619 static void
real_matrix_print_to_prn(const gretl_matrix * m,const char * msg,int plain,const char ** colnames,const DATASET * dset,int imin,int imax,PRN * prn)1620 real_matrix_print_to_prn (const gretl_matrix *m,
1621 const char *msg,
1622 int plain,
1623 const char **colnames,
1624 const DATASET *dset,
1625 int imin, int imax,
1626 PRN *prn)
1627 {
1628 const char **rownames = NULL;
1629 char numstr[32];
1630 char obs[OBSLEN];
1631 double x;
1632 int strwidth = 12;
1633 int rnamelen = 0;
1634 int dated = 0;
1635 int cmax = 0, cpad = 2;
1636 int mt1 = 0, mt2 = 0;
1637 int i, j;
1638
1639 if (prn == NULL) {
1640 return;
1641 }
1642
1643 if (m == NULL) {
1644 if (msg != NULL && *msg != '\0') {
1645 pprintf(prn, _("%s: matrix is NULL"), msg);
1646 } else {
1647 pputs(prn, _("matrix is NULL"));
1648 }
1649 pputc(prn, '\n');
1650 return;
1651 } else if (m->rows == 0 || m->cols == 0) {
1652 if (msg != NULL && *msg != '\0') {
1653 pprintf(prn, _("%s: matrix is empty (%d x %d)"),
1654 msg, m->rows, m->cols);
1655 } else {
1656 pprintf(prn, _("matrix is empty (%d x %d)"),
1657 m->rows, m->cols);
1658 }
1659 pputc(prn, '\n');
1660 return;
1661 }
1662
1663 dated = gretl_matrix_is_dated(m);
1664 if (dated) {
1665 mt1 = gretl_matrix_get_t1(m);
1666 mt2 = gretl_matrix_get_t2(m);
1667 }
1668
1669 /* @plain means skip the header stuff */
1670
1671 if (msg != NULL && *msg != '\0' && !plain) {
1672 pprintf(prn, "%s (%d x %d)", msg, m->rows, m->cols);
1673 if (imin > 0 && imax > 0) {
1674 pprintf(prn, " rows %d to %d\n\n", imin+1, imax);
1675 } else if (dated && dset == NULL) {
1676 pprintf(prn, " [t1 = %d, t2 = %d]\n\n", mt1 + 1, mt2 + 1);
1677 } else {
1678 pputs(prn, "\n\n");
1679 }
1680 }
1681
1682 cmax = max_numchars(m);
1683 if (cmax > 5) {
1684 cmax = 0;
1685 }
1686
1687 if (colnames == NULL) {
1688 colnames = gretl_matrix_get_colnames(m);
1689 }
1690
1691 rownames = gretl_matrix_get_rownames(m);
1692
1693 if (rownames != NULL) {
1694 rnamelen = max_label_length(rownames, m->rows);
1695 } else if (dated && dset != NULL) {
1696 rnamelen = max_obs_marker_length(dset);
1697 }
1698
1699 if (colnames != NULL) {
1700 int numeric_names = 1;
1701
1702 if (rnamelen > 0) {
1703 bufspace(rnamelen + 1, prn);
1704 }
1705 for (j=0; j<m->cols; j++) {
1706 pprintf(prn, "%*.*s ", strwidth, strwidth, colnames[j]);
1707 if (!numeric_string(colnames[j])) {
1708 numeric_names = 0;
1709 }
1710 }
1711 pputc(prn, '\n');
1712 if (numeric_names) {
1713 pputc(prn, '\n');
1714 }
1715 }
1716
1717 if (imin < 0) imin = 0;
1718 if (imax < 0) imax = m->rows;
1719
1720 for (i=imin; i<imax; i++) {
1721 if (rnamelen > 0) {
1722 if (rownames != NULL) {
1723 pprintf(prn, "%*s ", rnamelen, rownames[i]);
1724 } else {
1725 ntolabel(obs, i + mt1, dset);
1726 pprintf(prn, "%*s ", rnamelen, obs);
1727 }
1728 }
1729 for (j=0; j<m->cols; j++) {
1730 x = gretl_matrix_get(m, i, j);
1731 if (cmax > 0) {
1732 if (colnames != NULL) {
1733 bufspace(strwidth - cmax - cpad, prn);
1734 }
1735 if (isnan(x)) {
1736 strcpy(numstr, "nan");
1737 } else {
1738 sprintf(numstr, "%g", x);
1739 }
1740 pprintf(prn, "%*s ", cmax + cpad, numstr);
1741 } else {
1742 make_numstr(numstr, x);
1743 pprintf(prn, "%*s ", strwidth, numstr);
1744 }
1745 }
1746 pputc(prn, '\n');
1747 }
1748
1749 pputc(prn, '\n');
1750 }
1751
1752 /**
1753 * gretl_matrix_print_to_prn:
1754 * @m: matrix to print.
1755 * @msg: accompanying message text (or NULL if no message is wanted).
1756 * @prn: pointer to gretl printing struct.
1757 *
1758 * Prints the matrix @m to @prn.
1759 */
1760
1761 void
gretl_matrix_print_to_prn(const gretl_matrix * m,const char * msg,PRN * prn)1762 gretl_matrix_print_to_prn (const gretl_matrix *m,
1763 const char *msg, PRN *prn)
1764 {
1765 real_matrix_print_to_prn(m, msg, 0, NULL, NULL,
1766 -1, -1, prn);
1767 }
1768
1769 /**
1770 * gretl_matrix_print_range:
1771 * @m: matrix to print.
1772 * @rmin: first row to print.
1773 * @prn: last row to print.
1774 *
1775 * Prints the specified rows of matrix @m to @prn.
1776 */
1777
gretl_matrix_print_range(const gretl_matrix * m,const char * msg,int rmin,int rmax,PRN * prn)1778 void gretl_matrix_print_range (const gretl_matrix *m,
1779 const char *msg,
1780 int rmin, int rmax,
1781 PRN *prn)
1782 {
1783 if (m->is_complex) {
1784 gretl_cmatrix_print_range(m, msg, rmin, rmax, prn);
1785 } else {
1786 real_matrix_print_to_prn(m, msg, 0, NULL, NULL,
1787 rmin, rmax, prn);
1788 }
1789 }
1790
1791 /**
1792 * gretl_matrix_print_with_col_heads:
1793 * @m: T x k matrix to print.
1794 * @title: accompanying title (or NULL if no title is wanted).
1795 * @heads: array of k strings to identify the columns.
1796 * @prn: pointer to gretl printing struct.
1797 *
1798 * Prints the matrix @m to @prn, with column headings given
1799 * by @heads.
1800 */
1801
gretl_matrix_print_with_col_heads(const gretl_matrix * m,const char * title,const char ** heads,const DATASET * dset,PRN * prn)1802 void gretl_matrix_print_with_col_heads (const gretl_matrix *m,
1803 const char *title,
1804 const char **heads,
1805 const DATASET *dset,
1806 PRN *prn)
1807 {
1808 real_matrix_print_to_prn(m, title, 0, heads, dset,
1809 -1, -1, prn);
1810 }
1811
maybe_print_col_heads(const gretl_matrix * m,const char * fmt,int wid,int prec,int icast,int llen,PRN * prn)1812 static void maybe_print_col_heads (const gretl_matrix *m,
1813 const char *fmt,
1814 int wid, int prec,
1815 int icast, int llen,
1816 PRN *prn)
1817 {
1818 const char **heads;
1819
1820 heads = gretl_matrix_get_colnames(m);
1821
1822 if (heads != NULL) {
1823 int numeric = 1;
1824 char wtest[32];
1825 double x;
1826 int j, n;
1827
1828 x = gretl_matrix_get(m, 0, 0);
1829
1830 if (icast) {
1831 if (wid >= 0 && prec >= 0) {
1832 snprintf(wtest, 32, fmt, wid, prec, (int) x);
1833 } else if (wid >= 0 || prec >= 0) {
1834 n = (wid >= 0)? wid : prec;
1835 snprintf(wtest, 32, fmt, n, (int) x);
1836 } else {
1837 snprintf(wtest, 32, fmt, (int) x);
1838 }
1839 } else {
1840 if (wid >= 0 && prec >= 0) {
1841 snprintf(wtest, 32, fmt, wid, prec, x);
1842 } else if (wid >= 0 || prec >= 0) {
1843 n = (wid >= 0)? wid : prec;
1844 snprintf(wtest, 32, fmt, n, x);
1845 } else {
1846 snprintf(wtest, 32, fmt, x);
1847 }
1848 }
1849
1850 if (llen > 0) {
1851 bufspace(llen + 1, prn);
1852 }
1853
1854 n = strlen(wtest);
1855 for (j=0; j<m->cols; j++) {
1856 pprintf(prn, "%*s", n, heads[j]);
1857 if (!numeric_string(heads[j])) {
1858 numeric = 0;
1859 }
1860 }
1861 pputc(prn, '\n');
1862 if (numeric) {
1863 pputc(prn, '\n');
1864 }
1865 }
1866 }
1867
1868 /* Return an array of char holding 1s for columns in
1869 @m that are all integers, 0s otherwise. Except
1870 that if we find that either none or all the columns
1871 if @m are composed of integers we scratch the
1872 array and signal the result via the @intcast pointer.
1873 */
1874
get_intcols(const gretl_matrix * m,int * intcast)1875 static char *get_intcols (const gretl_matrix *m,
1876 int *intcast)
1877 {
1878 char *ret = malloc(m->cols);
1879
1880 if (ret != NULL) {
1881 double x;
1882 int i, j, n = 0;
1883
1884 for (j=0; j<m->cols; j++) {
1885 ret[j] = 1;
1886 for (i=0; i<m->rows; i++) {
1887 x = gretl_matrix_get(m, i, j);
1888 if (x != floor(x)) {
1889 ret[j] = 0;
1890 break;
1891 }
1892 }
1893 n += ret[j];
1894 }
1895 if (n == 0 || n == m->cols) {
1896 /* no, or all, integer columns */
1897 free(ret);
1898 ret = NULL;
1899 *intcast = n > 0;
1900 }
1901 }
1902
1903 return ret;
1904 }
1905
1906 /* Take a format which is designed for floating-point
1907 values and convert it for use with integers: preserve
1908 the width but discard the precision, and change the
1909 trailing 'f', 'g' or whatever to 'd'.
1910 */
1911
revise_integer_format(char * ifmt)1912 static void revise_integer_format (char *ifmt)
1913 {
1914 char *p = strchr(ifmt, '.');
1915
1916 if (p != NULL) {
1917 *p = '\0';
1918 strcat(p, "d");
1919 } else {
1920 ifmt[strlen(ifmt)-1] = 'd';
1921 }
1922 }
1923
1924 /**
1925 * gretl_matrix_print_with_format:
1926 * @m: matrix to print.
1927 * @fmt: a printf-type format string.
1928 * @wid: an integer width, or 0.
1929 * @prec: an integer precision, or 0.
1930 * @prn: pointer to gretl printing struct.
1931 *
1932 * Prints the matrix @m to @prn, with the elements of @m
1933 * formatted as specified by @fmt, @wid and @prec.
1934 * The arguments @wid and/or @prec are required only if
1935 * @fmt contains one or more placeholders ("*") to be filled out.
1936 * For example, if @fmt is "\%14.6f" then neither @wid nor
1937 * @prec is needed and both should be set to 0; if
1938 * @fmt is "\%*.6f" then a @wid value is needed but @prec
1939 * should be 0.
1940 */
1941
gretl_matrix_print_with_format(const gretl_matrix * m,const char * fmt,int wid,int prec,PRN * prn)1942 void gretl_matrix_print_with_format (const gretl_matrix *m,
1943 const char *fmt,
1944 int wid, int prec,
1945 PRN *prn)
1946 {
1947 if (prn == NULL) {
1948 return;
1949 }
1950
1951 if (gretl_is_null_matrix(m) || fmt == NULL || *fmt == '\0') {
1952 real_matrix_print_to_prn(m, NULL, 1, NULL, NULL, -1, -1, prn);
1953 } else if (m->is_complex) {
1954 gretl_cmatrix_printf(m, fmt, prn);
1955 } else {
1956 const char **rownames = NULL;
1957 char *ifmt = NULL, *xfmt = NULL;
1958 char *intcols = NULL;
1959 int llen = 0, intcast = 0;
1960 int cpos = strlen(fmt) - 1;
1961 int variant = 0;
1962 char c = fmt[cpos];
1963 int i, j, k;
1964 double x;
1965
1966 xfmt = gretl_strdup(fmt);
1967
1968 if (c == 'd' || c == 'u' || c == 'x' || c == 'l') {
1969 intcast = 1;
1970 } else if (fmt[cpos-1] == 'v') {
1971 /* new-style "variant" column format */
1972 intcols = get_intcols(m, &intcast);
1973 gretl_delchar('v', xfmt);
1974 variant = 1;
1975 } else if (c == 'v') {
1976 /* old-style "variant" column format */
1977 intcols = get_intcols(m, &intcast);
1978 xfmt[cpos] = 'g';
1979 variant = 1;
1980 }
1981
1982 if (intcast || intcols) {
1983 ifmt = gretl_strdup(xfmt);
1984 if (variant) {
1985 revise_integer_format(ifmt);
1986 }
1987 }
1988
1989 rownames = gretl_matrix_get_rownames(m);
1990 if (rownames != NULL) {
1991 llen = max_label_length(rownames, m->rows);
1992 }
1993
1994 maybe_print_col_heads(m, xfmt, wid, prec, intcast, llen, prn);
1995
1996 for (i=0; i<m->rows; i++) {
1997 if (rownames != NULL) {
1998 pprintf(prn, "%*s ", llen, rownames[i]);
1999 }
2000 for (j=0; j<m->cols; j++) {
2001 x = gretl_matrix_get(m, i, j);
2002 if (intcols != NULL) {
2003 intcast = intcols[j];
2004 }
2005 if (intcast) {
2006 if (wid > 0 && prec > 0) {
2007 pprintf(prn, ifmt, wid, prec, (int) x);
2008 } else if (wid > 0 || prec > 0) {
2009 k = (wid > 0)? wid : prec;
2010 pprintf(prn, ifmt, k, (int) x);
2011 } else {
2012 pprintf(prn, ifmt, (int) x);
2013 }
2014 } else {
2015 if (wid > 0 && prec > 0) {
2016 pprintf(prn, xfmt, wid, prec, x);
2017 } else if (wid > 0 || prec > 0) {
2018 k = (wid > 0)? wid : prec;
2019 pprintf(prn, xfmt, k, x);
2020 } else {
2021 pprintf(prn, xfmt, x);
2022 }
2023 }
2024 }
2025 pputc(prn, '\n');
2026 }
2027
2028 free(intcols);
2029 free(ifmt);
2030 free(xfmt);
2031 }
2032 }
2033
2034 /**
2035 * debug_print_matrix:
2036 * @m: matrix to print.
2037 * @msg: accompanying message text (or NULL if no message is wanted).
2038 *
2039 * Prints the matrix @m to stderr with high precision, along
2040 * with the address of the matrix struct.
2041 */
2042
debug_print_matrix(const gretl_matrix * m,const char * msg)2043 void debug_print_matrix (const gretl_matrix *m, const char *msg)
2044 {
2045 char full[64] = {0};
2046
2047 if (msg != NULL) {
2048 strncpy(full, msg, 32);
2049 sprintf(full + strlen(full), " (%p)", (void *) m);
2050 } else {
2051 sprintf(full, " (%p)", (void *) m);
2052 }
2053
2054 if (m != NULL) {
2055 int i, n = m->rows * m->cols;
2056 int d = (int) ceil(log10((double) n));
2057
2058 fprintf(stderr, "%s\n", full);
2059 for (i=0; i<n; i++) {
2060 fprintf(stderr, "val[%0*d] = % .10E\n", d, i, m->val[i]);
2061 }
2062 } else {
2063 PRN *prn;
2064 int err = 0;
2065
2066 prn = gretl_print_new(GRETL_PRINT_STDERR, &err);
2067 if (!err) {
2068 gretl_matrix_print_to_prn(m, full, prn);
2069 gretl_print_destroy(prn);
2070 }
2071 }
2072 }
2073
count_unmasked_elements(const char * mask,int n)2074 static int count_unmasked_elements (const char *mask, int n)
2075 {
2076 int i, c = 0;
2077
2078 for (i=0; i<n; i++) {
2079 c += (mask[i] == 0);
2080 }
2081
2082 return c;
2083 }
2084
2085 /**
2086 * gretl_matrix_cut_rows:
2087 * @m: matrix to process.
2088 * @mask: character array of length equal to the rows of @m,
2089 * with 1s indicating rows to be cut, 0s for rows to be
2090 * retained.
2091 *
2092 * In-place reduction of @m based on @mask: the masked rows
2093 * are cut out of @m.
2094 *
2095 * Returns: 0 on success, non-zero on error.
2096 */
2097
gretl_matrix_cut_rows(gretl_matrix * m,const char * mask)2098 int gretl_matrix_cut_rows (gretl_matrix *m, const char *mask)
2099 {
2100 int i, j, k, n;
2101 double x;
2102
2103 if (m == NULL || mask == NULL) {
2104 return E_DATA;
2105 }
2106
2107 n = count_unmasked_elements(mask, m->rows);
2108
2109 for (j=0; j<m->cols; j++) {
2110 k = 0;
2111 for (i=0; i<m->rows; i++) {
2112 if (!mask[i]) {
2113 x = gretl_matrix_get(m, i, j);
2114 m->val[j * n + k] = x;
2115 k++;
2116 }
2117 }
2118 }
2119
2120 m->rows = n;
2121
2122 return 0;
2123 }
2124
2125 /**
2126 * gretl_matrix_cut_cols:
2127 * @m: matrix to process.
2128 * @mask: character array of length equal to the cols of @m,
2129 * with 1s indicating cols to be cut, 0s for cols to be
2130 * retained.
2131 *
2132 * In-place reduction of @m based on @mask: the masked cols
2133 * are cut out of @m.
2134 *
2135 * Returns: 0 on success, non-zero on error.
2136 */
2137
gretl_matrix_cut_cols(gretl_matrix * m,const char * mask)2138 int gretl_matrix_cut_cols (gretl_matrix *m, const char *mask)
2139 {
2140 int i, j, k, n;
2141 double x;
2142
2143 if (m == NULL || mask == NULL) {
2144 return E_DATA;
2145 }
2146
2147 n = count_unmasked_elements(mask, m->cols);
2148
2149 k = 0;
2150 for (i=0; i<m->cols; i++) {
2151 if (!mask[i]) {
2152 for (j=0; j<m->rows; j++) {
2153 x = gretl_matrix_get(m, j, i);
2154 m->val[k++] = x;
2155 }
2156 }
2157 }
2158
2159 m->cols = n;
2160
2161 return 0;
2162 }
2163
2164 /**
2165 * gretl_matrix_cut_rows_cols:
2166 * @m: square matrix to process.
2167 * @mask: character array of length equal to the dimension
2168 * of @m, with 1s indicating rows and columns to be cut, 0s
2169 * for rows/columns to be retained.
2170 *
2171 * In-place reduction of @m based on @mask: the masked rows
2172 * and columns are cut out of @m. (The data array within @m
2173 * is not "physically" resized.)
2174 *
2175 * Returns: 0 on success, non-zero on error.
2176 */
2177
gretl_matrix_cut_rows_cols(gretl_matrix * m,const char * mask)2178 int gretl_matrix_cut_rows_cols (gretl_matrix *m, const char *mask)
2179 {
2180 gretl_matrix *tmp;
2181 double x;
2182 int i, j, k, l, n;
2183
2184 if (m == NULL || mask == NULL) {
2185 return E_DATA;
2186 } else if (m->rows != m->cols) {
2187 return E_NONCONF;
2188 }
2189
2190 n = count_unmasked_elements(mask, m->rows);
2191 if (n == 0) {
2192 gretl_matrix_reuse(m, 0, 0);
2193 return 0;
2194 }
2195
2196 /* create smaller temporary matrix */
2197 tmp = gretl_matrix_alloc(n, n);
2198 if (tmp == NULL) {
2199 return E_ALLOC;
2200 }
2201
2202 /* copy unmasked values into temp matrix */
2203 k = 0;
2204 for (i=0; i<m->rows; i++) {
2205 if (!mask[i]) {
2206 l = 0;
2207 for (j=0; j<m->cols; j++) {
2208 if (!mask[j]) {
2209 x = gretl_matrix_get(m, i, j);
2210 gretl_matrix_set(tmp, k, l++, x);
2211 }
2212 }
2213 k++;
2214 }
2215 }
2216
2217 /* redimension the original matrix, copy the values back,
2218 and free the temp matrix */
2219 gretl_matrix_reuse(m, n, n);
2220 gretl_matrix_copy_values(m, tmp);
2221 gretl_matrix_free(tmp);
2222
2223 return 0;
2224 }
2225
2226 /**
2227 * gretl_matrix_zero_row_mask:
2228 * @m: matrix to process.
2229 * @err: location to receive error code.
2230 *
2231 * Checks matrix @m for rows that are all zero. If there are
2232 * any such rows, constructs a mask of length equal to the
2233 * number of rows in @m, with 1s indicating zero rows, 0s
2234 * elsewhere. If there are no such rows, returns NULL.
2235 *
2236 * E_ALLOC is written to @err in case a mask should have
2237 * been constructed but allocation failed.
2238 *
2239 * Returns: allocated mask or NULL.
2240 */
2241
gretl_matrix_zero_row_mask(const gretl_matrix * m,int * err)2242 char *gretl_matrix_zero_row_mask (const gretl_matrix *m, int *err)
2243 {
2244 char *mask = NULL;
2245 int row0, any0 = 0;
2246 int i, j;
2247
2248 mask = calloc(m->rows, 1);
2249 if (mask == NULL) {
2250 *err = E_ALLOC;
2251 return NULL;
2252 }
2253
2254 for (i=0; i<m->rows; i++) {
2255 row0 = 1;
2256 for (j=0; j<m->cols; j++) {
2257 if (gretl_matrix_get(m, i, j) != 0.0) {
2258 row0 = 0;
2259 break;
2260 }
2261 }
2262 if (row0) {
2263 mask[i] = 1;
2264 any0 = 1;
2265 }
2266 }
2267
2268 if (!any0) {
2269 free(mask);
2270 mask = NULL;
2271 }
2272
2273 return mask;
2274 }
2275
2276 /**
2277 * gretl_matrix_zero_col_mask:
2278 * @m: matrix to process.
2279 * @err: location to receive error code.
2280 *
2281 * Checks matrix @m for columns that are all zero. If there are
2282 * any such columns, constructs a mask of length equal to the
2283 * number of columns in @m, with 1s indicating zero columns, 0s
2284 * elsewhere. If there are no such columns, returns NULL.
2285 *
2286 * E_ALLOC is written to @err in case a mask should have
2287 * been constructed but allocation failed.
2288 *
2289 * Returns: allocated mask or NULL.
2290 */
2291
gretl_matrix_zero_col_mask(const gretl_matrix * m,int * err)2292 char *gretl_matrix_zero_col_mask (const gretl_matrix *m, int *err)
2293 {
2294 char *mask = NULL;
2295 int col0, any0 = 0;
2296 int i, j;
2297
2298 mask = calloc(m->cols, 1);
2299 if (mask == NULL) {
2300 *err = E_ALLOC;
2301 return NULL;
2302 }
2303
2304 for (j=0; j<m->cols; j++) {
2305 col0 = 1;
2306 for (i=0; i<m->rows; i++) {
2307 if (gretl_matrix_get(m, i, j) != 0.0) {
2308 col0 = 0;
2309 break;
2310 }
2311 }
2312 if (col0) {
2313 mask[j] = 1;
2314 any0 = 1;
2315 }
2316 }
2317
2318 if (!any0) {
2319 free(mask);
2320 mask = NULL;
2321 }
2322
2323 return mask;
2324 }
2325
2326 /**
2327 * gretl_matrix_zero_diag_mask:
2328 * @m: matrix to process.
2329 * @err: location to receive error code.
2330 *
2331 * Checks square matrix @m for diagonal elements that are zero.
2332 * If there are any such, constructs a mask of length equal to
2333 * the number of rows (and columns) of @m, with 1s indicating
2334 * the zero diagonal entries. If there are no zero diagonal
2335 * elements, returns NULL.
2336 *
2337 * E_ALLOC is written to @err in case a mask should have
2338 * been constructed but allocation failed.
2339 *
2340 * Returns: allocated mask or NULL.
2341 */
2342
gretl_matrix_zero_diag_mask(const gretl_matrix * m,int * err)2343 char *gretl_matrix_zero_diag_mask (const gretl_matrix *m, int *err)
2344 {
2345 char *mask = NULL;
2346 int i, trim = 0;
2347
2348 if (gretl_is_null_matrix(m)) {
2349 return NULL;
2350 }
2351
2352 if (m->rows != m->cols) {
2353 *err = E_NONCONF;
2354 return NULL;
2355 }
2356
2357 for (i=0; i<m->rows; i++) {
2358 if (gretl_matrix_get(m, i, i) == 0.0) {
2359 trim = 1;
2360 break;
2361 }
2362 }
2363
2364 if (trim) {
2365 mask = calloc(m->rows, 1);
2366 if (mask == NULL) {
2367 *err = E_ALLOC;
2368 } else {
2369 for (i=0; i<m->rows; i++) {
2370 if (gretl_matrix_get(m, i, i) == 0.0) {
2371 mask[i] = 1;
2372 }
2373 }
2374 }
2375 }
2376
2377 return mask;
2378 }
2379
2380 /**
2381 * gretl_matrix_rank_mask:
2382 * @m: matrix to process.
2383 * @err: location to receive error code.
2384 *
2385 * Performs a QR decomposition of matrix @m and uses this
2386 * to assess the rank of @m. If @m is not of full rank,
2387 * constructs a mask of length equal to the numbers of
2388 * columns in @m, with 1s in positions corresponding
2389 * to diagonal elements of R that are effectively 0, and
2390 * 0s elsewhere. If @m is of full column rank, NULL is
2391 * returned.
2392 *
2393 * E_ALLOC is written to @err in case a mask should have
2394 * been constructed but allocation failed.
2395 *
2396 * Returns: allocated mask or NULL.
2397 */
2398
gretl_matrix_rank_mask(const gretl_matrix * m,int * err)2399 char *gretl_matrix_rank_mask (const gretl_matrix *m, int *err)
2400 {
2401 gretl_matrix *Q = NULL;
2402 gretl_matrix *R = NULL;
2403 char *mask = NULL;
2404 double Ri;
2405 int fullrank = 1;
2406 int i, n = m->cols;
2407
2408 Q = gretl_matrix_copy(m);
2409 if (Q == NULL) {
2410 *err = E_ALLOC;
2411 return NULL;
2412 }
2413
2414 R = gretl_matrix_alloc(n, n);
2415 if (R == NULL) {
2416 gretl_matrix_free(Q);
2417 *err = E_ALLOC;
2418 return NULL;
2419 }
2420
2421 *err = gretl_matrix_QR_decomp(Q, R);
2422
2423 if (!*err) {
2424 mask = calloc(n, 1);
2425 if (mask == NULL) {
2426 *err = E_ALLOC;
2427 }
2428 }
2429
2430 if (!*err) {
2431 for (i=0; i<n; i++) {
2432 Ri = gretl_matrix_get(R, i, i);
2433 if (fabs(Ri) < R_DIAG_MIN) {
2434 mask[i] = 1;
2435 fullrank = 0;
2436 }
2437 }
2438 }
2439
2440 if (*err || fullrank) {
2441 free(mask);
2442 mask = NULL;
2443 }
2444
2445 gretl_matrix_free(Q);
2446 gretl_matrix_free(R);
2447
2448 return mask;
2449 }
2450
2451 /**
2452 * gretl_matrix_mp_ols:
2453 * @y: dependent variable vector.
2454 * @X: matrix of independent variables.
2455 * @b: vector to hold coefficient estimates.
2456 * @vcv: matrix to hold the covariance matrix of the coefficients,
2457 * or NULL if this is not needed.
2458 * @uhat: vector to hold the regression residuals, or NULL if
2459 * these are not needed.
2460 * @s2: pointer to receive residual variance, or NULL. Note:
2461 * if @s2 is NULL, the "vcv" estimate will be plain (X'X)^{-1}.
2462 *
2463 * Computes OLS estimates using Cholesky factorization, via
2464 * the GMP multiple-precision library, and puts the
2465 * coefficient estimates in @b. Optionally, calculates the
2466 * covariance matrix in @vcv and the residuals in @uhat.
2467 *
2468 * Returns: 0 on success, non-zero error code on failure.
2469 */
2470
gretl_matrix_mp_ols(const gretl_vector * y,const gretl_matrix * X,gretl_vector * b,gretl_matrix * vcv,gretl_vector * uhat,double * s2)2471 int gretl_matrix_mp_ols (const gretl_vector *y, const gretl_matrix *X,
2472 gretl_vector *b, gretl_matrix *vcv,
2473 gretl_vector *uhat, double *s2)
2474 {
2475 int (*matrix_mp_ols) (const gretl_vector *,
2476 const gretl_matrix *,
2477 gretl_vector *,
2478 gretl_matrix *,
2479 gretl_vector *,
2480 double *);
2481
2482 matrix_mp_ols = get_plugin_function("matrix_mp_ols");
2483 if (matrix_mp_ols == NULL) {
2484 return 1;
2485 }
2486
2487 return (*matrix_mp_ols)(y, X, b, vcv, uhat, s2);
2488 }
2489
2490 /* following: quadrature sources based on John Burkhart's GPL'd
2491 code at http://people.sc.fsu.edu/~jburkardt
2492 */
2493
2494 #define sign(x) (x < 0.0 ? -1.0 : 1.0)
2495
2496 /* Diagonalize a Jacobi (symmetric tridiagonal) matrix. On entry, @d
2497 contains the diagonal and @e the sub-diagonal. On exit, @d contains
2498 quadrature nodes or abscissae and @z the square roots of the
2499 corresponding weights; @e is used as workspace.
2500
2501 This is a C adaptation of subroutine IMTQLX, by Sylvan Elhay,
2502 Jaroslav Kautsky and John Burkardt. See
2503 http://people.sc.fsu.edu/~jburkardt/f_src/quadrature_weights/
2504 specifically, qw_golub_welsch.f90
2505 */
2506
diag_jacobi(int n,double * d,double * e,double * z)2507 static int diag_jacobi (int n, double *d, double *e, double *z)
2508 {
2509 double b, c, f, g;
2510 double p, r, s;
2511 int j, k, l, m = 0;
2512 int maxiter = 30;
2513 int i, ii, mml;
2514 int err = 0;
2515
2516 if (n == 1) {
2517 return 0;
2518 }
2519
2520 e[n-1] = 0.0;
2521
2522 errno = 0;
2523
2524 for (l=1; l<=n; l++) {
2525 j = 0;
2526 for ( ; ; ) {
2527 for (m=l; m<=n; m++) {
2528 if (m == n) {
2529 break;
2530 }
2531 p = fabs(d[m-1]) + fabs(d[m]);
2532 if (fabs(e[m-1]) <= DBL_EPSILON * p) {
2533 break;
2534 }
2535 }
2536 p = d[l-1];
2537 if (m == l) {
2538 break;
2539 }
2540 if (j >= maxiter) {
2541 fprintf(stderr, "diag_jacobi: iteration limit exceeded\n");
2542 return E_NOCONV;
2543 }
2544 j++;
2545 g = (d[l] - p) / (2.0 * e[l-1]);
2546 r = sqrt(g * g + 1.0);
2547 g = d[m-1] - p + e[l-1] / (g + fabs(r) * sign(g));
2548 s = 1.0;
2549 c = 1.0;
2550 p = 0.0;
2551 mml = m - l;
2552
2553 for (ii=1; ii<=mml; ii++) {
2554 i = m - ii;
2555 f = s * e[i-1];
2556 b = c * e[i-1];
2557 if (fabs(g) <= fabs(f)) {
2558 c = g / f;
2559 r = sqrt(c * c + 1.0);
2560 e[i] = f * r;
2561 s = 1.0 / r;
2562 c = c * s;
2563 } else {
2564 s = f / g;
2565 r = sqrt(s * s + 1.0);
2566 e[i] = g * r;
2567 c = 1.0 / r;
2568 s = s * c;
2569 }
2570 g = d[i] - p;
2571 r = (d[i-1] - g) * s + 2.0 * c * b;
2572 p = s * r;
2573 d[i] = g + p;
2574 g = c * r - b;
2575 f = z[i];
2576 z[i] = s * z[i-1] + c * f;
2577 z[i-1] = c * z[i-1] - s * f;
2578 }
2579
2580 d[l-1] = d[l-1] - p;
2581 e[l-1] = g;
2582 e[m-1] = 0.0;
2583 }
2584 }
2585
2586 /* sorting */
2587
2588 for (ii=2; ii<=m; ii++) {
2589 i = ii - 1;
2590 k = i;
2591 p = d[i-1];
2592 for (j=ii; j<=n; j++) {
2593 if (d[j-1] < p) {
2594 k = j;
2595 p = d[j-1];
2596 }
2597 }
2598 if (k != i) {
2599 d[k-1] = d[i-1];
2600 d[i-1] = p;
2601 p = z[i-1];
2602 z[i-1] = z[k-1];
2603 z[k-1] = p;
2604 }
2605 }
2606
2607 if (errno) {
2608 /* some calculation went badly */
2609 err = E_NAN;
2610 errno = 0;
2611 }
2612
2613 return err;
2614 }
2615
2616 /* Construct the Jacobi matrix for a quadrature rule: @d
2617 holds the diagonal, @e the subdiagonal.
2618 */
2619
make_jacobi(int n,int kind,double * d,double * e)2620 static double make_jacobi (int n, int kind, double *d, double *e)
2621 {
2622 double z0 = 0.0;
2623 int i;
2624
2625 if (kind == QUAD_LEGENDRE) {
2626 double xi, xj;
2627
2628 z0 = 2.0;
2629 for (i=1; i<=n; i++) {
2630 xi = i;
2631 xj = 2.0 * i;
2632 e[i-1] = sqrt(xi * xi / (xj * xj - 1.0));
2633 }
2634 } else if (kind == QUAD_LAGUERRE) {
2635 z0 = 1.0; /* tgamma(1.0) */
2636 for (i=1; i<=n; i++) {
2637 d[i-1] = 2.0 * i - 1.0;
2638 e[i-1] = i;
2639 }
2640 } else if (kind == QUAD_GHERMITE) {
2641 z0 = tgamma(0.5);
2642 for (i=1; i<=n; i++) {
2643 e[i-1] = sqrt(i / 2.0);
2644 }
2645 }
2646
2647 return z0;
2648 }
2649
2650 /* compute basic Gauss quadrature formula */
2651
gauss_quad_basic(int n,int kind,double * x,double * w)2652 static int gauss_quad_basic (int n, int kind, double *x, double *w)
2653 {
2654 double z0, *tmp;
2655 int i, err;
2656
2657 tmp = malloc(n * sizeof *tmp);
2658 if (tmp == NULL) {
2659 return E_ALLOC;
2660 }
2661
2662 z0 = make_jacobi(n, kind, x, tmp);
2663 w[0] = sqrt(z0);
2664 err = diag_jacobi(n, x, tmp, w);
2665
2666 if (!err) {
2667 for (i=0; i<n; i++) {
2668 w[i] = w[i] * w[i];
2669 }
2670 }
2671
2672 free(tmp);
2673
2674 return err;
2675 }
2676
2677 /* Scale Legendre quadrature rule to the interval (a, b) */
2678
legendre_scale(int n,double * x,double * w,double a,double b)2679 static int legendre_scale (int n, double *x, double *w,
2680 double a, double b)
2681 {
2682 double shft, slp;
2683 int i;
2684
2685 if (na(a) || na(b)) {
2686 return E_INVARG;
2687 }
2688
2689 if (fabs(b - a) <= DBL_EPSILON) {
2690 fprintf(stderr, "legendre: |b - a| too small\n");
2691 return E_DATA;
2692 }
2693
2694 shft = (a + b) / 2.0;
2695 slp = (b - a) / 2.0;
2696
2697 for (i=0; i<n; i++) {
2698 x[i] = shft + slp * x[i];
2699 w[i] = w[i] * slp;
2700 }
2701
2702 return 0;
2703 }
2704
hermite_scale(int n,double * x,double * w,double mu,double sigma)2705 static int hermite_scale (int n, double *x, double *w,
2706 double mu, double sigma)
2707 {
2708 double rtpi = sqrt(M_PI);
2709 int i;
2710
2711 if (na(mu) || na(sigma) || sigma <= 0.0) {
2712 return E_INVARG;
2713 }
2714
2715 for (i=0; i<n; i++) {
2716 x[i] = mu + M_SQRT2 * sigma * x[i];
2717 w[i] /= rtpi;
2718 }
2719
2720 return 0;
2721 }
2722
2723 /**
2724 * gretl_quadrule_matrix_new:
2725 * @n: the order (i.e. the number of abscissae and weights).
2726 * @method: should be one of the values in #QuadMethod.
2727 * @a: optional parameter (see below).
2728 * @b: optional parameter (see below).
2729 * @err: location to receive error code.
2730 *
2731 * Calculates a quadrature "rule" (i.e. a set of abscissae or
2732 * nodes and associated weights) for use in numerical integration.
2733 * The three supported methods are Gauss-Hermite, Gauss-Legendre
2734 * and Gauss-Laguerre.
2735 *
2736 * If the optional parameters are not to be used, both should
2737 * be given the value NADBL.
2738 *
2739 * If the method is QUAD_GHERMITE, the default weights are W(x)
2740 * = exp(-x^2), but @a and @b can be used to apply a normal
2741 * scaling with mean @a and standard deviation @b. The limits
2742 * of integration are in principle minus infinity and plus
2743 * infinity.
2744 *
2745 * If the method is QUAD_LEGENDRE the weights are W(x) = 1 and
2746 * the default limits of integration are -1 and 1, but @a and @b
2747 * can be used to adjust the lower and upper limit respectively.
2748 *
2749 * In the QUAD_LAGUERRE case W(x) = exp(-x) and the limits
2750 * of integration are 0 and plus infinity; @a and @b are
2751 * ignored.
2752 *
2753 * Returns: an @n x 2 matrix containing the abscissae in the first
2754 * column and the weights in the second, or NULL on failure.
2755 */
2756
gretl_quadrule_matrix_new(int n,int method,double a,double b,int * err)2757 gretl_matrix *gretl_quadrule_matrix_new (int n, int method,
2758 double a, double b,
2759 int *err)
2760 {
2761 gretl_matrix *m;
2762
2763 if (method < QUAD_GHERMITE || method >= QUAD_INVALID) {
2764 *err = E_DATA;
2765 return NULL;
2766 }
2767
2768 if (n < 0) {
2769 *err = E_DATA;
2770 return NULL;
2771 } else if (n == 0) {
2772 return gretl_null_matrix_new();
2773 }
2774
2775 m = gretl_zero_matrix_new(n, 2);
2776
2777 if (m == NULL) {
2778 *err = E_ALLOC;
2779 } else {
2780 double *x = m->val;
2781 double *w = x + n;
2782
2783 *err = gauss_quad_basic(n, method, x, w);
2784
2785 if (!*err) {
2786 if (method == QUAD_LEGENDRE) {
2787 if (na(a) && na(b)) {
2788 ; /* default */
2789 } else if (a == -1.0 && b == 1.0) {
2790 ; /* default */
2791 } else {
2792 *err = legendre_scale(n, x, w, a, b);
2793 }
2794 } else if (method == QUAD_GHERMITE) {
2795 if (!na(a) || !na(b)) {
2796 *err = hermite_scale(n, x, w, a, b);
2797 }
2798 } else if (method == QUAD_LAGUERRE) {
2799 ; /* a and b are ignored */
2800 }
2801 }
2802 }
2803
2804 if (*err && m != NULL) {
2805 gretl_matrix_free(m);
2806 m = NULL;
2807 }
2808
2809 return m;
2810 }
2811
gretl_gauss_hermite_matrix_new(int n,int * err)2812 gretl_matrix *gretl_gauss_hermite_matrix_new (int n, int *err)
2813 {
2814 return gretl_quadrule_matrix_new(n, QUAD_GHERMITE, 0, 1, err);
2815 }
2816
2817 /**
2818 * gretl_matrix_2d_convolution:
2819 * @A: first matrix.
2820 * @B: second matrix.
2821 * @err: location to receive error code.
2822 *
2823 * Computes the 2-dimensional convolution of the matrices
2824 * @A and @B. The code borrows from Octave's conv2.m.
2825 *
2826 * Returns: newly allocated convolution matrix, or %NULL on
2827 * failure.
2828 */
2829
gretl_matrix_2d_convolution(const gretl_matrix * A,const gretl_matrix * B,int * err)2830 gretl_matrix *gretl_matrix_2d_convolution (const gretl_matrix *A,
2831 const gretl_matrix *B,
2832 int *err)
2833 {
2834 gretl_matrix *C;
2835 double sum;
2836 int ci, cj, bj, aj, bi, ai;
2837 #if 0
2838 /* Octave: "Convolution works fastest if we choose the A matrix to be
2839 * the largest" -- but we're not really seeing any difference.
2840 */
2841 int Am, An, B, Bn, Cm, Cn;
2842 int edgM, edgN;
2843
2844 if (A->rows * A->cols < B->rows * B->cols) {
2845 gretl_matrix *tmp = A;
2846 B = A;
2847 A = tmp;
2848 }
2849
2850 Am = A->rows;
2851 An = A->cols;
2852 Bm = B->rows;
2853 Bn = B->cols;
2854 Cm = Am + Bm - 1;
2855 Cn = An + Bn - 1;
2856 edgM = Bm - 1;
2857 edgN = Bn - 1;
2858 #else
2859 /* leave the input matrices alone */
2860 int Am = A->rows;
2861 int An = A->cols;
2862 int Bm = B->rows;
2863 int Bn = B->cols;
2864 int Cm = Am + Bm - 1;
2865 int Cn = An + Bn - 1;
2866 int edgM = Bm - 1;
2867 int edgN = Bn - 1;
2868 #endif
2869
2870 C = gretl_matrix_alloc(Cm, Cn);
2871
2872 for (ci=0; ci < Cm; ci++) {
2873 for (cj=0; cj < Cn; cj++) {
2874 sum = 0;
2875 for (bj = Bn - 1 - MAX(0, edgN-cj), aj = MAX(0, cj-edgN);
2876 bj >= 0 && aj < An; bj--, aj++) {
2877 bi = Bm - 1 - MAX(0, edgM-ci);
2878 ai = MAX(0, ci-edgM);
2879 for ( ; bi >= 0 && ai < Am; bi--, ai++) {
2880 sum += gretl_matrix_get(A, ai, aj) * gretl_matrix_get(B, bi, bj);
2881 }
2882 }
2883 gretl_matrix_set(C, ci, cj, sum);
2884 }
2885 }
2886
2887 return C;
2888 }
2889
2890 /* scan format string: should be "%<num>m" or just "%m" */
2891
check_matrix_scan_format(const char * fmt,int ns,int * n)2892 static int check_matrix_scan_format (const char *fmt, int ns, int *n)
2893 {
2894 int err = 0;
2895
2896 if (*fmt == '%') {
2897 char *p;
2898 long len = strtol(++fmt, &p, 10);
2899
2900 if (p == fmt) {
2901 /* no "max rows" given */
2902 *n = ns;
2903 } else if (len < 0) {
2904 err = E_INVARG;
2905 } else if (len > ns) {
2906 *n = ns;
2907 } else {
2908 *n = len;
2909 }
2910 if (!err && strcmp(p, "m")) {
2911 err = E_INVARG;
2912 }
2913 } else {
2914 err = E_INVARG;
2915 }
2916
2917 return err;
2918 }
2919
vector_from_strings(char ** S,int ns,const char * fmt,int * nvals,int * err)2920 gretl_matrix *vector_from_strings (char **S, int ns,
2921 const char *fmt,
2922 int *nvals,
2923 int *err)
2924 {
2925 gretl_matrix *m = NULL;
2926 int n = 0;
2927
2928 *nvals = 0;
2929 *err = check_matrix_scan_format(fmt, ns, &n);
2930
2931 if (!*err) {
2932 if (n == 0) {
2933 m = gretl_null_matrix_new();
2934 } else {
2935 m = gretl_column_vector_alloc(n);
2936 }
2937 if (m == NULL) {
2938 *err = E_ALLOC;
2939 }
2940 }
2941
2942 if (!*err && n > 0) {
2943 char *s, *p;
2944 double x;
2945 int i;
2946
2947 gretl_push_c_numeric_locale();
2948 for (i=0; i<n; i++) {
2949 s = S[i];
2950 if (s == NULL || *s == '\0') {
2951 m->val[i] = NADBL;
2952 } else {
2953 x = strtod(s, &p);
2954 if (*p == '\0') {
2955 m->val[i] = x;
2956 *nvals += 1;
2957 } else {
2958 m->val[i] = NADBL;
2959 }
2960 }
2961 }
2962 gretl_pop_c_numeric_locale();
2963 }
2964
2965 return m;
2966 }
2967