/* nag_fit_dim1_spline_knots (e02bac) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.2, 2023.
 *
 */

#include <nag.h>
#include <stdio.h>

int main(void) {
  Integer exit_status = 0, j, m, ncap, ncap7, r, wght;
  NagError fail;
  Nag_Spline spline;
  double fit, ss, *weights = 0, *x = 0, xarg, *y = 0;

  INIT_FAIL(fail);

  /* Initialize spline */
  spline.lamda = 0;
  spline.c = 0;

  printf("nag_fit_dim1_spline_knots (e02bac) Example Program Results\n");
  scanf("%*[^\n]"); /* Skip heading in data file */
  while (scanf("%" NAG_IFMT "", &m) != EOF)
  {
    if (m >= 4) {
      if (!(weights = NAG_ALLOC(m, double)) || !(x = NAG_ALLOC(m, double)) ||
          !(y = NAG_ALLOC(m, double))) {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }
    } else {
      printf("Invalid m.\n");
      exit_status = 1;
      goto END;
    }
    scanf("%" NAG_IFMT "%" NAG_IFMT "", &ncap, &wght);
    if (ncap > 0) {
      ncap7 = ncap + 7;
      spline.n = ncap7;
      if (!(spline.lamda = NAG_ALLOC(ncap7, double))) {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }
    } else {
      printf("Invalid ncap.\n");
      exit_status = 1;
      goto END;
    }
    for (j = 4; j < ncap + 3; ++j)
      scanf("%lf", &(spline.lamda[j]));
    for (r = 0; r < m; ++r) {
      if (wght == 1) {
        scanf("%lf%lf", &x[r], &y[r]);
        weights[r] = 1.0;
      } else
        scanf("%lf%lf%lf", &x[r], &y[r], &weights[r]);
    }
    /* nag_fit_dim1_spline_knots (e02bac).
     * Least squares curve cubic spline fit (including
     * interpolation), one variable
     */
    nag_fit_dim1_spline_knots(m, x, y, weights, &ss, &spline, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_fit_dim1_spline_knots (e02bac).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }

    printf("\nNumber of distinct knots = %" NAG_IFMT "\n\n", ncap + 1);
    printf("Distinct knots located at \n\n");
    for (j = 3; j < ncap + 4; j++)
      printf("%8.4f%s", spline.lamda[j],
             (j - 3) % 6 == 5 || j == ncap + 3 ? "\n" : " ");
    printf("\n\n    J    B-spline coeff c\n\n");
    for (j = 0; j < ncap + 3; ++j)
      printf("    %" NAG_IFMT "  %13.4f\n", j + 1, spline.c[j]);
    printf("\nResidual sum of squares = ");
    printf("%11.2e\n\n", ss);
    printf("Cubic spline approximation and residuals\n");
    printf("  r        Abscissa        Weight       Ordinate"
           "          Spline      Residual\n\n");
    for (r = 0; r < m; ++r) {
      /* nag_fit_dim1_spline_eval (e02bbc).
       * Evaluation of fitted cubic spline, function only
       */
      nag_fit_dim1_spline_eval(x[r], &fit, &spline, &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_fit_dim1_spline_eval (e02bbc).\n%s\n",
               fail.message);
        exit_status = 1;
        goto END;
      }

      printf("%3" NAG_IFMT "    %11.4f    %11.4f    %11.4f     %11.4f"
             "    %10.1e\n",
             r + 1, x[r], weights[r], y[r], fit, fit - y[r]);
      if (r < m - 1) {
        xarg = (x[r] + x[r + 1]) * 0.5;
        /* nag_fit_dim1_spline_eval (e02bbc), see above. */
        nag_fit_dim1_spline_eval(xarg, &fit, &spline, &fail);
        if (fail.code != NE_NOERROR) {
          printf("Error from nag_fit_dim1_spline_eval (e02bbc).\n%s\n",
                 fail.message);
          exit_status = 1;
          goto END;
        }
        printf("    %14.4f             %33.4f\n", xarg, fit);
      }
    }
    /* Free memory used by spline */
    NAG_FREE(spline.lamda);
    NAG_FREE(spline.c);
  END:
    NAG_FREE(weights);
    NAG_FREE(x);
    NAG_FREE(y);
  }
  return exit_status;
}