/* nag_1d_spline_fit(e02bec) Example Program
 *
 * Copyright 1991 Numerical Algorithms Group.
 *
 * Mark 2, 1991.
 *
 * Mark 6 revised, 2000.
 */

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

#define MMAX 50
#define NEST MMAX + 4

int main(void)
{
  Integer m, r, j;
  double weights[MMAX], x[MMAX], y[MMAX];
  double s, fp, sp[2*MMAX-1], txr;
  Nag_Start start;
  Nag_Comm warmstartinf;
  Nag_Spline spline;

  Vprintf("e02bec Example Program Results\n");
  Vscanf("%*[^\n]"); /* Skip heading in data file */
  /* Input the number of data points, followed by the data
   * points x, the function values y and the weights w.
   */
  Vscanf("%ld",&m);
  if (m>0 && m<=MMAX)
    {
      start = Nag_Cold;
      for (r=0; r<m; r++)
        Vscanf("%lf%lf%lf",&x[r], &y[r], &weights[r]);
      /* Read in successive values of s until end of data file. */
      while(scanf("%lf",&s) !=EOF)
        {
          /* Determine the spline approximation. */
          e02bec(start, m, x, y, weights, s, (Integer)(NEST), &fp,
                 &warmstartinf, &spline, NAGERR_DEFAULT);
          /* Evaluate the spline at each x point and midway
           * between x points, saving the results in sp.
           */
          for (r=0; r<m; r++)
            e02bbc(x[r], &sp[r*2], &spline, NAGERR_DEFAULT);
          for (r=0; r<m-1; r++)
            {
              txr = (x[r] + x[r+1]) / 2;
              e02bbc(txr, &sp[r*2], &spline, NAGERR_DEFAULT);
            }
          /* Output the results. */
          Vprintf("\nCalling with smoothing factor s = %12.3e\n",s);
          Vprintf("\nNumber of distinct knots = %ld\n\n", spline.n-6);
          Vprintf("Distinct knots located at \n\n");
          for (j=3; j<spline.n-3; j++)
            Vprintf("%8.4f%s",spline.lamda[j],
                    (j-3)%6==5 || j==spline.n-4 ? "\n" : " ");
          Vprintf("\n\n    J      B-spline coeff c\n\n");
          for (j=0;  j<spline.n-4; ++j)
            Vprintf("  %3ld  %13.4f\n",j+1,spline.c[j]);
          Vprintf("\nWeighted sum of squared residuals fp = %12.3e\n",fp);
          if (fp == 0.0)
            Vprintf("The spline is an interpolating spline\n");
          else if (spline.n == 8)
            Vprintf("The spline is the weighted least-squares cubic"
                    "polynomial\n");
          start = Nag_Warm;
        }
      /* Free memory allocated in spline and warmstartinf */
      NAG_FREE(spline.lamda);
      NAG_FREE(spline.c);
      NAG_FREE(warmstartinf.nag_w);
      NAG_FREE(warmstartinf.nag_iw);
      return EXIT_SUCCESS;
    }
  else
    {
      Vfprintf(stderr, "m is out of range: m = %5ld\n", m);
      return EXIT_FAILURE;
    }
}