/* nag_blgm_lm_formula (g22yac) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.1, 2024.
 */
/* Pre-processor includes */
#include <nag.h>
#include <stdio.h>
#include <string.h>

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

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

char *read_line(char formula[], Integer nchar);

int main(void) {
  /* Integer scalar and array declarations */
  Integer i, j, ivalue, lddat, ldx, lvnames = 0, mx, nobs, nvar, sddat, sdx,
                                    lcvalue;
  Integer exit_status = 0;
  Integer *levels = 0;

  /* Nag Types */
  NagError fail;
  Nag_VariableType optype;

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

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

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

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

  printf("nag_blgm_lm_formula (g22yac) Example Program Results\n\n");

  /* Skip heading in data file */
  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;
  }

  /* Call nag_blgm_optget (g22znc) to extract the parsed formula */
  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;
  }

  printf(" Formula: %s\n", cvalue);
  printf("\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] ");

  /* Call nag_blgm_optset (g22zmc) to set optional arguments */
  /* Want the design matrix to include an explicit term for the mean effect */
  nag_blgm_optset(hform, "Explicit Mean = Yes", &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 = nobs;
  sdx = mx;
  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;
  }

  /* Display the design matrix */
  printf(" Design Matrix (X)\n");
  for (i = 0; i < nobs; i++) {
    for (j = 0; j < mx; j++)
      printf(" %4.1f", X(i, j));
    printf("\n");
  }

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);
  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;
  }
}