/* nag_glm_binomial(g02gbc) Example Program.
 *
 * Copyright 1996 Numerical Algorithms Group.
 *
 * Mark 4, 1996.
 *
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <ctype.h>
#include <nagg02.h>

static void set_enum(char linkc, Nag_Link *link, char meanc,
                     Nag_IncludeMean *mean);

#define NMAX 3
#define MMAX 2
#define TDX MMAX
#define TDV MMAX+6

int main(void)
{
  char  linkc, meanc, weightc;
  Nag_Link  link;
  Nag_IncludeMean  mean;
  Integer i, j, m, n, nvar;
  Integer ivar[MMAX];
  double  beta[MMAX], binom[NMAX], v[NMAX][TDV], wt[NMAX],
  x[NMAX][MMAX], y[NMAX];
  double  *wtptr, *offsetptr=(double *)0;
  Integer max_iter, print_iter;
  double tol, eps;
  Integer rank;
  double df, dev, se[MMAX], cov[MMAX*(MMAX+1)/2];
  static  NagError fail;

  Vprintf("g02gbc Example Program Results\n");
  /* Skip heading in data file */
  Vscanf("%*[^\n]");
  Vscanf(" %c %c %c %ld %ld %ld", &linkc, &meanc, &weightc,
         &n, &m, &print_iter);

  /* Check and set control parameters */
  set_enum(linkc, &link, meanc, &mean);
  if (n<=NMAX && m<MMAX)
    {
      if (toupper(weightc)=='W')
        {
          wtptr = wt;
          for (i=0; i<n; i++)
            {
              for (j=0; j<m; j++)
                Vscanf("%lf", &x[i][j]);
              Vscanf("%lf%lf%lf", &y[i], &binom[i], &wt[i]);
            }
        }
      else
        {
          wtptr = (double *)0;
          for (i=0; i<n; i++)
            {
              for (j=0; j<m; j++)
                Vscanf("%lf", &x[i][j]);
              Vscanf("%lf%lf", &y[i], &binom[i]);
            }
        }
      for (j=0; j<m; j++)
        Vscanf("%ld", &ivar[j]);

      /* Calculate nvar */
      nvar = 0;
      for (i=0; i<m; i++)
        if (ivar[i]>0) nvar += 1;
      if (mean == Nag_MeanInclude)
        nvar += 1;

      /* Set other control parameters */
      max_iter = 10;
      tol = 5e-5;
      eps = 1e-6;

      g02gbc(link, mean, n, &x[0][0], (Integer)TDX, m,
             ivar, nvar, y, binom, wtptr, offsetptr, &dev, &df, beta, &rank,
             se, cov,&v[0][0], (Integer)TDV, tol, max_iter, print_iter, "",
             eps, &fail);

      if (fail.code == NE_NOERROR || fail.code == NE_SVD_NOT_CONV ||
          fail.code == NE_LSQ_ITER_NOT_CONV ||
          fail.code == NE_RANK_CHANGED || fail.code == NE_ZERO_DOF_ERROR)
        {
          Vprintf("\nDeviance = %12.4e\n", dev);
          Vprintf("Degrees of freedom = %3.1f\n\n", df);
          Vprintf("       Estimate     Standard error\n\n");
          for (i=0; i<nvar; i++)
            Vprintf("%14.4f%14.4f\n", beta[i], se[i]);
          Vprintf("\n");
          Vprintf("     binom     y   fitted value  Residual  Leverage\n\n");
          for (i = 0; i < n; ++i)
            {
              Vprintf("%10.1f%7.1f%10.2f%12.4f%10.3f\n", binom[i], y[i],
                      v[i][1], v[i][4], v[i][5]);
            }
        }
      else
        {
          Vprintf("%s\n",fail.message);
          return EXIT_FAILURE;
        }
    }
  else
    {
      Vfprintf(stderr, "One or both of m and n are out of range:"
               " m = %-3ld while  n = %-3ld\n", m, n);

      return EXIT_FAILURE;
    }
  return EXIT_SUCCESS;
}

static void set_enum(char linkc, Nag_Link *link, char meanc,
                     Nag_IncludeMean *mean)
{
  if (toupper(linkc) == 'G' || toupper(linkc) == 'P' || toupper(linkc) == 'C')
    {
      switch (toupper(linkc))
        {
        case ('G'):
          *link = Nag_Logistic;
          break;
        case ('P'):
          *link = Nag_Probit;
          break;
        case ('C'):
          *link = Nag_Compl;
          break;
        default:
          ;
        }
    }
  else
    {
      Vfprintf(stderr, "The parameter link has an invalid value: link = %c\n",
               linkc);
      exit(EXIT_FAILURE);
    }
  if (toupper(meanc)=='M')
    *mean = Nag_MeanInclude;
  else if (toupper(meanc)=='Z')
    *mean = Nag_MeanZero;
  else
    {
      Vfprintf(stderr, "The parameter mean has an invalid value: mean = %c\n",
               meanc);
      exit(EXIT_FAILURE);
    }
  return;
}