/* nag_blgm_lm_submodel (g22ydc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.3, 2023.
 */
/* Pre-processor includes */
#include <nag.h>
#include <stdio.h>
#include <string.h>

#define MAX_FORMULA_LEN 200
#define MAX_VNAME_LEN 200
#define MAX_PLAB_LEN 200
#define MAX_CVALUE_LEN 200

#define DAT(I, J) dat[(J)*lddat + (I)]

char *read_line(char formula[], Integer nchar);
Integer construct_labels(Integer ip, char **plab[], char *const vnames[],
                         Integer vinfo[]);
Integer fit_lm(void *hform, Nag_IncludeIntercept intcpt, Integer nobs,
               Integer mx, double x[], Integer ldx, Integer isx[], Integer ip,
               double y[], char *plab[]);

int main(void) {
  /* Integer scalar and array declarations */
  Integer i, j, ip = 0, lddat, ldx, lisx, lplab = 0, lvinfo, lvnames = 0, mx,
                nobs, nvar, sddat, sdx, lenlab;
  Integer exit_status = 0;
  Integer *isx = 0, *levels = 0, *vinfo = 0;
  Integer tvinfo[3];

  /* Nag Types */
  NagError fail;
  Nag_IncludeIntercept intcpt;

  /* Double scalar and array declarations */
  double *dat = 0, *x = 0, *y = 0;

  /* Character scalar and array declarations */
  char formula[MAX_FORMULA_LEN];
  char **vnames = 0, **plab = 0;

  /* Void pointers */
  void *hform = 0, *hddesc = 0, *hxdesc = 0;

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

  printf("nag_blgm_lm_submodel (g22ydc) Example Program Results\n\n");

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

  /* Read in size of the data matrix and number of variable labels supplied */
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &nobs, &nvar,
        &lvnames);

  /* Allocate memory */
  if (!(levels = NAG_ALLOC(nvar, Integer)) ||
      !(vnames = NAG_ALLOC(lvnames, char *))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  for (i = 0; i < lvnames; i++)
    if (!(vnames[i] = NAG_ALLOC(MAX_VNAME_LEN, char))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Read in number of levels and names for the variables */
  for (i = 0; i < nvar; i++) {
    scanf("%" NAG_IFMT "", &levels[i]);
  }
  scanf("%*[^\n] ");

  if (lvnames > 0) {
    for (i = 0; i < lvnames; i++)
      scanf("%50s", vnames[i]);
    scanf("%*[^\n] ");
  }

  /* Call nag_blgm_lm_describe_data (g22ybc) to get a description of */
  /* the data matrix */
  nag_blgm_lm_describe_data(&hddesc, nobs, nvar, levels, lvnames, vnames,
                            &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blgm_lm_describe_data (g22ybc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Read in the data matrix and response variable */
  lddat = nobs;
  sddat = nvar;
  if (!(dat = NAG_ALLOC(lddat * sddat, double)) ||
      !(y = NAG_ALLOC(nobs, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  for (i = 0; i < nobs; i++) {
    for (j = 0; j < nvar; j++)
      scanf("%lf", &DAT(i, j));
    scanf("%lf", &y[i]);
  }
  scanf("%*[^\n] ");

  /* Read in the formula for the full model, remove comments and */
  /* call nag_blgm_lm_formula (g22yac) to parse it */
  read_line(formula, MAX_FORMULA_LEN);
  nag_blgm_lm_formula(&hform, formula, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blgm_lm_formula (g22yac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Start of constructing the design matrix ... */

  /* Call nag_blgm_optset (g22zmc) to alter the storage order of X as */
  /* nag_correg_linregm_fit uses VAROBS storage  */
  nag_blgm_optset(hform, "Storage Order = VAROBS", &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blgm_optset (g22zmc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Call nag_blgm_lm_design_matrix (g22ycc) to get the size of */
  /* the design matrix */
  ldx = 0;
  sdx = 0;
  nag_blgm_lm_design_matrix(hform, hddesc, dat, lddat, sddat, &hxdesc, x, ldx,
                            sdx, &mx, &fail);
  if (fail.code != NW_ARRAY_SIZE && fail.code != NW_ALTERNATIVE) {
    printf("Error from nag_blgm_lm_design_matrix (g22ycc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Allocate design matrix */
  ldx = mx;
  sdx = nobs;
  if (!(x = NAG_ALLOC(ldx * sdx, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Call nag_blgm_lm_design_matrix (g22ycc) to generate the design matrix */
  nag_blgm_lm_design_matrix(hform, hddesc, dat, lddat, sddat, &hxdesc, x, ldx,
                            sdx, &mx, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blgm_lm_design_matrix (g22ycc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }
  /* ... End of constructing the design matrix */

  /* Start of getting the isx vector and information on parameter labels ... */
  /* Get size of output arrays used by nag_blgm_lm_submodel (g22ydc) */
  lvinfo = 3;
  lisx = lplab = lenlab = 0;
  nag_blgm_lm_submodel(Nag_X, hform, hxdesc, &intcpt, &ip, lisx, isx, lplab,
                       plab, lenlab, lvinfo, tvinfo, &fail);
  if (fail.code != NW_ARRAY_SIZE) {
    printf("Error from nag_blgm_lm_submodel (g22ydc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Allocate output arrays (we already know that lisx = mx, but */
  /* nag_blgm_lm_submodel returns it just in case) */
  lisx = tvinfo[0];
  lplab = tvinfo[1];
  /* We don't need vinfo as we are using labels in plab */
  lvinfo = 0;
  if (!(isx = NAG_ALLOC(lisx, Integer)) || !(plab = NAG_ALLOC(lplab, char *))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  lenlab = MAX_PLAB_LEN;
  for (i = 0; i < lplab; i++)
    if (!(plab[i] = NAG_ALLOC(lenlab, char))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Call nag_blgm_lm_submodel (g22ydc) to get the isx flag */
  /* and parameter labels */
  nag_blgm_lm_submodel(Nag_X, hform, hxdesc, &intcpt, &ip, lisx, isx, lplab,
                       plab, lenlab, lvinfo, vinfo, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blgm_lm_submodel (g22ydc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  /* ... End of getting the isx vector and information on parameter labels */

  /* Fit a regression model and print the results */
  if ((exit_status = fit_lm(hform, intcpt, nobs, mx, x, ldx, isx, ip, y, plab)))
    goto END;

  /* Read in the formula for the sub-model, remove comments and */
  /* call nag_blgm_lm_formula (g22yac) to parse it */
  read_line(formula, MAX_FORMULA_LEN);
  nag_blgm_lm_formula(&hform, formula, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blgm_lm_formula (g22yac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Call nag_blgm_lm_submodel (g22ydc) to get the isx flag and parameter */
  /* labels as the new model has to be a sub-model of the original one, the */
  /* output arrays, isx, plab and vinfo can be reused as they will be of */
  /* sufficient size */
  nag_blgm_lm_submodel(Nag_SubModelX, hform, hxdesc, &intcpt, &ip, lisx, isx,
                       lplab, plab, lenlab, lvinfo, vinfo, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blgm_lm_submodel (g22ydc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  printf("\n\n");

  /* Fit the sub-model and print the results */
  exit_status = fit_lm(hform, intcpt, nobs, mx, x, ldx, isx, ip, y, plab);

END:
  /* Call nag_blgm_handle_free (g22zac) to clean-up the g22 handles */
  nag_blgm_handle_free(&hform, &fail);
  nag_blgm_handle_free(&hddesc, &fail);
  nag_blgm_handle_free(&hxdesc, &fail);

  NAG_FREE(dat);
  NAG_FREE(x);
  NAG_FREE(y);
  NAG_FREE(levels);
  for (i = 0; i < lvnames; i++)
    NAG_FREE(vnames[i]);
  NAG_FREE(vnames);
  for (i = 0; i < lplab; i++)
    NAG_FREE(plab[i]);
  NAG_FREE(plab);
  NAG_FREE(isx);
  NAG_FREE(vinfo);
  return (exit_status);
}

char *read_line(char formula[], Integer nchar) {
  /* Read in a line from stdin and remove any comments */
  char *pch;

  /* Read in the model formula */
  if (fgets(formula, nchar, stdin)) {
    /* Strip comments from formula */
    pch = strstr(formula, "::");
    if (pch)
      *pch = '\0';
    return formula;
  } else {
    return 0;
  }
}

Integer fit_lm(void *hform, Nag_IncludeIntercept intcpt, Integer nobs,
               Integer mx, double x[], Integer ldx, Integer isx[], Integer ip,
               double y[], char *plab[]) {
  /* Perform a multiple linear regression using */
  /* nag_correg_linregm_fit (g02dac) */

  /* Integer scalar and array declarations */
  Integer i, rank, ldq, exit_status = 0, lcvalue, ivalue;

  /* NAG types */
  Nag_Boolean svd;
  NagError fail;
  Nag_IncludeMean mean;
  Nag_VariableType optype;

  /* Double scalar and array declarations */
  double rss, tol, df, rvalue;
  double *b = 0, *cov = 0, *h = 0, *p = 0, *q = 0, *res = 0, *se = 0, *wt = 0,
         *com_ar = 0;

  /* Character scalar and array declarations */
  char cvalue[MAX_CVALUE_LEN];

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

  /* We are assuming un-weighted data */
  ldq = (ip + 1);
  if (!(b = NAG_ALLOC(ip, double)) || !(se = NAG_ALLOC(ip, double)) ||
      !(cov = NAG_ALLOC((ip * ip + ip) / 2, double)) ||
      !(res = NAG_ALLOC(nobs, double)) || !(h = NAG_ALLOC(nobs, double)) ||
      !(q = NAG_ALLOC(nobs * ldq, double)) ||
      !(p = NAG_ALLOC(ip * (ip + 2), double)) ||
      !(com_ar = NAG_ALLOC(ip * ip + 5 * (ip - 1), double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Use suggested value for tolerance */
  tol = 0.000001;

  mean = (intcpt == Nag_Intercept) ? Nag_MeanInclude : Nag_MeanZero;

  /* Call nag_correg_linregm_fit (g02dac) to fit a regression model */
  nag_correg_linregm_fit(mean, nobs, x, ldx, mx, isx, ip, y, wt, &rss, &df, b,
                         se, cov, res, h, q, ldq, &svd, &rank, p, tol, com_ar,
                         &fail);

  /* Call nag_blgm_optget (g22znc) to get the formula for the model */
  /* being fit */
  lcvalue = MAX_CVALUE_LEN;
  nag_blgm_optget(hform, "Formula", &ivalue, &rvalue, cvalue, lcvalue, &optype,
                  &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blgm_optget (g22znc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Display the results */
  printf(" Model: %s\n", cvalue);
  printf("                                Parameter   Standard\n");
  printf(" Coefficients                   Estimate      Error\n");
  printf(" ---------------------------------------------------\n");
  for (i = 0; i < ip; i++)
    printf(" %-30s %7.3f     %7.3f\n", plab[i], b[i], se[i]);
  printf(" ---------------------------------------------------\n");
  printf(" Residual sum of squares = %9.4f\n", rss);
  printf(" Degrees of freedom      = %9.0f\n", df);

END:
  NAG_FREE(h);
  NAG_FREE(res);
  NAG_FREE(wt);
  NAG_FREE(b);
  NAG_FREE(cov);
  NAG_FREE(p);
  NAG_FREE(q);
  NAG_FREE(com_ar);
  NAG_FREE(se);
  return exit_status;
}