grace_dofit.c

Go to the documentation of this file.

/*
 * Grace - Graphics for Exploratory Data Analysis
 * 
 * Home page: http://plasma-gate.weizmann.ac.il/Grace/
 * 
 * Copyright (c) 1991-95 Paul J Turner, Portland, OR
 * Copyright (c) 1996-98 GRACE Development Team
 * 
 * Maintained by Evgeny Stambulchik <fnevgeny@plasma-gate.weizmann.ac.il>
 * 
 * 
 *                           All Rights Reserved
 * 
 *    This program is free software; you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation; either version 2 of the License, or
 *    (at your option) any later version.
 * 
 *    This program is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 * 
 *    You should have received a copy of the GNU General Public License
 *    along with this program; if not, write to the Free Software
 *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include "ds_extern.h"
#include "ds_proto.h"

/* ALGORITHM AS274.1  APPL. STATIST. (1992) VOL 41, NO. 2

   Calling this routine updates d, rbar, thetab and sserr by the
   inclusion of xrow, yelem with the specified weight.  This version has
   been modified to make it slightly faster when the early elements of
   XROW are not zeroes.  *** WARNING *** The elements of XROW are
   over-written. */
int gf_includ(np, nrbar, w, xrow, y, d, rbar, thetab, sserr)
int np, nrbar;
double w, *xrow, y, *d, *rbar, *thetab, *sserr;
{
  int ier=0, i, k, nextr;
  double cbar, sbar, di, xi, xk, dpi, wxi;


  --thetab; 
  --rbar; 
  --d; 
  --xrow;

    /* the function includ() is called for each obs, and is where all
       the time is spent.  To make it faster, we skip the checks about
       whether nrbar is big enough for np.  If you'd like to have the
       checks in anyway, then compile with -DSTRINGENT */
#ifdef STRINGENT /* Some checks. */
  if(np < 1)
  {
    ier = 1;
  }

  if(nrbar < np * (np - 1) / 2)
  {
     ier += 2;
  }

  if(ier)
  {
    return(ier);
  }

#endif /* STRINGENT */

  nextr = 1;

/* Skip unnecessary transformations.  Test on exact zeroes must be used or stability can be destroyed. */

  for(i=1;i<=np;++i)
  {
    if(w == 0)
    {
      return(ier);
    }

    xi = xrow[i];

    if(xi == 0)
    {
      nextr += np - i;
    }
    else
    {
      di   = d[i];
      wxi  = w * xi;
      dpi  = di + wxi * xi;
      cbar = di / dpi;
      sbar = wxi / dpi;
      w    = cbar * w;
      d[i] = dpi;

      if(i != np)
      {
        for(k=i+1;k<=np;++k)
        {
          xk          = xrow[k];
          xrow[k]     = xk - xi * rbar[nextr];
          rbar[nextr] = cbar * rbar[nextr] + sbar * xk;
          ++nextr;
        }
      }

      xk        = y;
      y         = xk - xi * thetab[i];
      thetab[i] = cbar * thetab[i] + sbar * xk;
    }
  }

/* Y * SQRT(W) is now equal to the Brown, Durbin & Evans recursive residual. */
  *sserr += w * y * y;

  return(ier);
}



/* ALGORITHM AS274.2  APPL. STATIST. (1992) VOL.41, NO. 2
   Sets arrays to zero prior to calling AS75.1 */
int gf_clear(np, nrbar, d, rbar, thetab, sserr)
int np, nrbar;
double *d, *rbar, *thetab, *sserr;
{
  int ier=0;

  if(np < 1)
  {
    ier = 1;
  }

  if(nrbar < (np * (np - 1) / 2))
  {
    ier += 2;
  }

  if(ier)
  {
    return(ier);
  }

  memset(d,      0, np    * sizeof(double));
  memset(thetab, 0, np    * sizeof(double));
  memset(rbar,   0, nrbar * sizeof(double));

  *sserr = 0;

  return(0);
}


