/* lad.c -- Least Absolute Deviation regression for gretl */

#include "libgretl.h"

static double toler = 1.0e-9;

static void col_ (double *v1, double *v2, double amlt, 
		  int m1, int iout)
{
    int i;

    /* Parameter adjustments */
    --v2;
    --v1;

    for (i=1; i<=m1; i++) {
	if (i != iout) {
	    v1[i] -= v2[i] * amlt;
	}
    }
}

/*    Based on SUBROUTINE L1(M,N,TOLER,X,A,B)
C
C  ***************************************************************
C  *  COMMUNICATIONS OF THE ASSOCIATION FOR COMPUTING MACHINERY  *
C  *  ALGORITHM 478                                              *
C  *  SOLUTION OF AN OVERDETERMINED SYSTEM OF EQUATIONS IN THE   *
C  *  L1 NORM -- I. BARRODALE AND F.D.K. ROBERTS                 *
C  ***************************************************************
C
C  ****************************************************************
C  * THIS SUBROUTINE USES A MODIFICATION OF THE SIMPLEX METHOD OF *
C  * LINEAR PROGRAMMING TO CALCULATE AN L1 SOLUTION TO THE OVER-  *
C  * DETERMINED SYSTEM OF LINEAR EQUATIONS, A*X = B.
C  *                                                              *
C  *                         PARAMETERS                           *
C  *                                                              *
C  * M     - NUMBER OF EQUATIONS.  NOTE THAT THIS IS THE N IN     *
C  *         SUBROUTINE LAD.                                      *
C  * N     - NUMBER OF UNKNOWNS (M.GE.N).                         *
C  * TOLER - A SMALL POSITIVE TOLERANCE. THE ROUTINE REGARDS ANY  *
C  *         QUANTITY AS ZERO UNLES ITS MAGNITUDE EXCEEDS TOLER.  *
C  *         EMPIRICAL EVIDENCE SUGGESTS TOLER = 10**(-D*2/3)     *
C  *         WHERE D REPRESENTS THE NUMBER OF DECIMAL DIGITS OF   *
C  *         ACCURACY AVAILABLE.                                  *
C  * X     - ONE DIMENSIONAL REAL ARRAY.  ON EXIT, THIS ARRAY     *
C  *         CONTAINS A SOLUTION TO THE L1 PROBLEM.               *
C  *                                                              *
C  * THE ORIGINAL CONTENTS OF THE ARRAYS A AND B ARE DESTROYED BY *
C  * THIS ROUTINE.                                                *
C  *                                                              *
C  * ON EXIT FROM THE SUBROUTINE, THE ARRAY A CONTAINS THE        *
C  * FOLLOWING INFORMATION.                                       *
C  *                                                              *
C  * A(M+1,N+1)  THE MINIMUM SUM OF THE ABSOLUTE VALUES OF THE    *
C  *             RESIDUALS.                                       *
C  * A(M+1,N+2)  THE RANK OF THE MATRIX OF COEFFICIENTS.          *
C  * A(M+2,N+1)  EXIT CODE WITH VALUES                            *
C  *             0 - OPTIMAL SOLUTION - PROBABLY NON-UNIQUE       *
C  *             1 - UNIQUE OPTIMAL SOLUTION                      *
C  *             2 - PREMATURE TERMINATION DUE TO ROUNDING ERRORS *
C  * A(M+2,N+2)  NUMBER OF SIMPLEX ITERATIONS PERFORMED.          *
C  *                                                              *
C  ****************************************************************
C */

/* Modifications to the above: "toler" is not entered as a parameter,
   it is "global" to this translation unit.  Regression residuals
   are found in the array e on exit.  The following function was
   translated from the Fortran by f2c, then rendered into slightly more
   idiomatic C by me.  Allin Cottrell, September 2002.
*/

