/* nag_correg_ridge_opt (g02kac) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.1, 2023.
 */
/* Pre-processor includes */
#include <math.h>
#include <nag.h>
#include <stdio.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[40], sorig[40], soptloo[40];
  /*NAG Types */
  Nag_OrderType order;
  Nag_PredictError opt;
  Nag_EstimatesOption orig;
  Nag_OptionLOO optloo;
  NagError fail;

  INIT_FAIL(fail);

  printf("%s\n", "nag_correg_ridge_opt (g02kac) Example Program Results");
  /* Skip heading in data file */
  scanf("%*[^\n] ");
  /* Read in data and check array limits */
  scanf("%" NAG_IFMT "%" NAG_IFMT "%lf%39s %lf%" NAG_IFMT "%39s %39s%*[^\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("%" NAG_IFMT " ", &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_correg_ridge_opt (g02kac)
   * Ridge regression
   */
  nag_correg_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_correg_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  %5" NAG_IFMT "\n", "Degrees of freedom: ", df);
  printf("%s   %10.4f\n", "Number of effective parameters:", nep);
  printf("\n");
  ip1 = ip + 1;
  /*
   * nag_file_print_matrix_real_gen (x04cac)
   * Print real general matrix (easy-to-use)
   */
  fflush(stdout);
  nag_file_print_matrix_real_gen(order, Nag_GeneralMatrix, Nag_NonUnitDiag, ip1,
                                 one, b, pdb, "Parameter estimates", 0, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_file_print_matrix_real_gen (x04cac).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

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

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

  return exit_status;
}