int gf_regcf(int np, int nrbar, double *d, double *rbar,
          double *thetab, double *tol, double *beta,
          int nreq)
{
    int ier=0, i, j, nextr;

/* ALGORITHM AS274.3  APPL. STATIST. (1992) VOL 41, NO.2
   Modified version of AS75.4 to calculate regression coefficients
   for the first NREQ variables, given an orthogonal reduction from
   AS75.1. */

    --beta; --tol; --thetab; --rbar; --d;

    /*     Some checks. */
    if (np < 1) ier = 1;
    if (nrbar < np * (np - 1) / 2) ier += 2;
    if (nreq < 1 || nreq > np) ier += 4;
    if (ier) return ier;

    for (i = nreq; i >= 1; --i)
      if (sqrt(d[i]) < tol[i]) {
          beta[i] = d[i] = 0;
      } else {
          beta[i] = thetab[i];
          nextr = (i - 1) * (np + np - i) / 2 + 1;
          for (j=i+1; j<=nreq; ++j) beta[i] -= rbar[nextr++] * beta[j];
      }

    return ier;
} /* regcf */


int gf_tolset(int np, int nrbar, double *d, double *rbar, double *tol)
/* ALGORITHM AS274.4  APPL. STATIST. (1992) VOL 41, NO.2
   Sets up array TOL for testing for zeroes in an orthogonal
   reduction formed using AS75.1.
   EPS is a machine-dependent constant.   For compilers which use
   the IEEE format for floating-point numbers, recommended values
   are 1.E-06 for single precision and 1.E-15 for double precision. */
{
    const double eps = 1e-15;
    double d__1, sum, *work;
    int ier=0, col, pos, row;

    --tol; --rbar; --d;

    /* Some checks. */
    if (np < 1) ier = 1;
    if (nrbar < np * (np - 1) / 2) ier += 2;
    if (ier) return ier;

/* Set TOL(I) = sum of absolute values in column I of RBAR after
   scaling each element by the square root of its row multiplier. */

    work = -1 + (double *) malloc(np*sizeof(double));
    for (col = 1; col <= np; ++col) work[col] = sqrt(d[col]);
    for (col = 1; col <= np; ++col) {
        pos = col - 1;
        sum = work[col];
        for (row = 1; row <= col-1; ++row) {
            sum += (d__1 = rbar[pos], fabs(d__1)) * work[row];
            pos += np - row - 1;
        }
        tol[col] = eps * sum;
    }
    free(work+1);
    return 0;
} /* tolset */


int gf_sing(int np, int nrbar, double *d,
         double *rbar, double *thetab, double *sserr,
         double *tol, int *lindep)
/* ALGORITHM AS274.5  APPL. STATIST. (1992) VOL 41, NO. 2
   Checks for singularities, reports, and adjusts orthogonal
   reductions produced by AS75.1. */
{
    double d__1, temp, *work;
    int ier=0, np2, col, pos, row, pos2;

    --lindep; --tol; --thetab; --rbar; --d;

    if (np <= 0) ier = 1;
    if (nrbar < np * (np - 1) / 2) ier += 2;
    if (ier) return ier;
    work = -1 + (double *) malloc(np*sizeof(double));

/* Set elements within RBAR to zero if they are less than TOL(COL) in
   absolute value after being scaled by the square root of their row
   multiplier. */
    for (col = 1; col <= np; ++col) work[col] = sqrt(d[col]);
    for (col = 1; col <= np; ++col) {
        temp = tol[col];
        pos = col - 1;
        for (row = 1; row < col; ++row) {
            if ((d__1 = rbar[pos], fabs(d__1)) * work[row] < temp)
              rbar[pos] = 0;
            pos += np - row - 1;
        }

/* If diagonal element is near zero, set it to zero, set appropriate
   element of LINDEP, and use INCLUD to augment the projections in
   the lower rows of the orthogonalization. */

        lindep[col] = 0;
        if (work[col] < temp) {
            lindep[col] = 1;
            --ier;
            if (col < np) {
                np2 = np - col;
                pos2 = pos + np2 + 1;
                ier = gf_includ(np2, np2*(np2-1)/2, d[col], rbar+pos+1,
                             thetab[col], d+col+1, rbar+pos2,
                             thetab+col+1, sserr);
            } else {
                d__1 = thetab[col]; /* for squaring */
                *sserr += d[col] * (d__1 * d__1);
            }
            d[col] = work[col] = thetab[col] = 0;
        }
    }
    free(work+1);
    return ier;
} /* sing */


