/* nag_1d_cheb_fit_constr (e02agc) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.2, 2017.
 */

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

int main(void)
{
  /* Scalars */
  double fiti, xmax, xmin;
  Integer exit_status, i, iy, j, k, h, m, mf, n, pda, stride;
  NagError fail;
  Nag_OrderType order;

  /* Arrays */
  double *a = 0, *s = 0, *w = 0, *resid = 0, *x = 0, *xf = 0, *y = 0, *yf = 0;
  Integer *p = 0;

#ifdef NAG_COLUMN_MAJOR
#define A(I, J) a[(J-1)*pda + I - 1]
  order = Nag_ColMajor;
#else
#define A(I, J) a[(I-1)*pda + J - 1]
  order = Nag_RowMajor;
#endif

  INIT_FAIL(fail);

  exit_status = 0;
  printf("nag_1d_cheb_fit_constr (e02agc) Example Program Results\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  while (scanf("%" NAG_IFMT "%*[^\n] ", &mf) != EOF)
  {
    if (mf > 0) {
      /* Allocate memory for p and xf. */
      if (!(p = NAG_ALLOC(mf, Integer)) || !(xf = NAG_ALLOC(mf, double)))
      {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }

      /* Read p, xf and yf arrays */
      iy = 1;
      n = mf;
      for (i = 0; i < mf; ++i) {
        scanf("%" NAG_IFMT "%lf", &p[i], &xf[i]);
        h = iy + p[i] + 1;
        /* We need to extend array yf as we go along */
        if (!(yf = NAG_REALLOC(yf, h - 1, double)))
        {
          printf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
        for (j = iy - 1; j < h - 1; ++j)
          scanf("%lf", &yf[j]);
        scanf("%*[^\n] ");
        n += p[i];
        iy = h;
      }
      scanf("%" NAG_IFMT "%*[^\n] ", &m);

      if (m > 0) {
        /* Allocate memory for x, y and w. */
        if (!(x = NAG_ALLOC(m, double)) ||
            !(y = NAG_ALLOC(m, double)) || !(w = NAG_ALLOC(m, double)))
        {
          printf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
        for (i = 0; i < m; ++i)
          scanf("%lf%lf%lf", &x[i], &y[i], &w[i]);
        scanf("%*[^\n] ");
        scanf("%" NAG_IFMT "%lf%lf%*[^\n] ", &k, &xmin, &xmax);
        pda = k + 1;

        /* Allocate arrays a, s and resid */
        if (!(a = NAG_ALLOC((k + 1) * (k + 1), double)) ||
            !(s = NAG_ALLOC((k + 1), double)) ||
            !(resid = NAG_ALLOC(m, double)))
        {
          printf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }

        /* nag_1d_cheb_fit_constr (e02agc).
         * Least squares polynomial fit, values and derivatives may
         * be constrained, arbitrary data points
         */
        nag_1d_cheb_fit_constr(order, m, k, xmin, xmax, x, y, w, mf, xf,
                               yf, p, a, s, &n, resid, &fail);
        if (fail.code != NE_NOERROR) {
          printf("Error from nag_1d_cheb_fit_constr (e02agc).\n%s\n",
                 fail.message);
          exit_status = 1;
          goto END;
        }

        printf("\n");
        printf("Degree  RMS residual\n");
        for (i = n; i <= k; ++i)
          printf("%4" NAG_IFMT "%15.2e\n", i, s[i]);
        printf("\n");

        printf("Details of the fit of degree %2" NAG_IFMT "\n", k);
        printf("\n");
        printf("  Index   Coefficient\n");
        for (i = 0; i < k + 1; ++i)
          printf("%6" NAG_IFMT "%11.4f\n", i, A(k + 1, i + 1));
        printf("\n");

        printf("     i      x(i)       y(i)       Fit     Residual\n");
        for (i = 0; i < m; ++i) {
          /* Note that the coefficients of polynomial are stored in the
           * rows of A hence when the storage is in Nag_ColMajor order
           * then stride is the first dimension of A, k + 1.
           * When the storage is in Nag_RowMajor order then stride is 1.
           */
#ifdef NAG_COLUMN_MAJOR
          stride = k + 1;
#else
          stride = 1;
#endif
          /* nag_1d_cheb_eval2 (e02akc).
           * Evaluation of fitted polynomial in one variable from
           * Chebyshev series form
           */
          nag_1d_cheb_eval2(k, xmin, xmax, &A(k + 1, 1), stride, x[i],
                            &fiti, &fail);
          if (fail.code != NE_NOERROR) {
            printf("Error from nag_1d_cheb_eval2 (e02akc).\n%s\n",
                   fail.message);
            exit_status = 1;
            goto END;
          }
          printf("%6" NAG_IFMT "%11.4f%11.4f%11.4f", i, x[i], y[i], fiti);
          printf("%11.2e\n", fiti - y[i]);
        }
      }
    }
  }

END:
  NAG_FREE(a);
  NAG_FREE(s);
  NAG_FREE(w);
  NAG_FREE(resid);
  NAG_FREE(x);
  NAG_FREE(xf);
  NAG_FREE(y);
  NAG_FREE(yf);
  NAG_FREE(p);

  return exit_status;
}