/* nag_regsn_mult_linear(g02dac) Example Program
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5 revised, 1998.
 * Mark 6 revised, 2000.
 */

#include <nag.h>
#include <math.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nagg02.h>


static void ex1(void);
static void ex2(void);


int main(void)
{
  ex1();
  ex2();
  return EXIT_SUCCESS;
}


#define NMAX 20
#define MMAX 20
#define TDX MMAX
#define TDQ MMAX+1

static void ex1(void)
{
  double rss, tol;
  Integer i, ip, rank, j, m, n;
  double df;
  Boolean  svd;
  char weight, meanc;
  Nag_IncludeMean  mean;
  double b[MMAX], cov[(MMAX*MMAX+MMAX)/2], h[NMAX], p[MMAX*(MMAX+2)],
  q[NMAX][MMAX+1], res[NMAX], se[MMAX], com_ar[MMAX*MMAX+5*(MMAX-1)],
  wt[NMAX], x[NMAX][MMAX], y[NMAX];
  double *wtptr;
  Integer  sx[MMAX];

  Vprintf("g02dac Example 1 Program Results\n");
  /* Skip heading in data file */
  Vscanf("%*[^\n]");
  Vscanf("%ld %ld %c %c", &n, &m, &weight, &meanc);
  if (meanc=='m')
    mean = Nag_MeanInclude;
  else
    mean = Nag_MeanZero;
  if (n<=NMAX && m<MMAX)
    {
      if (weight=='w')
        {
          wtptr = wt;
          for (i=0; i<n; i++)
            {
              for (j=0; j<m; j++)
                Vscanf("%lf", &x[i][j]);
              Vscanf("%lf%lf", &y[i], &wt[i]);
            }
        }
      else
        {
          wtptr = (double *)0;
          for (i=0; i<n; i++)
            {
              for (j=0; j<m; j++)
                Vscanf("%lf", &x[i][j]);
              Vscanf("%lf", &y[i]);
            }
        }
      for (j=0; j<m; j++)
        Vscanf("%ld", &sx[j]);
      /* Calculate ip */
      ip = 0;
      if (mean==Nag_MeanInclude)
        ip += 1;
      for (i=0; i<m; i++)
        if (sx[i]>0) ip += 1;
      /* Set tolerance */
      tol = 0.00001e0;
      g02dac(mean, n, &x[0][0],  (Integer)TDX,  m,  sx,  ip, y,
             wtptr,  &rss,  &df,  b,  se,  cov,  res,  h,  &q[0][0],
             (Integer)(TDQ),  &svd,  &rank,  p,  tol,  com_ar,  NAGERR_DEFAULT);

      if (svd)
        Vprintf("Model not of full rank, rank = %4ld\n\n", rank);
      Vprintf("Residual sum of squares = %12.4e\n", rss);
      Vprintf("Degrees of freedom = %3.1f\n\n", df);
      Vprintf("Variable    Parameter estimate   Standard error\n\n");
      for (j=0; j<ip; j++)
        Vprintf("%6ld%20.4e%20.4e\n", j+1, b[j], se[j]);
      Vprintf("\n");
      Vprintf("    Obs         Residuals              h\n\n");
      for (i=0; i<n; i++)
        Vprintf("%6ld%20.4e%20.4e\n", i+1, res[i], h[i]);
    }
  else
    {
      Vfprintf(stderr, "One or both of m and n are out of range:"
               " m = %-3ld while  n = %-3ld\n", m, n);
      exit(EXIT_FAILURE);
    }
  return;
}



static void ex2(void)
{
  double rss, tol;
  Integer i, ip, rank, j, m, n, degree, digits;
  double df;
  Boolean  svd;
  Nag_IncludeMean  mean;
  double b[MMAX], cov[(MMAX*MMAX+MMAX)/2], h[NMAX], p[MMAX*(MMAX+2)],
  q[NMAX][MMAX+1], res[NMAX], se[MMAX], com_ar[MMAX*MMAX+5*(MMAX-1)],
  x[NMAX][MMAX], y[NMAX];
  double *wtptr = (double *)0;   /* don't use weights */
  Integer  sx[MMAX];

  Vprintf("\n\n\ng02dac Example 2 Program Results\n");
  /* Skip heading in data file */
  Vscanf(" %*[^\n]");

  /* Use mean = Nag_MeanInclude */

  mean = Nag_MeanInclude;

  Vscanf("%ld%ld%ld",&degree,&n,&digits);

  if (n<=NMAX)
    {

      /* Set tolerance */
      tol = pow(10.0, -(double)digits);
      m = degree;
      ip = degree + 1;

      for  (i = 0; i <ip-1; ++i)
        sx[i] = 1;

      for (i=0; i<n; i++)
        {
          Vscanf("%lf%lf", &x[i][degree-1],&y[i]);
          for (j=0; j <degree; ++j)
            x[i][j] = pow(x[i][degree-1],(double)(degree-j));
        }

      g02dac(mean, n, &x[0][0],  (Integer)TDX,  m,  sx,  ip, y,
             wtptr,  &rss,  &df,  b,  se,  cov,  res,  h,  &q[0][0],
             (Integer)(TDQ),  &svd,  &rank,  p,  tol,  com_ar,  NAGERR_DEFAULT);


      Vprintf("Regression estimates  (mean = Nag_MeanInclude) \n\n");
      Vprintf("Coefficient    Estimate         Standard error\n\n");
      for (j=1; j<ip; j++)
        Vprintf("a(%ld)%20.4e%20.4e\n", degree+1-j, b[j], se[j]);
      Vprintf("a(0)%20.4e%20.4e\n", b[0], se[0]);
      Vprintf("\n\n");

      /* Use mean = Nag_MeanZero */

      mean = Nag_MeanZero;

      m = degree + 1;
      for  (i = 0; i <ip; ++i)
        sx[i] = 1;

      for (i=0; i<n; i++)
        x[i][m-1] = 1.0;

      g02dac(mean, n, &x[0][0],  (Integer)TDX,  m,  sx,  ip, y,
             wtptr,  &rss,  &df,  b,  se,  cov,  res,  h,  &q[0][0],
             (Integer)(TDQ),  &svd,  &rank,  p,  tol,  com_ar,  NAGERR_DEFAULT);

      Vprintf("Regression estimates  (mean = Nag_MeanZero) \n\n");
      Vprintf("Coefficient    Estimate         Standard error\n\n");
      for (j=0; j<ip; j++)
        Vprintf("a(%ld)%20.4e%20.4e\n", degree-j, b[j], se[j]);
      Vprintf("\n\n");
    }
  else
    {
      Vfprintf(stderr, "n is out of range, n = %-3ld\n", n);
      exit(EXIT_FAILURE);
    }
return;
}