/* nag_regsn_quant_linear (g02qgc) Example Program.
 *
 * Copyright 2011, Numerical Algorithms Group.
 *
 * Mark 23, 2011.
 */
/* 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 CH(i, j, k) ch[k*ip*ip + j*ip + i]

#define LOPTSTR 80
#define ENUMLEN 30

int main(void)
{
  /* Integer scalar and array declarations */
  Integer              lseed = 1, liopts = 100, lopts = 100, lcvalue = LOPTSTR;
  Integer              exit_status = 0;
  Integer              genid, 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 structures */
  NagError             fail;
  Nag_OrderType        order = Nag_RowMajor;
  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;
  char                 optstr[LOPTSTR], corder[ENUMLEN], cintcpt[ENUMLEN],
                       cweighted[ENUMLEN], cgenid[ENUMLEN];
  char                 *clabs=0, **clabsc=0;

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

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

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

  /* Read in the problem size */
  scanf("%s%*[^\n]", corder);
  scanf("%s%s%*[^\n]", cintcpt, cweighted);
  scanf("%ld%ld%ld%*[^\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 initialise if it is */
      if (strcmp(semeth, "BOOTSTRAP XY") == 0)
        {
          /* Read in the generator ID and a seed */
          scanf("%s%ld%ld%*[^\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;
            }

          /* Initialise 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("%ld ", 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(" %3ld%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(" %3ld%10.3f\n", j+1, b[l*ip+j]);
        }
      printf("\n\n");
      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");
      /* 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:
  
  if (info) NAG_FREE(info);
  if (iopts) NAG_FREE(iopts);
  if (isx) NAG_FREE(isx);
  if (state) NAG_FREE(state);
  if (b) NAG_FREE(b);
  if (bl) NAG_FREE(bl);
  if (bu) NAG_FREE(bu);
  if (ch) NAG_FREE(ch);
  if (dat) NAG_FREE(dat);
  if (opts) NAG_FREE(opts);
  if (res) NAG_FREE(res);
  if (tau) NAG_FREE(tau);
  if (wt) NAG_FREE(wt);
  if (y) NAG_FREE(y);
  if (cvalue) NAG_FREE(cvalue);
  if (clabs) NAG_FREE(clabs);
  if (clabsc) NAG_FREE(clabsc);
  
  return(exit_status);
}