/* nag_regsn_quant_linear (g02qgc) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 */
/* Pre-processor includes */
#include <stdio.h>
#include <string.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg02.h>
#include <nagg05.h>
#include <nagx04.h>

#define DAT(i,j)    dat[(order==Nag_RowMajor) ? (i*pddat+j) : (j*pddat+i)]

#define LOPTSTR 80

int main(void)
{
  /* Integer scalar and array declarations */
  Integer lseed = 1, liopts = 100, lopts = 100, lcvalue = LOPTSTR;
  Integer exit_status = 0;
  Integer i, ip, ivalue, j, l, lc, lstate, loptstr,
         m, n, ntau, subid, tdch, pddat;
  Integer *info = 0, *iopts = 0, *isx = 0, *state = 0;
  Integer seed[1];
  Nag_BaseRNG genid;

  /* NAG structures */
  NagError fail;
  Nag_OrderType order;
  Nag_IncludeIntercept intcpt;
  Nag_Boolean weighted;
  Nag_VariableType optype;

  /* Double scalar and array declarations */
  double df, rvalue;
  double *b = 0, *bl = 0, *bu = 0, *ch = 0, *dat = 0,
         *opts = 0, *res = 0, *tau = 0, *wt = 0, *y = 0;

  /* Character scalar and array declarations */
  char semeth[30], *poptstr, *cvalue = 0;
  char optstr[LOPTSTR], corder[40], cintcpt[40], cweighted[40], cgenid[40];
  char *clabs = 0, **clabsc = 0;

  /* Initialize the error structure to print out any error messages */
  INIT_FAIL(fail);

  printf("nag_regsn_quant_linear (g02qgc) Example Program Results\n\n");
  fflush(stdout);

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

  /* Read in the problem size */
  scanf("%39s%*[^\n]", corder);
  scanf("%39s%39s%*[^\n]", cintcpt, cweighted);
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n]", &n, &m, &ntau);
  order = (Nag_OrderType) nag_enum_name_to_value(corder);
  intcpt = (Nag_IncludeIntercept) nag_enum_name_to_value(cintcpt);
  /* weighted is a Nag_Boolean flag used in this example program to indicate
   * whether weights are being supplied (weighted=Nag_TRUE)
   * or not (weighted=Nag_FALSE)
   */
  weighted = (Nag_Boolean) nag_enum_name_to_value(cweighted);

  pddat = (order == Nag_RowMajor) ? m : n;

  /* Allocate memory for input arrays */
  if (!(y = NAG_ALLOC(n, double)) ||
      !(tau = NAG_ALLOC(ntau, double)) ||
      !(isx = NAG_ALLOC(m, Integer)) ||
      !(dat = NAG_ALLOC(m * n, double)) ||
      !(cvalue = NAG_ALLOC(lcvalue, char)) ||
      !(clabs = NAG_ALLOC(10 * 10, char)) || !(clabsc = NAG_ALLOC(10, char *))
         )
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  if (weighted) {
    /* Data includes a weight */
    if (!(wt = NAG_ALLOC(n, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
    for (i = 0; i < n; i++) {
      for (j = 0; j < m; j++)
        scanf("%lf", &DAT(i, j));
      scanf("%lf%lf", &y[i], &wt[i]);
    }
    scanf("%*[^\n] ");
  }
  else {
    /* No weights supplied */
    for (i = 0; i < n; i++) {
      for (j = 0; j < m; j++)
        scanf("%lf", &DAT(i, j));
      scanf("%lf", &y[i]);
    }
    scanf("%*[^\n] ");
  }

  /* Read in variable inclusion flags and calculate IP */
  ip = (intcpt == Nag_Intercept) ? 1 : 0;
  for (j = 0; j < m; j++) {
    scanf("%" NAG_IFMT, &isx[j]);
    if (isx[j] == 1)
      ip++;
  }
  scanf("%*[^\n] ");

  /* Read in the quantiles required */
  for (l = 0; l < ntau; l++)
    scanf("%lf", &tau[l]);
  scanf("%*[^\n] ");

  /* Allocate memory for option arrays */
  if (!(opts = NAG_ALLOC(lopts, double)) ||
      !(iopts = NAG_ALLOC(liopts, Integer)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Initialize the optional argument array with nag_g02_opt_set (g02zkc) */
  nag_g02_opt_set("INITIALIZE = G02QG", iopts, liopts, opts, lopts, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_g02_opt_set (g02zkc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Read in any optional arguments. Reads in to the end of
     the input data, or until a blank line is reached */
  for (;;) {
    if (!fgets(optstr, LOPTSTR, stdin))
      break;

    /* Left justify the option */
    poptstr = (optstr + strspn(optstr, " \n\t"));
    /* Get the string length */
    loptstr = strlen(poptstr);
    if (poptstr[loptstr - 1] == '\n') {
      /* Remove any trailing line breaks */
      poptstr[(--loptstr)] = '\0';
    }
    else {
      /* Clear the rest of the line */
      scanf("%*[^\n] ");
    }

    /* Break if read in a blank line */
    if (!*(poptstr))
      break;

    /* Set the supplied option (g02zkc) */
    nag_g02_opt_set(optstr, iopts, liopts, opts, lopts, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_g02_opt_set (g02zkc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  }

  /* Allocate memory for the output arrays */
  if (!(b = NAG_ALLOC(ip * ntau, double)) ||
      !(info = NAG_ALLOC(ntau, Integer)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Query optional arguments via nag_g02_opt_get (g02zlc) and calculate which
   * of the optional arrays are required and their sizes
   * ...
   */
  nag_g02_opt_get("INTERVAL METHOD", &ivalue, &rvalue, cvalue, lcvalue,
                  &optype, iopts, opts, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_g02_opt_get (g02zlc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  strcpy(semeth, cvalue);
  if (strcmp(semeth, "NONE") != 0) {
    /* Require the intervals to be output */
    if (!(bl = NAG_ALLOC(ip * ntau, double)) ||
        !(bu = NAG_ALLOC(ip * ntau, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

    /* Decide whether the state array is required, and initialize if it is */
    if (strcmp(semeth, "BOOTSTRAP XY") == 0) {
      /* Read in the generator ID and a seed */
      scanf("%39s%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", cgenid, &subid,
            &seed[0]);
      genid = (Nag_BaseRNG) nag_enum_name_to_value(cgenid);

      /* Query the length of the state array (g05kfc) */
      lstate = 0;
      nag_rand_init_repeatable(genid, subid, seed, lseed, state, &lstate,
                               &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n",
               fail.message);
        exit_status = 1;
        goto END;
      }

      /* Allocate memory to state */
      if (!(state = NAG_ALLOC(lstate, Integer)))
      {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }

      /* Initialize the RNG (g05kfc) */
      nag_rand_init_repeatable(genid, subid, seed, lseed, state, &lstate,
                               &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n",
               fail.message);
        exit_status = 1;
        goto END;
      }
    }

    /* Calculate the size of the covariance matrix, ch. */
    tdch = 0;
    nag_g02_opt_get("MATRIX RETURNED", &ivalue, &rvalue, cvalue, lcvalue,
                    &optype, iopts, opts, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_g02_opt_get (g02zlc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

    if (strcmp(cvalue, "COVARIANCE") == 0) {
      tdch = ntau;
    }
    else if (strcmp(cvalue, "H INVERSE") == 0) {
      /* NB: If we are using bootstrap or IID errors then any request for
         H INVERSE is ignored */
      if (strcmp(semeth, "BOOTSTRAP XY") != 0 && strcmp(semeth, "IID") != 0)
        tdch = ntau + 1;
    }
    if (tdch > 0) {
      /* Need to allocate ch */
      if (!(ch = NAG_ALLOC(ip * ip * tdch, double)))
      {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }
    }

    /* Calculate the size of the residual array, res */
    nag_g02_opt_get("RETURN RESIDUALS", &ivalue, &rvalue, cvalue, lcvalue,
                    &optype, iopts, opts, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_g02_opt_get (g02zlc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

    if (strcmp(cvalue, "YES") == 0) {
      /* Need to allocate res */
      if (!(res = NAG_ALLOC(n * ntau, double)))
      {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }
    }
  }
  /* ...
   * end of handling the optional arguments, and allocating optional arrays
   */

  /* Call the model fitting routine (nag_regsn_quant_linear (g02qgc)) */
  nag_regsn_quant_linear(order, intcpt, n, m, dat, pddat, isx, ip, y, wt,
                         ntau, tau, &df, b, bl, bu, ch, res, iopts, opts,
                         state, info, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_regsn_quant_linear (g02qgc).\n%s\n", fail.message);
    if (fail.code == NW_POTENTIAL_PROBLEM) {
      printf("Additional error information: ");
      for (i = 0; i < ntau; i++)
        printf("%" NAG_IFMT " ", info[i]);
      printf("\n");
    }
    else {
      printf("Error from nag_regsn_quant_linear (g02qgc).\n%s\n",
             fail.message);
      exit_status = -1;
      goto END;
    }
  }

  /* Display the parameter estimates */
  for (l = 0; l < ntau; l++) {
    printf(" Quantile: %6.3f\n\n", tau[l]);
    if (bl && bu) {
      printf("         Lower   Parameter   Upper\n");
      printf("         Limit   Estimate    Limit\n");
      for (j = 0; j < ip; j++)
        printf(" %3" NAG_IFMT "%10.3f%10.3f%10.3f\n", j + 1, bl[l * ip + j],
               b[l * ip + j], bu[l * ip + j]);
    }
    else {
      printf("       Parameter\n");
      printf("       Estimate\n");
      for (j = 0; j < ip; j++)
        printf(" %3" NAG_IFMT "%10.3f\n", j + 1, b[l * ip + j]);
    }
    printf("\n\n");
    fflush(stdout);
    if (ch) {
      lc = l * ip * ip;
      if (tdch == ntau) {
        /* nag_gen_real_mat_print_comp (x04cbc).
         * Print real general matrix (comprehensive).
         */
        nag_gen_real_mat_print_comp(Nag_ColMajor, Nag_UpperMatrix,
                                    Nag_NonUnitDiag, ip, ip, &ch[lc], ip,
                                    "%9.2e", "Covariance matrix",
                                    Nag_NoLabels, 0, Nag_NoLabels, 0, 80,
                                    0, 0, &fail);
      }
      else {
        if (l == 0) {
          nag_gen_real_mat_print_comp(Nag_ColMajor, Nag_UpperMatrix,
                                      Nag_NonUnitDiag, ip, ip, ch, ip,
                                      "%9.2e", "J", Nag_NoLabels, 0,
                                      Nag_NoLabels, 0, 80, 0, 0, &fail);
          printf("\n");
        }
        lc = lc + ip * ip;
        nag_gen_real_mat_print_comp(Nag_ColMajor, Nag_UpperMatrix,
                                    Nag_NonUnitDiag, ip, ip, &ch[lc], ip,
                                    "%9.2e", "H inverse",
                                    Nag_NoLabels, 0, Nag_NoLabels, 0, 80,
                                    0, 0, &fail);
      }
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_gen_real_mat_print_comp (x04cbc).\n%s\n",
               fail.message);
        exit_status = 1;
        goto END;
      }
      printf("\n");
    }
  }

  if (res) {
    printf(" First 10 Residuals\n");
    fflush(stdout);
    /* set up column labels for matrix printer */
    for (l = 0; l < ntau; l++)
      sprintf(&clabs[10 * l], "%6.3f", tau[l]);
    for (l = 0; l < ntau; l++)
      clabsc[l] = &clabs[l * 10];
    /* nag_gen_real_mat_print_comp (x04cbc).
     * Print real general matrix (comprehensive).
     */
    nag_gen_real_mat_print_comp(Nag_ColMajor, Nag_GeneralMatrix,
                                Nag_NonUnitDiag, MIN(10, n), ntau, res, n,
                                "%10.5f", "                        Quantile",
                                Nag_IntegerLabels, NULL, Nag_CharacterLabels,
                                (const char **) clabsc, 80, 2, NULL, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_gen_real_mat_print_comp (x04cbc).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }
  }
  else {
    printf(" Residuals not returned\n");
  }

END:

  NAG_FREE(info);
  NAG_FREE(iopts);
  NAG_FREE(isx);
  NAG_FREE(state);
  NAG_FREE(b);
  NAG_FREE(bl);
  NAG_FREE(bu);
  NAG_FREE(ch);
  NAG_FREE(dat);
  NAG_FREE(opts);
  NAG_FREE(res);
  NAG_FREE(tau);
  NAG_FREE(wt);
  NAG_FREE(y);
  NAG_FREE(cvalue);
  NAG_FREE(clabs);
  NAG_FREE(clabsc);

  return (exit_status);
}