static int l1_ (int m, int n, 
		double *a, double *b, 
		double *x, double *e)
{
    const double big = 1e200;
    double amin, amax;
    int test = 0, stage;
    double sum;
    int iout = 0;
    int i, j, k, l;
    int m1, n1, m2, n2, kount;
    double d, pivot;
    int kl, in = 0, kr, *ls;
    int nrows = m + 2;

    ls = malloc(m * sizeof *ls);

    /* Parameter adjustments */
    --e;
    --x;
    --b;
    a -= (nrows + 1);

    m1 = m + 1;
    n1 = n + 1;
    m2 = m + 2;
    n2 = n + 2;

    for (j = 1; j <= n; ++j) {
	a[m2 + j * nrows] = (double) j;
	x[j] = 0.0;
    }

    for (i = 1; i <= m; ++i) {
	a[i + n2 * nrows] = (double) (n + i);
	a[i + n1 * nrows] = b[i];
	if (b[i] < 0.0) {
	    for (j = 1; j <= n2; ++j) {
		a[i + j * nrows] = -a[i + j * nrows];
	    }
	}
	e[i] = 0.0;
    }

    for (j = 1; j <= n1; ++j) {
	sum = 0.0;
	for (i = 1; i <= m; ++i) {
	    sum += a[i + j * nrows];
	}
	a[m1 + j * nrows] = sum;
    }

    stage = 1;
    kount = 0;
    kr = 1;
    kl = 1;

L70:
    amax = -1.0;
    for (j = kr; j <= n; ++j) {
	if (fabs(a[m2 + j * nrows]) > (double) n) {
	    continue;
	}
	d = fabs(a[m1 + j * nrows]);
	if (d <= amax) {
	    continue;
	}
	amax = d;
	in = j;
    }

    if (a[m1 + in * nrows] < 0.0) {
	for (i = 1; i <= m2; ++i) {
	    a[i + in * nrows] = -a[i + in * nrows];
	}
    }

L100:
    k = 0;
    for (i = kl; i <= m; ++i) {
	d = a[i + in * nrows];
	if (d > toler) {
	    ++k;
	    b[k] = a[i + n1 * nrows] / d;
	    ls[k - 1] = i;
	    test = 1;
	}
    }

L120:
    if (k <= 0) {
	test = 0;
	goto L150;
    }

    amin = big;
    for (i = 1; i <= k; ++i) {
	if (b[i] < amin) {
	    j = i;
	    amin = b[i];
	    iout = ls[i - 1];
	}
    }
    b[j] = b[k];
    ls[j - 1] = ls[k - 1];
    --k;

L150:
    if (test || !stage) {
	goto L170;
    }
    for (i = 1; i <= m2; ++i) {
	d = a[i + kr * nrows];
	a[i + kr * nrows] = a[i + in * nrows];
	a[i + in * nrows] = d;
    }
    ++kr;
    goto L260;

L170:
    if (test) {
	goto L180;
    }
    a[m2 + n1 * nrows] = 2.0;
    goto L350;

L180:
    pivot = a[iout + in * nrows];
    if (a[m1 + in * nrows] - pivot - pivot <= toler) {
	goto L200;
    }
    for (j = kr; j <= n1; ++j) {
	d = a[iout + j * nrows];
	a[m1 + j * nrows] = a[m1 + j * nrows] - d - d;
	a[iout + j * nrows] = -d;
    }
    a[iout + n2 * nrows] = -a[iout + n2 * nrows];
    goto L120;

L200:
    for (j = kr; j <= n1; ++j) {
	if (j != in) {
	    a[iout + j * nrows] /= pivot;
	}
    }
    for (j = kr; j <= n1; ++j) {
	if (j != in) {
	    col_(&a[j * nrows + 1], &a[in * nrows + 1], a[iout + j * nrows], 
		 m1, iout);
	}
    }
    for (i = 1; i <= m1; ++i) {
	if (i != iout) {
	    a[i + in * nrows] = -a[i + in * nrows] / pivot;
	}
    }

    a[iout + in * nrows] = 1.0 / pivot;
    d = a[iout + n2 * nrows];
    a[iout + n2 * nrows] = a[m2 + in * nrows];
    a[m2 + in * nrows] = d;
    ++kount;

    if (!stage) {
	goto L270;
    }

    ++kl;
    for (j = kr; j <= n2; ++j) {
	d = a[iout + j * nrows];
	a[iout + j * nrows] = a[kount + j * nrows];
	a[kount + j * nrows] = d;
    }

L260:
    if (kount + kr != n1) {
	goto L70;
    }

    stage = 0;

L270:
    amax = -big;
    for (j = kr; j <= n; ++j) {
	d = a[m1 + j * nrows];
	if (d < 0.0) {
	    if (d > -2.0) {
		continue;
	    }
	    d = -d - 2.0;
	}
	if (d <= amax) {
	    continue;
	}
	amax = d;
	in = j;
    }

    if (amax > toler) {
	if (a[m1 + in * nrows] > 0.0) {
	    goto L100;
	}
	for (i = 1; i <= m2; ++i) {
	    a[i + in * nrows] = -a[i + in * nrows];
	}
	a[m1 + in * nrows] -= 2.0;
	goto L100;
    }

    l = kl - 1;
    for (i = 1; i <= l; ++i) {
	if (a[i + n1 * nrows] < 0.0) {
	    for (j = kr; j <= n2; ++j) {
		a[i + j * nrows] = -a[i + j * nrows];
	    }
	}
    }

    a[m2 + n1 * nrows] = 0.0;
    if (kr != 1) {
	goto L350;
    }

    for (j = 1; j <= n; ++j) {
	d = fabs(a[m1 + j * nrows]);
	if (d <= toler || 2.0 - d <= toler) {
	    goto L350;
	}
    }
    a[m2 + n1 * nrows] = 1.0;

L350:
    for (i = 1; i <= m; ++i) {
	k = (int) a[i + n2 * nrows];
	d = a[i + n1 * nrows];
	if (k <= 0) {
	    k = -k;
	    d = -d;
	}
	if (i >= kl) {
	    k -= n;
	    e[k] = d;
	} else {
	    x[k] = d;
	}
    }

    a[m2 + n2 * nrows] = (double) kount;
    a[m1 + n2 * nrows] = (double) (n1 - kr);

    sum = 0.0;
    for (i = kl; i <= m; ++i) {
	sum += a[i + n1 * nrows];
    }

    a[m1 + n1 * nrows] = sum;

    free(ls);

    return 0;
} 