int gf_ss(int np, double *d, double *thetab, double *sserr, double *rss)
/* ALGORITHM AS274.6  APPL. STATIST. (1992) VOL 41, NO. 2
   Calculates partial residual sums of squares from an orthogonal
   reduction from AS75.1. */
{
    double sum;
    int i;

    --rss; --thetab; --d;

    if (np < 1) return 1;
    sum = rss[np] = *sserr;
    for (i = np; i >= 2; --i) {
        sum += d[i] * thetab[i] * thetab[i];
        rss[i - 1] = sum;
    }
    return 0;
} /* ss */


int gf_cov(int np, int nrbar, double *d, 
         double *rbar, int nreq, double *rinv, double *var, 
         double *covmat, int dimcov, double *sterr)
/* ALGORITHM AS274.7  APPL. STATIST. (1992) VOL 41, NO.2 
   Calculate covariance matrix for regression coefficients for the 
   first NREQ variables, from an orthogonal reduction produced from 
   AS75.1. */
{
    int ier=0, k, start, col, pos, row, pos1, pos2;
    double sum;

    --sterr; --covmat; --rinv; --rbar; --d;

    if (np < 1) ier = 1;
    if (nrbar < np * (np - 1) / 2) ier += 2;
    if (dimcov < nreq * (nreq + 1) / 2) ier += 4;

    for (row = 1; row <= nreq; ++row) if (d[row] == 0) ier = -row;
    if (ier) return ier;

    gf_inv(np, nrbar, rbar+1, nreq, rinv+1);
    start = pos = 1;
    for (row = 1; row <= nreq; ++row) {
        pos2 = start;
        for (col = row; col <= nreq; ++col) {
            pos1 = start + col - row;
            sum = (row == col ? 1.0 : rinv[pos1-1])/d[col];
            for (k = col + 1; k <= nreq; ++k)
              sum += rinv[pos1++] * rinv[pos2++] / d[k];
            covmat[pos] = sum * *var;
            if (row == col) sterr[row] = sqrt(covmat[pos]);
            ++pos;
        }
        start += nreq - row;
    }
    return 0;
} /* cov */


void gf_inv(int np, int nrbar, double *rbar, int nreq, double *rinv)
/* ALGORITHM AS274.8  APPL. STATIST. (1992) VOL 41, NO. 2
   Invert first NREQ rows and columns of Cholesky factorization
   produced by AS75.1. */
{
    int pos, row;

    --rinv; --rbar;

    /* Invert RBAR ignoring row multipliers, from the bottom up. */
    pos = nreq * (nreq - 1) / 2;
    for (row = nreq - 1; row >= 1; --row) {
        int col, start = (row - 1) * (np + np - row) / 2 + 1;
        for (col = nreq; col >= row+1; --col) {
            double sum=0;
            int pos1=start, pos2=pos, k;
            for (k = row + 1; k <= col-1; ++k) {
                pos2 += nreq - k;
                sum -= rbar[pos1++] * rinv[pos2];
            }
            rinv[pos] = sum - rbar[pos1];
            --pos;
        }
    }
} /* inv */


int gf_pcorr(int np, int nrbar, double *d, 
          double *rbar, double *thetab, double *sserr, int in, 
          double *cormat, int dimc, double *ycorr)
/* ALGORITHM AS274.9  APPL. STATIST. (1992) VOL 41, NO. 2
   Calculate partial correlations after the first IN variables
   have been forced into the regression. */
{
    int ier=0, i, start, in1;
    double *work;

    --cormat; --thetab; --rbar; --d;

    if (np < 1) ier = 1;
    if (nrbar < np * (np - 1) / 2) ier += 2;
    if (in < 0 || in > np - 1) ier += 4;
    if (dimc < (np - in) * (np - in - 1) / 2) ier += 8;
    if (ier) return ier;

    start = in * (np + np - in - 1) / 2 + 1;
    in1 = in + 1;
    work = (double *) malloc(np*sizeof(double));
    gf_cor(np-in, d+in1, rbar+start, thetab+in1, sserr, work, cormat+1, ycorr);

    /* Check for zeroes. */
    for (i=0; i<(np-in); ++i) if (work[i] <= 0) ier = -1-i;
    free(work);
    return ier;
} /* pcorr */


