/* nag_lars (g02mac) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 25, 2014.
 */
/* Pre-processor includes */
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg02.h>

int main(void)
{
  /* Integer scalar and array declarations */
  Integer i, j, k, ip, ldb, ldd, m, mnstep, n, nstep, lropt;
  Integer *isx = 0;
  Integer exit_status = 0;

  /* NAG structures and types */
  NagError fail;
  Nag_LARSModelType mtype;
  Nag_LARSPreProcess pred, prey;

  /* Double scalar and array declarations */
  double *b = 0, *d = 0, *fitsum = 0, *y = 0, *ropt = 0;

  /* Character scalar and array declarations */
  char cmtype[40], cpred[40], cprey[40];

  /* Initialise the error structure */
  INIT_FAIL(fail);

  printf("nag_lars (g02mac) Example Program Results\n\n");

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

  /* Read in the problem size */
  scanf("%ld%ld%*[^\n] ",&n, &m);

  /* Read in the model specification */
  scanf("%39s%39s%39s%ld%*[^\n] ", cmtype, cpred, cprey, &mnstep);
  mtype = (Nag_LARSModelType) nag_enum_name_to_value(cmtype);
  pred = (Nag_LARSPreProcess) nag_enum_name_to_value(cpred);
  prey = (Nag_LARSPreProcess) nag_enum_name_to_value(cprey);

  /* Using all variables */
  isx = 0;

  /* Optional arguments (using defaults) */
  lropt = 0;
  ropt = 0;

  /* Allocate memory for the data */
  ldd = n;
  if (!(y = NAG_ALLOC(n, double)) ||
      !(d = NAG_ALLOC(ldd*m, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Read in the data */
  for (i = 0; i < n; i++)
    {
      for (j = 0; j < m; j++)
        {
          scanf("%lf",&d[j*ldd + i]);
        }
      scanf("%lf",&y[i]);
    }
  scanf("%*[^\n] ");

  /* Allocate output arrays */
  ldb = m;
  if (!(b = NAG_ALLOC(ldb*(mnstep+2), double)) ||
      !(fitsum = NAG_ALLOC(6*(mnstep+1), double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Call nag_lars (g02mac) to fit the model */
  nag_lars(mtype, pred, prey, n, m, d, ldd, isx, y, mnstep, &ip, &nstep,
           b, ldb, fitsum, ropt, lropt, &fail);
  if (fail.code != NE_NOERROR)
    {
      if (fail.code != NW_OVERFLOW_WARN && fail.code != NW_POTENTIAL_PROBLEM
          && fail.code != NW_LIMIT_REACHED)
        {
          printf("Error from nag_lars (g02mac).\n%s\n", fail.message);
          exit_status = 11;
          goto END;
        }
      else
        {
          printf("Warning from nag_lars (g02mac).\n%s\n", fail.message);
          exit_status = 2;
        }
    }

  /* Display the parameter estimates */
  printf("  Step ");
  for (i = 0; i < MAX(ip-2,0)*5;i++) printf(" ");
  printf(" Parameter Estimate\n ");
  for (i = 0; i < 5+ip*10; i++) printf("-");
  printf("\n");
  for (k = 0; k < nstep; k++)
    {
      printf("  %3ld",k + 1);
      for (j = 0; j < ip; j++)
        {
          printf(" %9.3f",b[k*ldb + j]);
        }
      printf("\n");
    }
  printf("\n");
  printf(" alpha: %9.3f\n", fitsum[6*nstep]);
  printf("\n");
  printf("  Step     Sum      RSS       df       Cp       Ck     Step Size\n ");
  for (i = 0; i < 64; i++) printf("-");
  printf("\n");
  for (k = 0; k < nstep; k++)
    {
      printf("  %3ld %9.3f %9.3f %6.0f  %9.3f %9.3f %9.3f\n",
        k+1,fitsum[k*6],fitsum[k*6 + 1],fitsum[k*6 + 2],
        fitsum[k*6 + 3],fitsum[k*6 + 4],fitsum[k*6 + 5]);
    }
  printf("\n");
  printf(" sigma^2: %9.3f\n", fitsum[nstep*6+4]);

 END:
  NAG_FREE(y);
  NAG_FREE(d);
  NAG_FREE(b);
  NAG_FREE(fitsum);
  NAG_FREE(ropt);

  return(exit_status);
}