/* nag_robust_m_estim_1var_usr (g07dcc) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 */

#include <math.h>
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg07.h>

#ifdef __cplusplus
extern "C"
{
#endif
  static double NAG_CALL chi(double t, Nag_Comm *comm);
  static double NAG_CALL psi(double t, Nag_Comm *comm);

#ifdef __cplusplus
}
#endif

int main(void)
{

  /* Scalars */
  double beta, sigma, sigsav, thesav, theta, tol;
  Integer exit_status, i, isigma, maxit, n, nit;
  NagError fail;
  Nag_Comm comm;

  /* Arrays */
  static double ruser[2] = { -1.0, -1.0 };
  double *rs = 0, *x = 0;

  INIT_FAIL(fail);

  exit_status = 0;
  printf("nag_robust_m_estim_1var_usr (g07dcc) Example Program Results\n");

  /* For communication with user-supplied functions: */
  comm.user = ruser;

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  scanf("%" NAG_IFMT "%*[^\n] ", &n);
  /* Allocate memory */
  if (!(rs = NAG_ALLOC(n, double)) || !(x = NAG_ALLOC(n, double)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  printf("\n");

  for (i = 1; i <= n; ++i) {
    scanf("%lf", &x[i - 1]);
  }
  scanf("%*[^\n] ");
  scanf("%lf%" NAG_IFMT "%*[^\n] ", &beta, &maxit);
  printf("          Input parameters     Output parameters\n");
  printf("isigma   sigma   theta   tol    sigma  theta\n");
  while (scanf("%" NAG_IFMT "%lf%lf%lf%*[^\n] ", &isigma, &sigma,
               &theta, &tol) != EOF) {
    sigsav = sigma;
    thesav = theta;

    /* nag_robust_m_estim_1var_usr (g07dcc).
     * Robust estimation, M-estimates for location and scale
     * parameters, user-defined weight functions
     */
    nag_robust_m_estim_1var_usr(chi, psi, isigma, n, x, beta, &theta,
                                &sigma, maxit, tol, rs, &nit, &comm, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_robust_m_estim_1var_usr (g07dcc).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }

    printf("%3" NAG_IFMT "%3s%8.4f%8.4f%7.4f", isigma, "", sigsav,
           thesav, tol);
    printf("%8.4f%8.4f\n", sigma, theta);
  }
END:
  NAG_FREE(rs);
  NAG_FREE(x);
  return exit_status;
}

static double NAG_CALL psi(double t, Nag_Comm *comm)
{
  /* Scalars */
  double abst;
  double ret_val;

  /* Hampel's Piecewise Linear Function. */
  if (comm->user[0] == -1.0) {
    printf("(User-supplied callback psi, first invocation.)\n");
    comm->user[0] = 0.0;
  }
  abst = fabs(t);
  if (abst < 4.5) {
    if (abst <= 3.0) {
      ret_val = MIN(1.5, abst);
    }
    else {
      ret_val = (4.5 - abst) * 1.5 / 1.5;
    }
    if (t < 0.0) {
      ret_val = -ret_val;
    }
  }
  else {
    ret_val = 0.0;
  }
  return ret_val;
} /* psi */

double NAG_CALL chi(double t, Nag_Comm *comm)
{
  /* Scalars */
  double abst, ps;
  double ret_val;

  /* Huber's chi function. */
  if (comm->user[1] == -1.0) {
    printf("(User-supplied callback chi, first invocation.)\n");
    comm->user[1] = 0.0;
  }
  abst = fabs(t);
  ps = MIN(1.5, abst);
  ret_val = ps * ps / 2;
  return ret_val;
}