/* nag_surviv_cox_model (g12bac) Example Program.
 *
 * Copyright 2000 Numerical Algorithms Group.
 *
 *
 * NAG C Library
 *
 * Mark 6, 2000.
 * Mark 7, revised, 2001.
 *
 */

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

/* Table of constant values */

static double c_b42 = 0.0;
static Integer c__0 = 0;

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

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

  INIT_FAIL(fail);
  Vprintf("g12bac Example Program Results\n");

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

  Vscanf("%ld %ld %ld %ld %ld ", &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)))
    {
      Vprintf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  if (ns > 0)
    {
      for (i = 1; i <= n; ++i)
        {
          Vscanf("%lf", &t[i - 1]);
          for (j = 1; j <= m; ++j)
            Vscanf("%lf", &Z(i,j));
          Vscanf("%ld", &ic[i - 1]);
          Vscanf("%ld", &isi[i - 1]);
        }
    }
  else
    {
      for (i = 1; i <= n; ++i)
        {
          Vscanf("%lf", &t[i - 1]);
          for (j = 1; j <= m; ++j)
            Vscanf("%lf", &Z(i,j));
          Vscanf("%ld", &ic[i - 1]);
        }
    }
  for (i = 1; i <= m; ++i)
    Vscanf("%ld", &sz[i - 1]);
  Vscanf("%ld", &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)))
    {
      Vprintf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  tol = 5e-5;
  i__1 = n;
  for (i = 1; i <= n; ++i)
    {
      y[i - 1] = 1.0 - (double) ic[i - 1];
      offset[i-1] = log (t[i - 1]);
    }
  g02gcc(Nag_Log, Nag_MeanInclude, n, z, m, m, sz, ip1, y, 0, offset,
         c_b42, &dev, &df, b, &irank, se, cov, v, tdv, tol,
         maxit, c__0, 0, c_b42, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from 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)
    Vprintf(" WARNING: covariates not of full rank\n");

  g12bac(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)
    {
      Vprintf("Error from g12bac.\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  Vprintf("\n");
  Vprintf(" Parameter      Estimate       Standard Error\n");
  Vprintf("\n");
  for (i = 1; i <= ip; ++i)
    Vprintf("%6ld          %8.4f          %8.4f\n",
            i, b[i - 1], se[i - 1]);
  Vprintf("\n");
  Vprintf(" Deviance = %13.4e\n", dev);
  Vprintf("\n");
  Vprintf("    Time     Survivor Function\n");
  Vprintf("\n");
  ns = MAX(ns, 1);
  for (i = 1; i <= nd; ++i)
    {
      Vprintf("%10.0f", tp[i - 1]);
      for (j = 1; j <= ns; ++j)
        Vprintf("     %8.4f%s", sur[(i-1)*tdsur + j-1], j%3?"":"\n");
      Vprintf("\n");
    }
 END:
  if (z) NAG_FREE(z);
  if (sz) NAG_FREE(sz);
  if (t) NAG_FREE(t);
  if (ic) NAG_FREE(ic);
  if (omega) NAG_FREE(omega);
  if (isi) NAG_FREE(isi);
  if (res) NAG_FREE(res);
  if (sur) NAG_FREE(sur);
  if (tp) NAG_FREE(tp);
  if (b) NAG_FREE(b);
  if (se) NAG_FREE(se);
  if (sc) NAG_FREE(sc);
  if (cov) NAG_FREE(cov);
  if (v) NAG_FREE(v);
  if (y) NAG_FREE(y);
  if (offset) NAG_FREE(offset);
  return exit_status;
}