/* nag_1d_spline_fit_knots (e02bac) Example Program.

 *

 * Copyright 1996 Numerical Algorithms Group.

 *

 * Mark 4, 1996.

 *

 * 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, ncap, ncap7, r, wght;

  NagError fail;

  Nag_Spline spline;

  double fit, ss, *weights=0, *x=0, xarg, *y=0;



  INIT_FAIL(fail);



  Vprintf("nag_1d_spline_fit_knots (e02bac) Example Program Results\n");

  Vscanf("%*[^\n]");  /* Skip heading in data file */

  while(scanf("%ld",&m) != EOF)

    {

      if (m>=4)

        {

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

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

               !( y = NAG_ALLOC(m, double))

               )

            {

              Vprintf("Allocation failure\n");

              exit_status = -1;

              goto END;

            }

        }

      else

        {

          Vprintf("Invalid m.\n");

          exit_status = 1;

          return exit_status;

        }

      Vscanf("%ld%ld",&ncap,&wght);

      if (ncap>0)

        {

          ncap7 = ncap+7;

          spline.n = ncap7;

          if ( !( spline.lamda = NAG_ALLOC(ncap7, double)) )

            {

              Vprintf("Allocation failure\n");

              exit_status = -1;

              goto END;

            }

        }

      else

        {

          Vprintf("Invalid ncap.\n");

          exit_status = 1;

          return exit_status;

        }

      for (j=4; j<ncap+3; ++j)

        Vscanf("%lf",&(spline.lamda[j]));

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

        {

          if (wght==1)

            {

              Vscanf("%lf%lf",&x[r],&y[r]);

              weights[r] = 1.0;

            }

          else

            Vscanf("%lf%lf%lf",&x[r],&y[r],&weights[r]);

        }

      /* nag_1d_spline_fit_knots (e02bac).

       * Least-squares curve cubic spline fit (including

       * interpolation), one variable

       */

      nag_1d_spline_fit_knots(m, x, y, weights, &ss, &spline, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_1d_spline_fit_knots (e02bac).\n%s\n",

                  fail.message);

          exit_status = 1;

          goto END;

        }



      Vprintf("\nNumber of distinct knots = %ld\n\n", ncap+1);

      Vprintf("Distinct knots located at \n\n");

      for (j=3; j<ncap+4; j++)

        Vprintf("%8.4f%s",spline.lamda[j],

                (j-3)%6==5 || j==ncap+3 ? "\n" : " ");

      Vprintf("\n\n    J    B-spline coeff c\n\n");

      for (j=0;  j<ncap+3; ++j)

        Vprintf("    %ld  %13.4f\n",j+1,spline.c[j]);

      Vprintf("\nResidual sum of squares = ");

      Vprintf("%9.2e\n\n",ss);

      Vprintf("Cubic spline approximation and residuals\n");

      Vprintf("  r        Abscissa        Weight       Ordinate"

              "          Spline      Residual\n\n");

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

        {

          /* nag_1d_spline_evaluate (e02bbc).

           * Evaluation of fitted cubic spline, function only

           */

          nag_1d_spline_evaluate(x[r], &fit, &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;

            }



          Vprintf("%3ld    %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_1d_spline_evaluate (e02bbc), see above. */

              nag_1d_spline_evaluate(xarg, &fit, &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;

                }

              Vprintf("    %14.4f             %33.4f\n",xarg,fit);

            }

        }

      /* Free memory used by spline */

      NAG_FREE(spline.lamda);

      NAG_FREE(spline.c);

    END:

      if (weights) NAG_FREE(weights);

      if (x) NAG_FREE(x);

      if (y) NAG_FREE(y);

    }

  return exit_status;

}