/* nag_step_regsn (g02eec) Example Program.
 *
 * Copyright 2002 Numerical Algorithms Group.
 *
 * Mark 7, 2002.
 */

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

int main(void)
{

  /* Scalars */
  Nag_Boolean addvar;
  double chrss, f, fin, rss;
  Integer exit_status, i, idf, ifr, istep, j, m,
    maxip, n, nterm, pdq, pdx;
  NagError fail;
  Nag_OrderType order;
  Nag_IncludeMean mean_enum;
  char mean, weight;

  /* Arrays */
  char *newvar;
  double *exss=0, *p=0, *q=0, *wt=0, *x=0, *y=0, *wtptr=0;
  Integer *sx=0;

  char **free_vars, **model;

  char *vname[] = { "DAY", "BOD", "TKN", "TS", "TVS", "COD"};

#ifdef NAG_COLUMN_MAJOR
#define X(I,J) x[(J-1)*pdx + I - 1]
  order = Nag_ColMajor;
#else
#define X(I,J) x[(I-1)*pdx + J - 1]
  order = Nag_RowMajor;
#endif

  INIT_FAIL(fail);
  exit_status = 0;
  Vprintf("nag_step_regsn (g02eec) Example Program Results\n");

  /* Skip heading in data file */
  Vscanf("%*[^\n] ");

  Vscanf("%ld%ld ' %c ' ' %c '%*[^\n] ", &n, &m, &mean,
         &weight);
  maxip = m;

  /* Allocate memory */
  if ( !(exss = NAG_ALLOC(maxip, double)) ||
       !(p = NAG_ALLOC(maxip+1, double)) ||
       !(q = NAG_ALLOC(n * (maxip+2), double)) ||
       !(wt = NAG_ALLOC(n, double)) ||
       !(x = NAG_ALLOC(n * m, double)) ||
       !(y = NAG_ALLOC(n, double)) ||
       !(sx = NAG_ALLOC(m, Integer)) ||
       !(free_vars =  NAG_ALLOC(maxip, char *)) ||
       !(model =  NAG_ALLOC(maxip, char *))
       )
    {
      Vprintf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

#ifdef NAG_COLUMN_MAJOR
  pdx = n;
  pdq = n;
#else
  pdx = m;
  pdq = maxip+2;
#endif

  if (weight == 'W' || weight == 'w')
    {
      for (i = 1; i <= n; ++i)
        {
          for (j = 1; j <= m; ++j)
            Vscanf("%lf", &X(i,j));
          Vscanf("%lf%lf%*[^\n] ", &y[i - 1], &wt[i - 1]);
          wtptr = wt;
        }
    }
  else
    {
      for (i = 1; i <= n; ++i)
        {
          for (j = 1; j <= m; ++j)
            Vscanf("%lf", &X(i,j));
          Vscanf("%lf%*[^\n] ", &y[i - 1]);
        }
    }

  for (j = 1; j <= m; ++j)
    Vscanf("%ld", &sx[j - 1]);
  Vscanf("%*[^\n] ");

  Vscanf("%lf%*[^\n] ", &fin);

  if (mean == 'M')
    mean_enum = Nag_MeanInclude;
  else if (mean == 'Z')
    mean_enum = Nag_MeanZero;
  else
    {
      Vprintf("Incorrect value for mean: '%c'\n", mean);
      exit_status = -1;
      goto END;
    }

  Vprintf("\n");

  istep = 0;
  for (i = 1; i <= m; ++i)
    {
      /* nag_step_regsn (g02eec).
       * Fits a linear regression model by forward selection
       */
      nag_step_regsn(order, &istep, mean_enum, n, m, x, pdx, vname, sx, maxip,
                     y, wtptr, fin, &addvar, &newvar, &chrss, &f, model, &nterm,
                     &rss, &idf, &ifr, free_vars, exss, q, pdq, p, &fail);

      if (fail.code != NE_NOERROR)
        {
          Vprintf("Error from nag_step_regsn (g02eec).\n%s\n", fail.message);
          exit_status = 1;
          goto END;
        }

      Vprintf("Step %ld\n", istep);
      if (!addvar)
        {
          Vprintf("No further variables added maximum F =%7.2f\n", f);
          Vprintf("Free variables:        ");

          for (j = 1; j <= ifr; ++j)
            Vprintf("%3.3s %s", free_vars[j-1], j%6 == 0 || j == ifr ?"\n":" ");

          Vprintf("\n");
          Vprintf("Change in residual sums of squares for free variables:\n");

          Vprintf("                 ");
          for (j = 1; j <= ifr; ++j)
            {
              Vprintf("%9.4f", exss[j - 1]);
              Vprintf("%s", j%6 == 0 || j == ifr ?"\n":" ");
            }
          goto END;
        }
      else
        {
          Vprintf("Added variable is %3s\n", newvar);
          Vprintf("Change in residual sum of squares =%13.4e\n", chrss);
          Vprintf("F Statistic = %7.2f\n\n", f);
          Vprintf("Variables in model: ");

          for (j = 1; j <= nterm; ++j)
            Vprintf("%3s %s", model[j-1], j%6 == 0 || j == nterm ?"\n":" ");
          Vprintf("Residual sum of squares = %13.4e\n", rss);
          Vprintf("Degrees of freedom = %ld\n\n", idf);
          if (ifr == 0)
            {
              Vprintf("No free variables remaining\n");
              goto END;
            }

          Vprintf("%s%6s", "Free variables:  ", "");
          for (j = 1; j <= ifr; ++j)
            {
              Vprintf("%3.3s ", free_vars[j-1]);
              Vprintf(j%6 == 0 || j == ifr ?"\n":" ");
            }
          Vprintf("Change in residual sums of squares for free variables:\n");
          Vprintf("                 ");

          for (j = 1; j <= ifr; ++j)
            Vprintf("%9.4f%s", exss[j - 1], j%6 == 0 || j == ifr ?"\n":" ");
          Vprintf("\n");
        }
    }
 END:
  if (model) NAG_FREE(model);
  if (free_vars) NAG_FREE(free_vars);
  if (exss) NAG_FREE(exss);
  if (p) NAG_FREE(p);
  if (q) NAG_FREE(q);
  if (wt) NAG_FREE(wt);
  if (x) NAG_FREE(x);
  if (y) NAG_FREE(y);
  if (sx) NAG_FREE(sx);

  return exit_status;
}