extern "C" {
#include <gretl/libgretl.h>
#include <gretl/matrix_extra.h>
}

// perform OLS using matrix functionality 

static int matrix_test_1 (double **Z, const DATAINFO *pdinfo,
			  PRN *prn)
{
    gretl_matrix *X, *y, *b;
    int err = 0;

    X = gretl_matrix_alloc(pdinfo->n, pdinfo->v - 1);
    y = gretl_matrix_alloc(pdinfo->n, 1);
    b = gretl_matrix_alloc(pdinfo->v - 1, 1);

    if (X == NULL || y == NULL || b == NULL) {
	err = E_ALLOC;
    } else {
	int i, t, j = 0;

	for (i=0; i<pdinfo->v; i++) {
	    for (t=0; t<pdinfo->n; t++) {
		if (i == 1) {
		    gretl_vector_set(y, t, Z[i][t]);
		} else {
		    gretl_matrix_set(X, t, j, Z[i][t]);
		}
	    }
	    if (i != 1) j++;
	}

	gretl_matrix_print_to_prn(X, "X", prn);
	gretl_matrix_print_to_prn(y, "y", prn);

	gretl_matrix_ols(y, X, b, NULL, NULL, NULL);
	gretl_matrix_print_to_prn(b, "b", prn);
    }

    gretl_matrix_free(X);
    gretl_matrix_free(y);
    gretl_matrix_free(b);

    return err;
}

// perform OLS via the built-in function lsq() and retrieve
// the coefficient vector

static int matrix_test_2 (double **Z, DATAINFO *pdinfo, PRN *prn)
{
    MODEL m;
    int *list;
    int err = 0;

    list = gretl_list_new(5);
    list[1] = 1;
    list[2] = 0;
    list[3] = 2;
    list[4] = 3;
    list[5] = 4;

    m = lsq(list, Z, pdinfo, OLS, OPT_NONE);

    if (m.errcode != 0) {
	pprintf(prn, "lsq: got error %d\n", m.errcode);
	err = m.errcode;
    } else {
	gretl_matrix *b;

	printmodel(&m, pdinfo, OPT_NONE, prn);
	b = gretl_model_get_matrix(&m, M_COEFF, &err);
	if (err) {
	    pprintf(prn, "error %d getting model matrix\n", err);
	} else {
	    gretl_matrix_print_to_prn(b, "b (from model)", prn);
	    gretl_matrix_free(b);
	}
	clear_model(&m);
    }

    return err;
}

int main (int argc, char **argv)
{
    char *fname;
    double **Z = NULL;
    DATAINFO *pdinfo;
    GretlFileType ftype;
    PRN *prn;
    int err;

    if (argc >= 2) {
        fname = argv[1];
    } else {
	fprintf(stderr, "%s: please give the name of a data file\n",
		argv[0]);
        exit(EXIT_FAILURE);
    }

    libgretl_init();
    prn = gretl_print_new(GRETL_PRINT_STDOUT, NULL);
    pdinfo = datainfo_new();

    ftype = gretl_detect_filetype(fname);

    if (ftype == GRETL_XML_DATA) {
	err = gretl_read_native_data(fname, &Z, pdinfo);
    } else {
	err = gretl_read_foreign_data(fname, ftype, &Z, pdinfo, prn);
    }

    if (err) {
        pprintf(prn, "Got error %d reading data from %s\n", err, fname);
        errmsg(err, prn);
    } else {
        pprintf(prn, "Read data from %s OK\n", fname);
	print_smpl(pdinfo, 0, prn);
        varlist(pdinfo, prn);
	err = matrix_test_1(Z, pdinfo, prn);
	pprintf(prn, "matrix_test_1: err = %d\n", err);
	err = matrix_test_2(Z, pdinfo, prn);
	pprintf(prn, "matrix_test_2: err = %d\n", err);
    }

    destroy_dataset(Z, pdinfo);
    gretl_print_destroy(prn);
    libgretl_cleanup();

    return 0;
}