static int missobs_before (const MODEL *pmod, int t)
{
    int i, c = 0;

    for (i=0; i<t; i++) {
	if (pmod->missmask[i + pmod->t1] == '1') {
	    c++;
	}
    }

    return c;
}

static void 
adjust_sample_for_missing (int *sample, int n, const MODEL *pmod)
{
    int i;

    for (i=0; i<n; i++) {
	sample[i] += missobs_before(pmod, sample[i]);
    }
}

#define ITERS 500

/* populate dependent var using resampled residuals */

static void
make_data_arrays (MODEL *pmod, double **Z,
		  double *a, double *b,
		  const int *sample,
		  int nrows, int k, int m)
{
    int i, j, v, t;

    /* we need to do this on each iteration because the "a" array is
       overwritten by the LAD calculations
    */
    for (j=0; j<k; j++) {
	v = pmod->list[j+2];
	t = pmod->t1;
	for (i=0; i<m; i++) {
	    while (model_missing(pmod, t)) {
		t++;
	    }
	    a[i + j * nrows] = Z[v][t++];
	}
    }

    t = pmod->t1;

    for (i=0; i<m; i++) {
	while (model_missing(pmod, t)) {
	    t++;
	}
	b[i] = a[i + k * nrows] = 
	    pmod->yhat[t++] + pmod->uhat[sample[i]];
    }
}

/* obtain bootstrap estimates of LAD covariance matrix */

static int bootstrap_vcv (MODEL *pmod, double **Z,
			  double *a, double *b, double *e, double *x,
			  int m, int n, int dim)
{
    double **coeffs = NULL;
    double *meanb = NULL;
    int *sample = NULL;
    double xi, xj;
    int i, j, k;
    int nvcv, nrows = m + 2;
    int err = 0;

    /* note: new_vcv sets all entries to zero */
    err = gretl_model_new_vcv(pmod, &nvcv);
    if (err) {
	return err;
    }

    /* an array for each coefficient */
    coeffs = malloc(pmod->ncoeff * sizeof *coeffs);
    if (coeffs == NULL) {
	err = E_ALLOC;
	goto bailout;
    }

    /* a scalar for each coefficient mean */
    meanb = malloc(pmod->ncoeff * sizeof *meanb);
    if (meanb == NULL) {
	err = E_ALLOC;
	goto bailout;
    }    

    /* each array has length ITERS */
    for (i=0; i<pmod->ncoeff; i++) {
	coeffs[i] = malloc(ITERS * sizeof **coeffs);
	if (coeffs[i] == NULL) {
	    for (k=0; k<i; k++) {
		free(coeffs[k]);
	    }
	    free(coeffs);
	    coeffs = NULL;
	    err = E_ALLOC;
	    goto bailout;
	}
    }

    /* sample array has length pmod->nobs */
    sample = malloc(m * sizeof *sample);
    if (sample == NULL) {
	err = E_ALLOC;
	goto bailout;
    }

    for (k=0; k<ITERS; k++) {

	/* create random sample index array */
	for (i=0; i<m; i++) {
	    sample[i] = pmod->t1 + gretl_rand_int_max(m);
	}

	if (pmod->missmask != NULL) {
	    adjust_sample_for_missing(sample, m, pmod);
	}

	/* initialize arrays */
	for (i=0; i<dim; i++) {
	    a[i] = 0.0;
	}
	for (i=0; i<m; i++) {
	    e[i] = b[i] = 0.0;
	}
	for (i=0; i<n; i++) {
	    x[i] = 0.0;
	}

	make_data_arrays(pmod, Z, a, b, sample, nrows, n, m);

	/* estimate LAD model and store coeffs */
	l1_(m, n, a, b, x, e);

	for (i=0; i<n; i++) {
	    coeffs[i][k] = x[i];
	}
    }

    /* find means of coeff estimates */
    for (i=0; i<pmod->ncoeff; i++) {
	double bbar = 0.0;

	for (k=0; k<ITERS; k++) {
	   bbar += coeffs[i][k];
	} 
	meanb[i] = bbar / ITERS;
    }    

    /* find variances and covariances */
    for (i=0; i<pmod->ncoeff; i++) {
	double vi = 0.0;

	for (k=0; k<ITERS; k++) {
	    xi = coeffs[i][k] - meanb[i];
	    vi += xi * xi;
	    for (j=0; j<=i; j++) {
		xj = coeffs[j][k] - meanb[j];
		pmod->vcv[ijton(i, j, pmod->ncoeff)] += xi * xj;
	    }
	}
	pmod->sderr[i] = sqrt(vi / ITERS);
    }

    for (i=0; i<nvcv; i++) {
	pmod->vcv[i] /= ITERS;
    }

 bailout:

    free(sample);

    if (coeffs != NULL) {
	for (i=0; i<pmod->ncoeff; i++) {
	    free(coeffs[i]);
	}
	free(coeffs);
    }

    free(meanb);

    return 0;
}

