/* nag_glm_predict (g02gpc) Example Program.
 *
 * Copyright 2008, Numerical Algorithms Group.
 *
 * Mark 9, 2009.
 */
/* Pre-processor includes */
#include <stdio.h>
#include <math.h>
#include <ctype.h>
#include <nag.h>
#include <nagx04.h>
#include <nag_stdlib.h>
#include <nagg02.h>

#define T_X(I, J) t_x[(I) *t_tdx + J]
#define X(I, J)   x[(I) *tdx + J]
int main(int argc, char *argv[])
{
  FILE              *fpin, *fpout;
  /* Integer scalar and array declarations */
  Integer           i, ip, j, m, n, t_n, tdx, t_tdx, print_iter;
  Integer           exit_status = 0, tdv, rank, lx, lt_x, lv;
  Integer           *sx = 0;
  /* NAG structures */
  Nag_Link          link;
  Nag_IncludeMean   mean;
  Nag_Boolean       vfobs, weight, t_weight, ioffset, t_ioffset;
  Nag_Distributions errfn;
  NagError          fail;
  /* Character scalar and array declarations */
  char              sioffset[10], st_ioffset[10], sweight[10], st_weight[10];
  char              slink[13], smean[16], svfobs[10];
  /* Double scalar and array declarations */
  double            rss, scale, ex_power, df;
  double            *b = 0, *cov = 0, *eta = 0, *offset = 0, *t_offset = 0;
  double            *pred = 0, *se = 0, *seeta = 0, *sepred = 0, *binom_t = 0;
  double            *v = 0, *wt = 0, *x = 0, *y = 0, *t_x = 0, *t_wt = 0;
  /* Set control parameters */
  double            eps = 0.000001;
  double            tol = 0.00005;
  Integer           max_iter = 10;

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

  /* Check for command-line IO options */
  fpin = nag_example_file_io(argc, argv, "-data", NULL);
  fpout = nag_example_file_io(argc, argv, "-results", NULL);

  fprintf(fpout, "nag_glm_predict (g02gpc) Example Program Results\n");

  /* Skip headings in data file */
  fscanf(fpin, "%*[^\n] ");
  fscanf(fpin, "%*[^\n] ");
  /* Read in training data for model that will be used for prediction */
  fscanf(fpin, "%s %s %s %s %ld %ld %lf %ld%*[^\n] ",
         slink, smean, st_ioffset, st_weight, &t_n, &m, &scale, &print_iter);
  /*
   * nag_enum_name_to_value (x04nac).
   * Converts NAG enum member name to value
   */
  link = (Nag_Link) nag_enum_name_to_value(slink);
  mean = (Nag_IncludeMean) nag_enum_name_to_value(smean);
  t_ioffset = (Nag_Boolean) nag_enum_name_to_value(st_ioffset);
  t_weight = (Nag_Boolean) nag_enum_name_to_value(st_weight);

  t_tdx = m;
  lt_x = t_tdx * t_n;

  /* Allocate memory */
  if (t_weight)
    {
      if (!(t_wt = NAG_ALLOC(t_n, double)))
        {
          fprintf(fpout, "Allocation failure\n");
          exit_status = -1;
          goto END;
        }
    }
  if (t_ioffset)
    {
      if (!(t_offset = NAG_ALLOC(t_n, double)))
        {
          fprintf(fpout, "Allocation failure\n");
          exit_status = -1;
          goto END;
        }
    }
  if (!(t_x = NAG_ALLOC(lt_x, double)) ||
      !(y = NAG_ALLOC(t_n, double)) ||
      !(sx = NAG_ALLOC(m, Integer)))
    {
      fprintf(fpout, "Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Read in the data */
  for (i = 0; i < t_n; i++)
    {
      for (j = 0; j < m; j++)
        fscanf(fpin, "%lf", &T_X(i, j));
      fscanf(fpin, "%lf", &y[i]);
      if (t_ioffset)
        fscanf(fpin, "%lf", &t_offset[i]);
      if (t_weight)
        fscanf(fpin, "%lf", &t_wt[i]);
      fscanf(fpin, "%*[^\n] ");
    }

  for (j = 0; j < m; j++)
    fscanf(fpin, "%ld%*[^\n] ", &sx[j]);

  if (link == Nag_Expo)
    fscanf(fpin, "%lf%*[^\n] ", &ex_power);
  else
    ex_power = 0.0;

  /* Calculate ip */
  ip = 0;
  for (j = 0; j < m; j++)
    if (sx[j] > 0) ip++;
  if (mean == Nag_MeanInclude)
    ip++;

  tdv = ip+6;
  lv = tdv * t_n;

  if (!(b = NAG_ALLOC(ip, double)) ||
      !(v = NAG_ALLOC(lv, double)) ||
      !(se = NAG_ALLOC(ip, double)) ||
      !(cov = NAG_ALLOC(ip*(ip+1)/2, double)))
    {
      fprintf(fpout, "Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Call nag_glm_normal (g02gac) to fit model to training data */
  nag_glm_normal(link, mean, t_n, t_x, t_tdx, m, sx, ip, y, t_wt, t_offset,
                 &scale, ex_power, &rss, &df, b, &rank, se, cov, v, tdv,
                 tol, max_iter, print_iter, "", eps, &fail);
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_glm_normal (g02gac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  /* Display parameter estimates for training data */
  fprintf(fpout,
          "\nResidual sum of squares = %12.4g, Degrees of freedom = %2f\n\n",
          rss, df);
  fprintf(fpout, "       Estimate     Standard error\n\n");
  for (i = 0; i < ip; i++)
    fprintf(fpout, " %14.4f %14.4f\n", b[i], se[i]);
  fprintf(fpout, "\n");

  /* Skip second lot of headings in data file */
  fscanf(fpin, "%*[^\n] ");

  /* Read in data to predict from and check array sizes */
  fscanf(fpin, "%ld %s %s %s%*[^\n] ", &n, svfobs, sioffset, sweight);
  /*
   * nag_enum_name_to_value (x04nac).
   * Converts NAG enum member name to value
   */
  vfobs = (Nag_Boolean) nag_enum_name_to_value(svfobs);
  ioffset = (Nag_Boolean) nag_enum_name_to_value(sioffset);
  weight = (Nag_Boolean) nag_enum_name_to_value(sweight);

  if (weight)
    {
      if (!(wt = NAG_ALLOC(n, double)))
        {
          fprintf(fpout, "Allocation failure\n");
          exit_status = -1;
          goto END;
        }
    }
  if (ioffset)
    {
      if (!(offset = NAG_ALLOC(n, double)))
        {
          fprintf(fpout, "Allocation failure\n");
          exit_status = -1;
          goto END;
        }
    }

  tdx = m;
  lx = tdx * n;

  if (!(x = NAG_ALLOC(lx, double)) ||
      !(eta = NAG_ALLOC(n, double)) ||
      !(seeta = NAG_ALLOC(n, double)) ||
      !(pred = NAG_ALLOC(n, double)) ||
      !(sepred = NAG_ALLOC(n, double)))
    {
      fprintf(fpout, "Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  for (i = 0; i < n; i++)
    {
      for (j = 0; j < m; j++)
        fscanf(fpin, "%lf", &X(i, j));
      if (offset)
        fscanf(fpin, "%lf", &offset[i]);
      if (weight)
        fscanf(fpin, "%lf", &wt[i]);
      fscanf(fpin, "%*[^\n] ");
    }

  /* Using nag_glm_normal (g02gac) to fit training model, so error structure
     is normal */
  errfn = Nag_Normal;

  /* Call nag_glm_predict (g02gpc) to calculate predictions */
  nag_glm_predict(errfn, link, mean, n, x, tdx, m, sx, ip, binom_t, offset,
                  wt, scale, ex_power, b, cov, vfobs, eta, seeta, pred,
                  sepred, &fail);
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_glm_predict (g02gpc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  /* Display predicted values */
  fprintf(fpout,
          "   I      ETA         SE(ETA)      Predicted  SE(Predicted)\n");
  fprintf(fpout, "\n");
  for (i = 0; i < n; i++)
    {
      fprintf(fpout, " %3ld) %10.5f    %10.5f    %10.5f    %10.5f\n", i+1,
              eta[i], seeta[i], pred[i], sepred[i]);
    }

 END:
  if (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  if (t_wt) NAG_FREE(t_wt);
  if (t_x) NAG_FREE(t_x);
  if (y) NAG_FREE(y);
  if (sx) NAG_FREE(sx);
  if (b) NAG_FREE(b);
  if (v) NAG_FREE(v);
  if (se) NAG_FREE(se);
  if (cov) NAG_FREE(cov);
  if (wt) NAG_FREE(wt);
  if (x) NAG_FREE(x);
  if (offset) NAG_FREE(offset);
  if (eta) NAG_FREE(eta);
  if (seeta) NAG_FREE(seeta);
  if (pred) NAG_FREE(pred);
  if (sepred) NAG_FREE(sepred);

  return exit_status;
}