/* nag_blgm_lm_design_matrix (g22ycc) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.2, 2016.
 */
/* Pre-processor includes */
#include <stdio.h>
#include <string.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg02.h>
#include <nagg22.h>
#include <nagx04.h>

#define MAX_FORMULA_LEN 200
#define MAX_VNAME_LEN 200
#define MAX_PLAB_LEN 200
#define MAX_OPTION_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);
void print_x(Nag_IncludeIntercept intcpt,char *const plab[],Integer nobs,
             Integer mx,const double x[],Integer ldx,const char *text);

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, mncat;
  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;

  /* Character scalar and array declarations */
  char formula[MAX_FORMULA_LEN], tcontrast[MAX_OPTION_LEN],
    optstr[MAX_OPTION_LEN];
  char *pch;
  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_design_matrix (g22ycc) 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;
  }

  /* Read in the contrast to use for all parameters */
  read_line(tcontrast,MAX_OPTION_LEN);

  /* Set up the option string */
  mncat = MAX_OPTION_LEN;
  strcpy(optstr, "Contrast = ");
  mncat -= strlen(optstr) + 1;
  strncat(optstr, tcontrast, mncat);

  /* Call nag_blgm_optset (g22zmc) to set the contrast optional argument */
  nag_blgm_optset(hform,optstr,&fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blgm_optset (g22zmc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* 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 */
  lddat = nobs;
  sddat = nvar;
  if (!(dat = NAG_ALLOC(lddat*sddat, 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("%*[^\n] ");

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

  /* 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;
  }
  /* ... End of constructing the design matrix */

  /* Start of getting the column labels for X ... */
  /* Get size of output arrays used by nag_blgm_lm_submodel (g22ydc) */
  lvinfo = 3;
  lisx = lplab = lenlab = 0;
  nag_blgm_lm_submodel(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 only want plab) */
  lisx = lvinfo = 0;
  lplab = tvinfo[1];
  if (!(plab = NAG_ALLOC(ip, char *))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  lenlab = MAX_PLAB_LEN;
  for (i = 0; i < ip; i++)
    if (!(plab[i] = NAG_ALLOC(lenlab, char))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Call nag_blgm_lm_submodel (g22ydc) to generate the labels */
  nag_blgm_lm_submodel(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 column labels for X */

  /* Display the design matrix */
  print_x(intcpt,plab,nobs,mx,x,ldx,"First");
  printf("\n");

  /* Read in the name of the variable whose contrasts need to be changed, */
  /* the value to change them to, remove comments and set the contrast */
  /* optional argument for the specified variable */
  for (;read_line(tcontrast,MAX_OPTION_LEN);) {
    /* Add an equals sign between variable name and contrast name */
    pch = strstr(tcontrast," ");
    *pch = '=';

    /* Set up the option string */
    mncat = MAX_OPTION_LEN;
    strncpy(optstr, "Contrast:",mncat);
    mncat -= strlen(optstr) + 1;
    strncat(optstr, tcontrast, mncat);

    /* Call nag_blgm_optset (g22zmc) to set the contrast optional argument */
    nag_blgm_optset(hform,optstr,&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 regenerate */
  /* 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;
  }

  /* Call nag_blgm_lm_submodel (g22ydc) to regenerate the column labels */
  nag_blgm_lm_submodel(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;
  }

  /* Display the design matrix */
  print_x(intcpt,plab,nobs,mx,x,ldx,"Second");

 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(levels);
  for (i = 0; i < lvnames; i++)
    NAG_FREE(vnames[i]);
  NAG_FREE(vnames);
  for (i = 0; i < ip; 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;
  }
}

void print_x(Nag_IncludeIntercept intcpt,char *const plab[],Integer nobs,
             Integer mx,const double x[],Integer ldx,const char *text) {
  /* Print the transpose of the first 10 rows of the design matrix */

  /* Integer scalar and array declarations */
  Integer max_rows = 10, i, j, pnobs, si;

  /* plab holds the labels for the model parameters, so includes a label */
  /* for the mean effect, if one is present. As the mean effect is not */
  /* being explicitly included in the design matrix, we may need to skip */
  /* the first element of plab (which will always be the label for the */
  /* mean effect if one is present) */
  si = (intcpt == Nag_Intercept) ? 1 : 0;

  /* Printing the first max_rows rows of the design matrix */
  pnobs = MIN(max_rows,nobs);

  /* Display the design matrix */
  printf(" Transpose of First %ld Rows of the %s Design Matrix (X)\n",
         pnobs, text);
  printf(" Column Name   ");
  for (i = 0; i < pnobs; i++)
    printf("   %2"NAG_IFMT, i+1);
  printf("\n ");
  for (i = 0; i < 15+pnobs*5; i++)
    printf("-");
  printf("\n");
  for (j = 0; j < mx; j++) {
    printf(" %-15s", plab[j+si]);
    for (i = 0; i < pnobs; i++)
      printf(" %4.1f", X(i,j));
    printf("\n");
  }

  /* Display the intercept flag using nag_enum_value_to_name (x04nbc) */
  printf(" Intercept flag = %s\n", nag_enum_value_to_name(intcpt));
}