/* nag_correg_lars_param (g02mcc) Example Program.
 *
 * Copyright 2019 Numerical Algorithms Group.
 *
 * Mark 27.0, 2019.
 */
/* Pre-processor includes */
#include <stdio.h>
#include <nag.h>

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

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

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

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

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

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

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

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

  /* Read in the model specification */
  scanf("%39s%39s%39s%" NAG_IFMT "%*[^\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_correg_lars (g02mac) to fit the model */
  nag_correg_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_correg_lars (g02mac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    else {
      printf("Warning from nag_correg_lars (g02mac).\n%s\n", fail.message);
      exit_status = 2;
    }
  }

  /* Read in the number of additional parameter estimates required and the
     way they will be specified */
  scanf("%39s%" NAG_IFMT "%*[^\n] ", cktype, &lnk);
  ktype = (Nag_LARSTargetType) nag_enum_name_to_value(cktype);

  ldnb = ip;
  if (!(nk = NAG_ALLOC(lnk, double)) || !(nb = NAG_ALLOC(ip * lnk, double)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read in the target values */
  for (i = 0; i < lnk; i++) {
    scanf("%lf", &nk[i]);
  }
  scanf("%*[^\n] ");

  /* nag_correg_lars_param (g02mcc) Calculate the additional parameter estimates */
  nag_correg_lars_param(nstep, ip, b, ldb, fitsum, ktype, nk, lnk, nb, ldnb, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_correg_lars_param (g02mcc).\n%s\n", fail.message);
    exit_status = 3;
    goto END;
  }

  printf("Parameter Estimates from nag_correg_lars (g02mac)\n");
  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("  %3" NAG_IFMT "", k + 1);
    for (j = 0; j < ip; j++) {
      printf(" %9.3f", b[k * ldb + j]);
    }
    printf("\n");
  }
  printf("\n");

  printf("Additional Parameter Estimates from nag_correg_lars_param (g02mcc)\n");
  printf("   nk  ");
  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 < lnk; k++) {
    printf("  %3.1f", nk[k]);
    for (j = 0; j < ip; j++) {
      printf(" %9.3f", nb[k * ldnb + j]);
    }
    printf("\n");
  }

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

  return (exit_status);
}