static int record_y_median (MODEL *pmod, const double *y)
{
    int T = pmod->t2 - pmod->t1 + 1;
    double *sy, m;
    int t, n, n2p;

    sy = malloc(T * sizeof *sy);

    if (sy == NULL) {
	return E_ALLOC;
    }

    n = 0;
    for (t=pmod->t1; t<=pmod->t2; t++) {
	if (!model_missing(pmod, t)) {
	    sy[n++] = y[t];
	}
    }

    if (n == 0) {
	free(sy);
	return E_DATA;
    }

    qsort(sy, n, sizeof *sy, gretl_compare_doubles); 

    n2p = (T = n / 2) + 1;
    m = (n % 2)? sy[n2p - 1] : 0.5 * (sy[T - 1] + sy[n2p - 1]);

    gretl_model_set_double(pmod, "ymedian", m);

    free(sy);

    return 0;
}

int lad_driver (MODEL *pmod, double **Z, DATAINFO *pdinfo)
{
    double *a = NULL, *b = NULL, *e = NULL, *x = NULL;
    int i, j, t, m, n, nrows, dim;
    int yno = pmod->list[1];
    int ladcode;

    m = pmod->nobs;
    n = pmod->list[0] - 1;

    nrows = m + 2;
    dim = nrows * (n + 2);

    a = malloc(dim * sizeof *a);
    x = malloc(n * sizeof *x);
    e = malloc(m * sizeof *e);
    b = malloc(m * sizeof *b);

    if (a == NULL || x == NULL || e == NULL || b == NULL) {
	free(a);
	free(x);
	free(e);
	free(b);
	return 1;
    }

    /* initialize arrays */
    for (i=0; i<dim; i++) {
	a[i] = 0.0;
    }
    for (i=0; i<m; i++) {
	e[i] = b[i] = 0.0;
    }
    for (i=0; i<n; i++) {
	x[i] = 0.0;
    }

    /* populate data array */
    for (j=0; j<n; j++) {
	int v = pmod->list[j+2];

	t = pmod->t1;
	for (i=0; i<m; i++) {
	    while (model_missing(pmod, t)) {
		t++;
	    }
	    a[i + j * nrows] = Z[v][t++];
	}
    }

    t = pmod->t1;
    for (i=0; i<m; i++) {
	while (model_missing(pmod, t)) {
	    t++;
	}
	b[i] = a[i + n * nrows] = Z[yno][t++];
    }

    l1_(m, n, a, b, x, e);

    /* handle case where exit code indicates numeric error */
    ladcode = (int) a[m + 1 + n * nrows];
    if (ladcode == 2) {
	pmod->errcode = E_SINGULAR;
    } else {
	gretl_model_set_int(pmod, "ladcode", ladcode);
    }

    if (pmod->errcode == 0) {

	for (i=0; i<n; i++) {
	    pmod->coeff[i] = x[i];
	}

	pmod->ess = 0.0;
	for (i=0; i<m; i++) {
	    t = i + pmod->t1;
	    pmod->yhat[t] = Z[yno][t] - e[i];
	    pmod->uhat[t] = e[i];
	    pmod->ess += e[i] * e[i];
	}

	/* sum of absolute residuals */
	gretl_model_set_double(pmod, "ladsum", a[m + n * nrows]);

	/* median of dependent variable */
	record_y_median(pmod, Z[yno]);

	/* set ess-based stats to missing value */
	pmod->rsq = NADBL;
	pmod->adjrsq = NADBL;
	pmod->fstt = NADBL;

	/* LaPlace errors: equivalent of standard error is sum of
	   absolute residuals over nobs */
	pmod->sigma = pmod->rho / pmod->nobs; 

	if (bootstrap_vcv(pmod, Z, a, b, e, x, m, n, dim)) {
	    pmod->errcode = E_ALLOC;
	}

    }

    pmod->ci = LAD;

    free(a);
    free(x);
    free(e);
    free(b);

    return 0;
}