void gf_cor(int np, double *d, double *rbar, 
         double *thetab, double *sserr, double *work,
         double *cormat, double *ycorr)
/* ALGORITHM AS274.10  APPL. STATIST. (1992) VOL 41, NO.2
   Calculate correlations from an orthogonal reduction.   This
   routine will usually be called from PCORR, which will have
   removed the appropriate number of rows at the start. */
{
    int i__1, i__2, diff, pos, row, col1, col2, pos1, pos2;
    double sumy, sum;

    --ycorr; --cormat; --work; --thetab; --rbar; --d;

/* process by columns, including the projections of the dependent
   variable (THETAB). */
    sumy = *sserr;
    i__1 = np;
    for (row = 1; row <= i__1; ++row)
      sumy += d[row] * thetab[row] * thetab[row];
    sumy = sqrt(sumy);

    pos = np * (np - 1) / 2;
    for (col1 = np; col1 >= 1; --col1) { /* find length of column COL1. */
        sum = d[col1];
        i__1 = pos1 = col1 - 1;
        for (row = 1; row <= i__1; ++row) {
            sum += d[row] * rbar[pos1] * rbar[pos1];
            pos1 += np - row - 1;
        }
        work[col1] = sqrt(sum);

/* If SUM = 0, set all correlations with this variable to zero. */
        if (sum == 0) {
            ycorr[col1] = 0;
            i__1 = col1 + 1;
            for (col2 = np; col2 >= i__1; --col2) {
                cormat[pos] = 0; --pos;
            }
            goto L70;
        }

/* Form cross-products, then divide by product of column lengths. */
        sum = d[col1] * thetab[col1];
        i__1 = pos1 = col1 - 1;
        for (row = 1; row <= i__1; ++row) {
            sum += d[row] * rbar[pos1] * thetab[row];
            pos1 += np - row - 1;
        }
        ycorr[col1] = sum / (sumy * work[col1]);

        i__1 = col1 + 1;
        for (col2 = np; col2 >= i__1; --col2) {
            if (work[col2] > 0) {
                pos1 = col1 - 1;
                pos2 = col2 - 1;
                diff = col2 - col1;
                i__2 = col1 - 1;
                for (sum=0, row = 1; row <= i__2; ++row) {
                    sum += d[row] * rbar[pos1] * rbar[pos2];
                    pos1 += np - row - 1;
                    pos2 = pos1 + diff;
                }
                sum += d[col1] * rbar[pos2];
                cormat[pos] = sum / (work[col1] * work[col2]);
            } else cormat[pos] = 0;
            --pos;
        }
L70:
        ;
    }
} /* cor */




/* ALGORITHM AS274.11 APPL. STATIST. (1992) VOL 41, NO. 2
 * Move variable from position FROM to position TO in an
 * orthogonal reduction produced by AS75.1.
 */
