/* nag_glm_gamma(g02gdc) 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 10
#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, ip;
  double  ex_power, scale;
  Integer sx[MMAX];
  double  b[MMAX], v[NMAX][TDV], wt[NMAX], x[NMAX][MMAX], y[NMAX];
  double  *wtptr, *offsetptr=(double *)0;
  Integer max_iter;
  Integer print_iter;
  double eps;
  double tol;
  double df;
  double dev;
  Integer rank;
  double se[MMAX], cov[MMAX*(MMAX+1)/2];
  static  NagError fail;

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

  /* 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", &y[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", &y[i]);
            }
        }
      for (j=0; j<m; j++)
        Vscanf("%ld", &sx[j]);

      /* Calculate ip */
      ip = 0;
      for (j=0; j<m; j++)
        if (sx[j]>0) ip += 1;
      if (mean == Nag_MeanInclude)
        ip += 1;
      if (link == Nag_Expo)
        Vscanf("%lf", &ex_power);
      else
        ex_power = 0.0;

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

      g02gdc(link, mean, n, &x[0][0], (Integer)TDX, m, sx, ip, y,
             wtptr, offsetptr, &scale, ex_power, &dev, &df, b, &rank,
             se, cov, &v[0][0], (Integer)TDV, tol, max_iter,
             print_iter, "", eps, &fail);

      if (fail.code == NE_NOERROR || 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<ip; i++)
            Vprintf("%14.4f%14.4f\n", b[i], se[i]);
          Vprintf("\n");
          Vprintf("     y   fitted value  Residual  Leverage\n\n");
          for (i = 0; i < n; ++i)
            {
              Vprintf("%7.1f%10.2f%12.4f%10.3f\n", 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) == 'E' || toupper(linkc) == 'I' || toupper(linkc) == 'L'
      || toupper(linkc) == 'S' || toupper(linkc) == 'R')
    {
      switch (toupper(linkc))
        {
        case ('E'):
          *link = Nag_Expo;
          break;
        case ('I'):
          *link = Nag_Iden;
          break;
        case ('L'):
          *link = Nag_Log;
          break;
        case ('S'):
          *link = Nag_Sqrt;
          break;
        case ('R'):
          *link = Nag_Reci;
          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;
}