/* nag_surviv_cox_model (g12bac) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 *
 * NAG C Library
 *
 * Mark 26.1, 2017.
 */

#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg02.h>
#include <nagg12.h>

int main(void)
{
  Integer exit_status = 0, i, *ic = 0, ip, ip1, iprint, irank, *isi = 0, j, m,
         maxit;
  Integer n, nd, ndmax, ns, *sz = 0, tdsur, tdv;
  double dev, df, tol;
  double *b = 0, *cov = 0, *offset = 0, *omega = 0, *res = 0, *sc = 0;
  double *se = 0, *sur = 0, *t = 0, *tp = 0, *v = 0, *y = 0, *z = 0;
  NagError fail;

#define Z(I, J) z[((I) -1)*m +(J) -1]

  INIT_FAIL(fail);

  printf("nag_surviv_cox_model (g12bac) Example Program Results\n");

  /* Skip heading in data file */
  scanf("%*[^\n]");
  scanf("%" NAG_IFMT " %" NAG_IFMT " %" NAG_IFMT " %" NAG_IFMT " %" NAG_IFMT
        " ", &n, &m, &ns, &maxit, &iprint);
  ndmax = 42;
  tdsur = MAX(1, ns);
  if (!(z = NAG_ALLOC(n * m, double))
      || !(sz = NAG_ALLOC(m, Integer))
      || !(t = NAG_ALLOC(n, double))
      || !(ic = NAG_ALLOC(n, Integer))
      || !(omega = NAG_ALLOC(n, double))
      || !(isi = NAG_ALLOC(n, Integer))
      || !(res = NAG_ALLOC(n, double))
      || !(sur = NAG_ALLOC(ndmax * tdsur, double))
      || !(tp = NAG_ALLOC(ndmax, double)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  if (ns > 0) {
    for (i = 1; i <= n; ++i) {
      scanf("%lf", &t[i - 1]);
      for (j = 1; j <= m; ++j)
        scanf("%lf", &Z(i, j));
      scanf("%" NAG_IFMT "", &ic[i - 1]);
      scanf("%" NAG_IFMT "", &isi[i - 1]);
    }
  }
  else {
    for (i = 1; i <= n; ++i) {
      scanf("%lf", &t[i - 1]);
      for (j = 1; j <= m; ++j)
        scanf("%lf", &Z(i, j));
      scanf("%" NAG_IFMT "", &ic[i - 1]);
    }
  }
  for (i = 1; i <= m; ++i)
    scanf("%" NAG_IFMT "", &sz[i - 1]);
  scanf("%" NAG_IFMT "", &ip);
  ip1 = ip + 1;
  if (!(b = NAG_ALLOC(ip1, double))
      || !(se = NAG_ALLOC(ip1, double))
      || !(sc = NAG_ALLOC(ip1, double))
      || !(cov = NAG_ALLOC(ip1 * (ip1 + 1) / 2, double))
      || !(tdv = ip1 + 6)
      || !(v = NAG_ALLOC(n * tdv, double))
      || !(y = NAG_ALLOC(n, double))
      || !(offset = NAG_ALLOC(n, double)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  tol = 5e-5;
  for (i = 1; i <= n; ++i) {
    y[i - 1] = 1.0 - (double) ic[i - 1];
    offset[i - 1] = log(t[i - 1]);
  }
  /* nag_glm_poisson (g02gcc).
   * Fits a generalized linear model with Poisson errors
   */
  nag_glm_poisson(Nag_Log, Nag_MeanInclude, n, z, m, m, sz, ip1, y, 0, offset,
                  0.0, &dev, &df, b, &irank, se, cov, v, tdv, tol,
                  maxit, 0, 0, 0.0, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_glm_poisson (g02gcc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  for (i = 1; i <= ip; ++i)
    b[i - 1] = b[i];
  if (irank != ip + 1)
    printf(" WARNING: covariates not of full rank\n");

  /* nag_surviv_cox_model (g12bac).
   * Fits Cox's proportional hazard model
   */
  nag_surviv_cox_model(n, m, ns, z, m, sz, ip, t, ic, (double *) 0,
                       isi, &dev, b, se, sc, cov, res, &nd, tp, sur, tdsur,
                       ndmax, tol, maxit, iprint, "", &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_surviv_cox_model (g12bac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  printf("\n");
  printf(" Parameter      Estimate       Standard Error\n");
  printf("\n");
  for (i = 1; i <= ip; ++i)
    printf("%6" NAG_IFMT "          %8.4f          %8.4f\n",
           i, b[i - 1], se[i - 1]);
  printf("\n");
  printf(" Deviance = %13.4e\n", dev);
  printf("\n");
  printf("    Time     Survivor Function\n");
  printf("\n");
  ns = MAX(ns, 1);
  for (i = 1; i <= nd; ++i) {
    printf("%10.0f", tp[i - 1]);
    for (j = 1; j <= ns; ++j)
      printf("     %8.4f%s", sur[(i - 1) * tdsur + j - 1], j % 3 ? "" : "\n");
    printf("\n");
  }

END:
  NAG_FREE(z);
  NAG_FREE(sz);
  NAG_FREE(t);
  NAG_FREE(ic);
  NAG_FREE(omega);
  NAG_FREE(isi);
  NAG_FREE(res);
  NAG_FREE(sur);
  NAG_FREE(tp);
  NAG_FREE(b);
  NAG_FREE(se);
  NAG_FREE(sc);
  NAG_FREE(cov);
  NAG_FREE(v);
  NAG_FREE(y);
  NAG_FREE(offset);

  return exit_status;
}