int gf_vmove(np, nrbar, vorder, d, rbar, thetab, rss, from, to, tol)
int np, nrbar, *vorder;
double *d, *rbar, *thetab, *rss;
int from, to;
double *tol;
{
  int ier=0, last, m;
  double d__1, cbar, sbar, d1new, d2new, x, y, d1, d2;
  int first, m1, m2, mp1, inc, col, pos, row;

  --tol;
  --rss;
  --thetab;
  --rbar;
  --d;
  --vorder;

  if(np <= 0)
  {
    ier = 1;
  }

  if(nrbar < np * (np - 1) / 2)
  {
    ier += 2;
  }

  if(from < 1 || from > np)
  {
    ier += 4;
  }

  if(to < 1 || to > np)
  {
    ier += 8;
  }

  if(ier)
  {
    return ier;
  }

  if(from == to)
  {
    return ier;
  }

  if(from < to)
  {
    first = from;
    last  = to - 1;
    inc   = 1;
  }
  else
  {
    first = from - 1;
    last  = to;
    inc   = -1;
  }

/* Find addresses of first elements of RBAR in rows M and (M+1). */    
  for (m = first; inc < 0 ? m >= last : m <= last; m += inc)
  {
    m1  = (m - 1) * (np + np - m) / 2 + 1;
    m2  = m1 + np - m;
    mp1 = m + 1;
    d1  = d[m];
    d2  = d[mp1];

        /* Special cases. */
    if(d1 == 0 && d2 == 0)
    {
      goto L40;
    }

    x = rbar[m1];
    if(fabs(x) * sqrt(d1) < tol[mp1])
    {
      x = 0;
    }

    if(d1 == 0 || x == 0)
    {
      d[m] = d2;
      d[mp1] = d1;
      rbar[m1] = 0;

      for (col = m + 2; col <= np; ++col)
      {
        x          = rbar[++m1];
        rbar[m1]   = rbar[m2];
        rbar[m2++] = x;
      }

      x           = thetab[m];
      thetab[m]   = thetab[mp1];
      thetab[mp1] = x;

      goto L40;
    }
    else
    if(d2 == 0)
    {
      d__1 = x;
      d[m] = d1 * d__1 * d__1;
      rbar[m1] = 1.0/x;

      for(col = m + 2; col <= np; ++col)
      {
        rbar[++m1] /= x;
        thetab[m] /= x;
      }

      goto L40;
    }

/*  Planar rotation in regular case. */
    d__1     = x;
    d1new    = d2 + d1 * (d__1 * d__1);
    cbar     = d2 / d1new;
    sbar     = x * d1 / d1new;
    d2new    = d1 * cbar;
    d[m]     = d1new;
    d[mp1]   = d2new;
    rbar[m1] = sbar;

    for(col = m + 2; col <= np; ++col)
    {
      y        = rbar[++m1];
      rbar[m1] = cbar * rbar[m2] + sbar * y;
      rbar[m2] = y - x * rbar[m2];
      ++m2;
    }

    y           = thetab[m];
    thetab[m]   = cbar * thetab[mp1] + sbar * y;
    thetab[mp1] = y - x * thetab[mp1];

/*     Swap columns M and (M+1) down to row (M-1). */
L40:
    if(m == 1)
    {
      goto L60;
    }

    pos = m;

    for(row = 1; row < m; ++row)
    {
      x             = rbar[pos];
      rbar[pos]     = rbar[pos - 1];
      rbar[pos - 1] = x;
      pos          += np - row - 1;
    }

/* Adjust variable order (VORDER), the tolerances (TOL) and the vector
   of residual sums of squares (RSS). */
L60:
    m1          = vorder[m];
    vorder[m]   = vorder[mp1];
    vorder[mp1] = m1;
    x           = tol[m];
    tol[m]      = tol[mp1];
    tol[mp1]    = x;
    d__1        = thetab[mp1];
    rss[m]      = rss[mp1] + d[mp1] * (d__1 * d__1);
  }

  return(0);
} /* gf_vmove */


/* ALGORITHM AS274.12 APPL. STATIST. (1992) VOL 41, NO.2 Re-order the
 * variables in an orthogonal reduction produced by AS75.1 so that the N
 * variables in LIST start at position POS1, though will not necessarily
 * be in the same order as in LIST.  Any variables in VORDER before
 * position POS1 are not moved.  Auxiliary routine called: VMOVE
 */
