/* nag_robust_m_regsn_estim(g02hac) 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 regtypec, Nag_RegType *regtype, char psifunc,
                     Nag_PsiFun *psifun, char sigestc, Nag_SigmaEst *sigma_est);

#define NMAX 8
#define MMAX 3
#define TDC MMAX
#define TDX MMAX

int main(void)
{
  char  regtypec, psifunc, sigestc, covmat_estc;
  Nag_RegType regtype;
  Nag_PsiFun psifun;
  Nag_SigmaEst sigma_est;
  Nag_CovMatrixEst covmat_est;
  Integer i, j, m, n;
  double  sigma;
  double  c[MMAX][MMAX], theta[MMAX], x[NMAX][MMAX], y[NMAX];
  double cpsi;
  double hpsi[3];
  double cucv;
  double dchi;
  double rs[NMAX];
  double wt[NMAX];
  double tol;
  Integer max_iter;
  Integer print_iter;
  double info[4];
  static NagError fail;

  Vprintf("g02hac Example Program Results\n\n");
  /* Skip heading in data file */
  Vscanf("%*[^\n]");
  Vscanf("%ld %ld", &n, &m);
  if ((n> 0 && n<=NMAX) && (m>0 && m<=MMAX))
    {
      /* Read in x and y */
      for (i=0; i<n; i++)
        {
          for (j=0; j<m; j++)
            Vscanf("%lf", &x[i][j]);
          Vscanf("%lf", &y[i]);
        }
      /* Read in control parameters */
      Vscanf(" %c %c %c", &regtypec, &psifunc, &sigestc);

      /* Check control parameters */
      set_enum(regtypec, &regtype, psifunc, &psifun, sigestc, &sigma_est);

      /* Read in appropriate weight function parameters. */
      if (regtype !=  Nag_HuberReg)
        Vscanf(" %c %lf", &covmat_estc, &cucv);

      if (covmat_estc == 'A')
        covmat_est = Nag_CovMatAve;
      else
        covmat_est = Nag_CovMatObs;

      if (psifun != Nag_Lsq)
        {
          if (psifun == Nag_HuberFun)
            Vscanf("%lf", &cpsi);
          else
            cpsi = 0.0;
          if (psifun == Nag_HampelFun)
            for (j=0; j<3; j++)
              Vscanf("%lf", &hpsi[j]);
          if (sigma_est == Nag_SigmaChi)
            Vscanf("%lf", &dchi);
        }
      /* Set values of remaining parameters */
      tol = 5e-5;
      max_iter = 50;
      /* Change print_iter to a positive value if monitoring information
       * is required
       */
      print_iter = 1;
      sigma = 1.0e0;
      for (i = 0; i < m; ++i)
        theta[i] = 0.0e0;

      g02hac(regtype, psifun, sigma_est, covmat_est, n, m, &x[0][0],
             (Integer)TDX, y, cpsi, hpsi, cucv, dchi, theta, &sigma,
             &c[0][0], (Integer)TDC, rs, wt, tol, max_iter, print_iter,
             "", info, &fail);

      if ((fail.code == NE_NOERROR) || (fail.code == NE_THETA_ITER_EXCEEDED) ||
          (fail.code == NE_LSQ_FAIL_CONV) || (fail.code == NE_MAT_SINGULAR) ||
          (fail.code == NE_WT_LSQ_NOT_FULL_RANK) ||
          (fail.code == NE_REG_MAT_SINGULAR) ||
          (fail.code == NE_COV_MAT_FACTOR_ZERO) ||
          (fail.code == NE_VAR_THETA_LEQ_ZERO) ||
          (fail.code == NE_VAR_THETA_LEQ_ZERO) ||
          (fail.code == NE_ERR_DOF_LEQ_ZERO) ||
          (fail.code == NE_ESTIM_SIGMA_ZERO))
        {
          if (fail.code != NE_NOERROR)
            {
              Vprintf("%s\n", fail.message);
              Vprintf("    Some of the following results may be unreliable\n");
            }
          Vprintf("Sigma = %10.4f\n\n", sigma);
          Vprintf("       Theta      Standard errors\n\n");
          for (j = 0; j < m; ++j)
            Vprintf("%12.4f %13.4f\n", theta[j], c[j][j]);
          Vprintf("\n     Weights      Residuals\n\n");
          for (i = 0; i < n; ++i)
            Vprintf("%12.4f %13.4f\n", wt[i], rs[i]);
        }
      else
        {
          Vprintf("%s\n",fail.message);
        }
    }
  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 regtypec, Nag_RegType *regtype, char psifunc,
                     Nag_PsiFun *psifun, char sigestc, Nag_SigmaEst *sigma_est)
{
  if (toupper(regtypec) == 'H' || toupper(regtypec) == 'M' ||
      toupper(regtypec) == 'S')
    {
      switch (toupper(regtypec))
        {
        case ('H'):
          *regtype = Nag_HuberReg;
          break;
        case ('M'):
          *regtype = Nag_MallowsReg;
          break;
        case ('S'):
          *regtype = Nag_SchweppeReg;
          break;
        default:
          ;
        }
    }
  else
    {
      Vfprintf(stderr, "The parameter regtype has an invalid value:"
               " regtype = %c\n", regtypec);

      exit(EXIT_FAILURE);
    }

  if (toupper(psifunc) == 'L' || toupper(psifunc) == 'H' ||
      toupper(psifunc) == 'P' || toupper(psifunc) == 'A'
      || toupper(psifunc) == 'T')
    {
      switch (toupper(psifunc))
        {
        case ('L'):
          *psifun = Nag_Lsq;
          break;
        case ('H'):
          *psifun = Nag_HuberFun;
          break;
        case ('P'):
          *psifun = Nag_HampelFun;
          break;
        case ('A'):
          *psifun = Nag_AndrewFun;
          break;
        case ('T'):
          *psifun = Nag_TukeyFun;
          break;
        default:
          ;
        }
    }
  else
    {
      Vfprintf(stderr, "The parameter psifun has an invalid value: "
               "psifun = %c\n", psifunc);
      exit(EXIT_FAILURE);
    }
  if (toupper(sigestc) == 'R' || toupper(sigestc) == 'C' ||
      toupper(sigestc) == 'X')
    {
      switch (toupper(sigestc))
        {
        case ('R'):
          *sigma_est = Nag_SigmaRes;
          break;
        case ('C'):
          *sigma_est = Nag_SigmaConst;
          break;
        case ('X'):
          *sigma_est = Nag_SigmaChi;
          break;
        default:
          ;
        }
    }
  else
    {
      Vfprintf(stderr, "The parameter sigma_est has an invalid value:"
               " sigma_est = %c\n", sigestc);

      exit(EXIT_FAILURE);
    }
  return;
}