# include <stdbool.h>
# include <stdlib.h>
# include <stdio.h>
# include <math.h>

# include "pwl_approx_1d.h"
# include "qr_solve.h"
# include "r8lib.h"

/******************************************************************************/

double *pwl_approx_1d ( int nd, double xd[], double yd[], int nc, double xc[] )

/******************************************************************************/
/*
  Purpose:

    PWL_APPROX_1D determines the control values for a PWL approximant.

  Discussion:

    The piecewise linear approximant is defined by NC control pairs 
    (XC(I),YC(I)) and approximates ND data pairs (XD(I),YD(I)).

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    10 October 2012

  Author:

    John Burkardt

  Parameters:

    Input, int ND, the number of data points.
    ND must be at least 1.

    Input, double XD[ND], the data points.

    Input, double YD[ND], the data values.

    Input, int NC, the number of control points.
    NC must be at least 1.

    Input, double XC[NC], the control points.  Set these with a
    command like 
      xc = r8vec_linspace_new ( nc, xmin, xmax );

    Output, double PWL_APPROX_1D[NC], the control values.
*/
{
  double *a;
  double *yc;
/*
  Define the NDxNC linear system that determines the control values.
*/
  a = pwl_approx_1d_matrix ( nd, xd, yd, nc, xc );
/*
  Solve the system.
*/
  yc = qr_solve ( nd, nc, a, yd );

  free ( a );

  return yc;
}
/******************************************************************************/

double *pwl_approx_1d_matrix ( int nd, double xd[], double yd[], int nc, 
  double xc[] )

/******************************************************************************/
/*
  Purpose:

    PWL_APPROX_1D_MATRIX returns the matrix for the PWL approximant controls.

  Discussion:

    The value of the piecewise linear approximant, using control points XC
    and control values YC, evaluated at the point XD, can be represented by

      YD = A * YC

    where A is a matrix whose values depend on XC and XD.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    09 October 2012

  Author:

    John Burkardt

  Parameters:

    Input, int ND, the number of data points.
    ND must be at least 1.

    Input, double XD[ND], the data points.

    Input, double YD[ND], the data values.

    Input, int NC, the number of control points.
    NC must be at least 1.

    Input, double XC[NC], the control points.

    Output, double PWL_APPROX_1D_MATRIX[ND*NC], the matrix.
*/
{
  double *a;
  int i;
  int j;
  int k;
  double t;

  a = ( double * ) malloc ( nd * nc * sizeof ( double ) );

  for ( j = 0; j < nc; j++ )
  {
    for ( i = 0; i < nd; i++ )
    {
      a[i+j*nd] = 0.0;
    }
  }

  for ( i = 0; i < nd; i++ )
  {
    k = nc - 2;
    for ( j = 1; j < nc - 1; j++ )
    {
      if ( xd[i] < xc[j] )
      {
        k = j - 1;
        break;
      }
    }
    t = ( xd[i] - xc[k] ) / ( xc[k+1] - xc[k] );
    a[i+k*nd]     = 1.0 - t;
    a[i+(k+1)*nd] =       t;
  }

  return a;
}
/******************************************************************************/

double *pwl_interp_1d ( int nd, double xd[], double yd[], int ni, double xi[] )

/******************************************************************************/
/*
  Purpose:

    PWL_INTERP_1D evaluates the piecewise linear interpolant.

  Discussion:

    The piecewise linear interpolant L(ND,XD,YD)(X) is the piecewise
    linear function which interpolates the data (XD(I),YD(I)) for I = 1
    to ND.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    09 October 2012

  Author:

    John Burkardt

  Parameters:

    Input, int ND, the number of data points.
    ND must be at least 1.

    Input, double XD[ND], the data points.

    Input, double YD[ND], the data values.

    Input, int NI, the number of interpolation points.

    Input, double XI[NI], the interpolation points.

    Output, double PWL_INTERP_1D[NI], the interpolated values.
*/
{
  int i;
  int k;
  double t;
  double *yi;

  yi = ( double * ) malloc ( ni * sizeof ( double ) );

  if ( nd == 1 )
  {
    for ( i = 0; i < ni; i++ )
    {
      yi[i] = yd[0];
    }
    return yi;
  }

  for ( i = 0; i < ni; i++ )
  {
    if ( xi[i] <= xd[0] )
    {
      t = ( xi[i] - xd[0] ) / ( xd[1] - xd[0] );
      yi[i] = ( 1.0 - t ) * yd[0] + t * yd[1];
    }
    else if ( xd[nd-1] <= xi[i] )
    {
      t = ( xi[i] - xd[nd-2] ) / ( xd[nd-1] - xd[nd-2] );
      yi[i] = ( 1.0 - t ) * yd[nd-2] + t * yd[nd-1];
    }
    else
    {
      for ( k = 1; k < nd; k++ )
      {
        if ( xd[k-1] <= xi[i] && xi[i] <= xd[k] )
        {
          t = ( xi[i] - xd[k-1] ) / ( xd[k] - xd[k-1] );
          yi[i] = ( 1.0 - t ) * yd[k-1] + t * yd[k];
          break;
        }
      }
    }
  }
  
  return yi;
}