/* nag_robust_m_regsn_estim (g02hac) Example Program.
 *
 * Copyright 1996 Numerical Algorithms Group.
 *
 * Mark 4, 1996.
 * Mark 8 revised, 2004.
 *
 */

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

static int set_enum(char regtypec, Nag_RegType *regtype, char psifunc,
                    Nag_PsiFun *psifun, char sigestc, Nag_SigmaEst *sigma_est);

#define C(I,J) c[(I)*tdc + J]
#define X(I,J) x[(I)*tdx + J]
int main(void)
{
  Integer exit_status=0, i, j, m, max_iter, n, print_iter, tdc, tdx;
  NagError fail;
  Nag_CovMatrixEst covmat_est;
  Nag_PsiFun psifun;
  Nag_RegType regtype;
  Nag_SigmaEst sigma_est;
  char covmat_estc, psifunc, regtypec, sigestc;
  double *c=0, cpsi, cucv, dchi, *hpsi=0, *info=0, *rs=0, sigma, *theta=0;
  double tol, *wt=0, *x=0, *y=0;

  INIT_FAIL(fail);
  Vprintf("nag_robust_m_regsn_estim (g02hac) Example Program Results\n\n");
  /* Skip heading in data file */
  Vscanf("%*[^\n]");
  Vscanf("%ld %ld", &n, &m);
  if (n>1 && (m>=1 && m<=n))
    {
      if ( !( c = NAG_ALLOC(m*m, double)) ||
           !( theta = NAG_ALLOC(m, double)) ||
           !( x = NAG_ALLOC(n*m, double)) ||
           !( y = NAG_ALLOC(n, double)) ||
           !( rs = NAG_ALLOC(n, double)) ||
           !( wt = NAG_ALLOC(n, double)) ||
           !( info = NAG_ALLOC(4, double)) ||
           !( hpsi = NAG_ALLOC(3, double)) )
        {
          Vprintf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
      tdc = m;
      tdx = m;
    }
  else
    {
      Vprintf("Invalid n or m.\n");
      exit_status = 1;
      return exit_status;
    }
  /* 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 */
  if (set_enum(regtypec, &regtype, psifunc, &psifun, sigestc, &sigma_est))
    {
      Vprintf("Invalid enum type.\n");
      exit_status = 1;
      goto END;
    }

  /* 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;

  /* nag_robust_m_regsn_estim (g02hac).
   * Robust regression, standard M-estimates
   */
  nag_robust_m_regsn_estim(regtype, psifun, sigma_est, covmat_est, n, m, x,
                           tdx, y, cpsi, hpsi, cucv, dchi, theta, &sigma,
                           c, 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("Error from nag_robust_m_regsn_estim (g02hac).\n%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("Error from nag_robust_m_regsn_estim (g02hac).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }
 END:
  if (c) NAG_FREE(c);
  if (theta) NAG_FREE(theta);
  if (x) NAG_FREE(x);
  if (y) NAG_FREE(y);
  if (rs) NAG_FREE(rs);
  if (wt) NAG_FREE(wt);
  if (info) NAG_FREE(info);
  if (hpsi) NAG_FREE(hpsi);
  return exit_status;
}

static int set_enum(char regtypec, Nag_RegType *regtype, char psifunc,
                    Nag_PsiFun *psifun, char sigestc, Nag_SigmaEst *sigma_est)
{
  Integer exit_status=0;
  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
    {
      Vprintf("The parameter regtype has an invalid value:"
              " regtype = %c\n", regtypec);

      exit_status = 1;
    }

  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
    {
      Vprintf("The parameter psifun has an invalid value: "
              "psifun = %c\n", psifunc);
      exit_status = 1;
    }
  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
    {
      Vprintf("The parameter sigma_est has an invalid value:"
              " sigma_est = %c\n", sigestc);

      exit_status = 1;
    }
  return exit_status;
}