int gf_reordr(np, nrbar, vorder, d, rbar, thetab, rss, tol, list, n, pos1)
int np, nrbar, *vorder;
double *d, *rbar, *thetab, *rss, *tol;
int *list, n, pos1;
{
  int ier=0, next, i, j;

  --list;
  --tol;
  --rss;
  --thetab;
  --rbar;
  --d;
  --vorder;

  if(np < 1)
  {
    ier = 1;
  }

  if(nrbar < np * (np - 1) / 2)
  {
    ier += 2;
  }

  if(n < 1 || n >= np + 1 - pos1)
  {
    ier += 4;
  }

  if(ier)
  {
    return(ier);
  }

/* Work through VORDER finding variables which are in LIST. */
  i = next = pos1;
  do
  {
    int l=vorder[i];

    for (j=1; j<=n; j++)
    {
      if (l==list[j])
      {
        if (i > next) 
        {
          gf_vmove(np, nrbar, vorder+1, d+1, rbar+1, thetab+1, rss+1, i, next, tol+1);
        }

        if(++next >= n + pos1)
        {
          return 0;
        }

        break;
      }
    }
  }while(++i <= np);

  return(8);
} /* gf_reordr */



/* ALGORITHM AS274.13  APPL. STATIST. (1992) VOL.41, NO.2 */     
int gf_hdiag(double *xrow, int np, int nrbar, double *d, double *rbar, double *tol, int nreq, double *hii)
{
  int i__1, i__2, ier=0, col, pos, row;
  double d__1, sum, *wk;

  --tol;
  --rbar;
  --d;
  --xrow;

  if(np < 1)
  {
    ier = 1;
  }

  if(nrbar < np * (np - 1) / 2)
  {
    ier += 2;
  }

  if(nreq > np)
  {
    ier += 4;
  }

  if(ier)
  {
    return ier;
  }


/* The elements of XROW.inv(RBAR).sqrt(D) are calculated and stored in WK. */
  wk = -1 + (double *) malloc(np*sizeof(double));
  *hii = 0;
  i__1 = nreq;

  for(col=1;col<= i__1;++col)
  {
    if(sqrt(d[col]) <= tol[col])
    {
      wk[col] = 0;
    }
    else
    {
      i__2 = pos = col - 1;
      sum = xrow[col];

      for(row = 1; row <= i__2; ++row)
      {
        sum -= wk[row] * rbar[pos];
        pos += np - row - 1;
      }

      d__1 = wk[col] = sum;
      *hii += d__1 * d__1 / d[col];
    }
  }

  free(wk+1);

  return(0);
} /* gf_hdiag */




/* "print a upper triangular double matrix stored as a vector"
 *  The matrix is N x N, the vector has N*(N+1)/2 elements.
 *  Each element is formatted using width and precision.
 *  There are no sanity checks at all.
 */
void gf_pr_utdm_v(double *x, int N, int width, int precision)
{
  int pos=0, i, j, leavespace;
  char s[100], fmt[100];      /* will be used in making printf() formats */
  char buf[256];

  sprintf(fmt, "%%%d.%dg", width, precision);

  for(i=0;i<N;i++)
  {
    leavespace = i*width;
    sprintf(s, "%%%ds", leavespace); 
    printf(buf, s, "");

    for(j=i;j<N;j++)
    {
      printf(buf, fmt, x[pos++]);
    }

    printf("\n");
  }
}



double *gf_dvector(int l, int h)
{
  double *block;
  int size;

  size = h - l + 1;
  block = (double *) malloc(sizeof(double)*size);

  if(block == NULL)
  {
    fprintf(stderr, "malloc failure in dvector()\n");

    return NULL;
  }

  return(block-l);
}



int *gf_ivector(int l, int h)
{
  int *block, size;

  size  = h - l + 1;
  block = (int *) malloc(sizeof(int)*size);

  if(block == NULL)
  {
    fprintf(stderr, "malloc failure in ivector()\n");

    return NULL;
  }

  return(block-l);
}



double **gf_dmatrix(int rl, int rh, int cl, int ch)
{
  double *block;
  double **m;
  int size, i, rowsize, numrows;

  rowsize = ch - cl + 1;        /* #locations consumed by 1 row */
  numrows = rh -rl + 1;
  size    = numrows*rowsize;
  block   = (double *) malloc((unsigned) sizeof(double)*size);

  if(block == NULL)
  {
    fprintf(stderr, "malloc failure in matrix allocation\n");

    return(NULL);
  }
    /* so we have the matrix. */

  /* Now for the row pointers */
  m = (double **) malloc((unsigned) sizeof(double *)*numrows);

  if(m == NULL)
  {
    fprintf(stderr, "malloc failure in matrix allocation\n");

    return(NULL);
  }
  m -= rl;      /* fixup m pointer so m[rl] == old m[0] */

    /* Finally, setup pointers to rows */
  block -= cl;

  for(i=rl;i<=rh;i++)
  {
    m[i] = block; block += rowsize;
  }

  return(m);
}



