/* nag_1d_spline_fit (e02bec) Example Program.

 *

 * Copyright 1991 Numerical Algorithms Group.

 *

 * Mark 2, 1991.

 *

 * Mark 6 revised, 2000.

 * Mark 8 revised, 2004.

 */



#include <nag.h>

#include <stdio.h>

#include <nag_stdlib.h>

#include <nage02.h>



int main(void)

{

  Integer exit_status=0, j, m, nest, r;

  NagError fail;

  Nag_Comm warmstartinf;

  Nag_Spline spline;

  Nag_Start start;

  double fp, s, *sp=0, txr, *weights=0, *x=0, *y=0;



  INIT_FAIL(fail);



  Vprintf("nag_1d_spline_fit (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);

  nest = m + 4;

  if (m>=4)

    {

      if ( !( weights = NAG_ALLOC(m, double)) ||

           !( x = NAG_ALLOC(m, double)) ||

           !( y = NAG_ALLOC(m, double)) ||

           !( sp = NAG_ALLOC(2*m-1, double)) )

        {

          Vprintf("Allocation failure\n");

          exit_status = -1;

          goto END;

        }

    }

  else

    {

      Vprintf("Invalid m.\n");

      exit_status = 1;

      return exit_status;

    }

  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. */

      /* nag_1d_spline_fit (e02bec).

       * Least-squares cubic spline curve fit, automatic knot

       * placement, one variable

       */

      nag_1d_spline_fit(start, m, x, y, weights, s, nest, &fp,

                        &warmstartinf, &spline, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_1d_spline_fit (e02bec).\n%s\n", fail.message);

          exit_status = 1;

          goto END;

        }



      /* Evaluate the spline at each x point and midway

       * between x points, saving the results in sp.

       */

      for (r=0; r<m; r++)

        {

          /* nag_1d_spline_evaluate (e02bbc).

           * Evaluation of fitted cubic spline, function only

           */

          nag_1d_spline_evaluate(x[r], &sp[r*2], &spline, &fail);

          if (fail.code != NE_NOERROR)

            {

              Vprintf("Error from nag_1d_spline_fit (e02bec).\n%s\n",

                      fail.message);

              exit_status = 1;

              goto END;

            }

        }



      for (r=0; r<m-1; r++)

        {

          txr = (x[r] + x[r+1]) / 2;

          /* nag_1d_spline_evaluate (e02bbc), see above. */

          nag_1d_spline_evaluate(txr, &sp[r*2], &spline, &fail);

          if (fail.code != NE_NOERROR)

            {

              Vprintf("Error from nag_1d_spline_evaluate (e02bbc).\n%s\n",

                      fail.message);

              exit_status = 1;

              goto END;

            }

        }

      /* 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 */

 END:

  NAG_FREE(spline.lamda);

  NAG_FREE(spline.c);

  NAG_FREE(warmstartinf.nag_w);

  NAG_FREE(warmstartinf.nag_iw);

  if (weights) NAG_FREE(weights);

  if (x) NAG_FREE(x);

  if (y) NAG_FREE(y);

  if (sp) NAG_FREE(sp);

  return exit_status;

}