/* nag_regsn_ridge_opt (g02kac) Example Program.
 *
 * Copyright 2008, Numerical Algorithms Group.
 *
 * Mark 9, 2009.
 */
/* Pre-processor includes */
#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg02.h>
#include <nagx04.h>

int main(void)
{
  /*Integer scalar and array declarations */
  Integer             exit_status = 0;
  Integer             df, i, ip, ip1, j, m, n, niter, one = 1;
  Integer             pdb, pdres, pdvif, pdx;
  Integer             *isx = 0;
  /*Double scalar and array declarations */
  double              h, nep, rss, tau, tol;
  double              *b = 0, *perr = 0, *res = 0, *vif = 0, *x = 0, *y = 0;
  /*Character scalar and array declarations */
  char                sopt[10], sorig[24], soptloo[12];
  /*NAG Types */
  Nag_OrderType       order;
  Nag_PredictError    opt;
  Nag_EstimatesOption orig;
  Nag_OptionLOO       optloo;
  NagError            fail;

  INIT_FAIL(fail);

  printf("%s\n",
          "nag_regsn_ridge_opt (g02kac) Example Program Results");
  /* Skip heading in data file*/
  scanf("%*[^\n] ");
  /* Read in data and check array limits*/
  scanf("%ld%ld%lf%s %lf%ld%s %s%*[^\n] ",
         &n, &m, &h, sopt, &tol, &niter, sorig, soptloo);
  opt = (Nag_PredictError) nag_enum_name_to_value(sopt);
  orig = (Nag_EstimatesOption) nag_enum_name_to_value(sorig);
  optloo = (Nag_OptionLOO) nag_enum_name_to_value(soptloo);

  #ifdef NAG_COLUMN_MAJOR
  pdx = n;
   #define X(I, J) x[(J-1)*pdx + I-1]
  order = Nag_ColMajor;
  #else
  pdx = m;
    #define X(I, J) x[(I-1)*pdx + J-1]
  order = Nag_RowMajor;
  #endif
  if (!(b = NAG_ALLOC(m+1, double)) ||
      !(perr = NAG_ALLOC(5, double)) ||
      !(res = NAG_ALLOC(n, double)) ||
      !(vif = NAG_ALLOC(m, double)) ||
      !(x = NAG_ALLOC(pdx*(order == Nag_RowMajor?n:m), double)) ||
      !(y = NAG_ALLOC(n, double)) ||
      !(isx = NAG_ALLOC(m, Integer)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  for (i = 1; i <= n; i++)
    {
      for (j = 1; j <= m; j++)
        scanf("%lf ", &X(i, j));
      scanf("%lf ", &y[i-1]);
    }
  scanf("%*[^\n] ");
  for (j = 0; j < m; j++)
    scanf("%ld ", &isx[j]);
  scanf("%*[^\n] ");

  /* Total number of variables.*/
  ip = 0;
  for (j = 0; j < m; j++)
    {
      if (isx[j] == 1)
        ip = ip+1;
    }
    #ifdef NAG_COLUMN_MAJOR
  pdb = n;
  pdres = n;
  pdvif = ip;
    #else
  pdb = one;
  pdres = one;
  pdvif = one;
    #endif
  /* Tolerance for setting singular values of H to zero.*/
  tau = 0.00e0;
  df = 0;
  /* Call function.*/
  /*
   * nag_regsn_ridge_opt (g02kac)
   * Ridge regression
   */
  nag_regsn_ridge_opt(order, n, m, x, pdx, isx, ip, tau, y, &h, opt, &niter,
                      tol, &nep, orig, b, vif, res, &rss, &df, optloo, perr,
                      &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_regsn_ridge_opt (g02kac).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }
  /* Print results:*/
  printf("\n");
  printf("%s   %10.4f\n", "Value of ridge parameter:", h);
  printf("\n");
  printf("%s   %13.4e\n", "Sum of squares of residuals:", rss);
  printf("%s  %5ld\n", "Degrees of freedom: ", df);
  printf("%s   %10.4f\n", "Number of effective parameters:", nep);
  printf("\n");
  ip1 = ip + 1;
  /*
   * nag_gen_real_mat_print (x04cac)
   * Print real general matrix (easy-to-use)
   */
  fflush(stdout);
  nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, ip1, one,
                         b, pdb, "Parameter estimates", 0, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_gen_real_mat_print (x04cac).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  printf("\n");
  printf("%s%ld\n", "Number of iterations:            ", niter);
  printf("\n");
  if (opt == Nag_GCV)
    {
      printf("%s\n", "Ridge parameter minimises GCV");
    }
  else if (opt == Nag_UEV)
    {
      printf("%s\n", "Ridge parameter minimises UEV");
    }
  else if (opt == Nag_FPE)
    {
      printf("%s\n", "Ridge parameter minimises FPE");
    }
  else if (opt == Nag_BIC)
    {
      printf("%s\n", "Ridge parameter minimises BIC");
    }
  printf("\n");
  printf("%s\n", "Estimated prediction errors:");
  printf("%s   %10.4f\n", "GCV    =", perr[0]);
  printf("%s   %10.4f\n", "UEV    =", perr[1]);
  printf("%s   %10.4f\n", "FPE    =", perr[2]);
  printf("%s   %10.4f\n", "BIC    =", perr[3]);
  if (optloo == Nag_WantLOO)
    {
      printf("%s   %10.4f\n", "LOO CV =", perr[4]);
    }
  printf("\n");

  /*
   * nag_gen_real_mat_print (x04cac)
   * Print real general matrix (easy-to-use)
   */
  fflush(stdout);
  nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, n, one,
                         res, pdres, "Residuals", 0, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_gen_real_mat_print (x04cac).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }
  printf("\n");
  /*
   * nag_gen_real_mat_print (x04cac)
   * Print real general matrix (easy-to-use)
   */
  fflush(stdout);
  nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, ip, one,
                         vif, pdvif, "Variance inflation factors", 0,
                         &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_gen_real_mat_print (x04cac).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

 END:
  if (b) NAG_FREE(b);
  if (perr) NAG_FREE(perr);
  if (res) NAG_FREE(res);
  if (vif) NAG_FREE(vif);
  if (x) NAG_FREE(x);
  if (y) NAG_FREE(y);
  if (isx) NAG_FREE(isx);

  return exit_status;
}