void gf_putdvec(const char *s, double *x, int l, int h)
{
  int i;
  char buf[512];

  printf(buf, "Vector %-10s: \n", s);

  for(i=l;i<=h;i++)
  {
    printf(buf, " %d: %.4g \n", i, x[i]);
  }
}


int dofitcurve(int cnt, double *xd, double *yd, int nd, double *c)
{
  double **x, *y, *d, *rbar, *thetab, *xrow, sserr, *beta;
  double *tol, *rss;
  int i, j, nvars = nd + 1, nrbar = 200, *vorder, *lindep, error;
  double xval;
  char buf[256];
  extern int debug;

  x      = gf_dmatrix(0, cnt - 1, 0, nvars - 1);
  y      = gf_dvector(0, cnt - 1);
  d      = gf_dvector(0, nvars - 1);
  rbar   = gf_dvector(0, nrbar);
  thetab = gf_dvector(0, nvars - 1);
  xrow   = gf_dvector(0, nvars - 1);
  beta   = gf_dvector(0, nvars - 1);
  tol    = gf_dvector(0, nvars - 1);
  rss    = gf_dvector(0, nvars - 1);
  vorder = gf_ivector(0, nvars - 1);
  lindep = gf_ivector(0, nvars - 1);

  for(i=0;i<nvars;i++)
  {
    vorder[i] = i;
  }

  for(i=0;i<cnt;i++)
  {
    xval = xd[i];

    for (j=1;j<nvars;j++)
    {
      x[i][j - 1] = pow(xval, (double) j);
    }

    y[i] = yd[i];
  }

  error = gf_clear(nvars, nrbar, d, rbar, thetab, &sserr);

  if(error)
  {
    printf("as274c: clear() returned %d\n", error);
    goto bustout;
  }

  for (i=0; i<cnt; i++)
  {
    xrow[0] = 1; /* include constant */

    for(j=0;j<nvars-1;j++)
    {
      xrow[j+1] = x[i][j];
    }

    error = gf_includ(nvars, nrbar, 1.0, xrow, y[i], d, rbar, thetab, &sserr);

    if(error)
    {
      printf("as274c: includ() returned %d at row %d\n", error, i);
      goto bustout;
    }
  }

  error = gf_tolset(nvars, nrbar, d, rbar, tol);

  if(error)
  {
    printf("as274c: tolset() returned %d\n", error);
    goto bustout;
  }

  error = gf_sing(nvars, nrbar, d, rbar, thetab, &sserr, tol, lindep);

  if(error)
  {
    printf("as274c: sing() returned %d\n", error);
    goto bustout;
  }

  if(debug)
  {
    printf(buf, "\nSSerr = %17g\n", sserr);
  }

  error = gf_regcf(nvars, nrbar, d, rbar, thetab, tol, beta, nvars);

  if(error)
  {
    printf("as274c: regcf() returned %d\n", error);
    goto bustout;
  }

  if(debug)
  {
    printf("\nVariable order:\n ");

    for(j=0;j<nvars;j++)
    {
      printf(buf, "   %d ", vorder[j]);
    }

    printf("\n");

    gf_putdvec("Beta", beta, 0, nvars - 1);
    gf_putdvec("d", d, 0, nvars - 1);

    printf(buf, "rbar matrix:\n"); 

    gf_pr_utdm_v(rbar, nvars - 1, 13, 5);
    gf_putdvec("thetab", thetab, 0, nvars - 1);
  }

  for(j=0;j<nvars;j++)
  {
    c[j] = beta[j];
  }

bustout:;
  free(x[0]);
  free(x);
  free(y);
  free(d);
  free(rbar);
  free(thetab);
  free(xrow);
  free(beta);
  free(tol);
  free(rss);
  free(vorder);
  free(lindep);

  return(0);
}