"Fossies" - the Fresh Open Source Software Archive

Member "gretl-2020b/lib/src/matrix_extra.c" (7 Apr 2020, 62280 Bytes) of package /linux/misc/gretl-2020b.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "matrix_extra.c" see the Fossies "Dox" file reference documentation and the last Fossies "Diffs" side-by-side code changes report: 2020a_vs_2020b.

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