#! /usr/bin/env python
#
def chebyshev_coef_1d ( nd, xd, yd ):

#*****************************************************************************80
#
## CHEBYSHEV_COEF_1D determines the Chebyshev interpolant coefficients.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    26 July 2017
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer ND, the number of data points.
#    ND must be at least 1.
#
#    Input, real XD(ND), the data locations.
#
#    Input, real YD(ND), the data values.
#
#    Output, real C(ND), the Chebyshev coefficients.
#
#    Output, real XMIN, XMAX, the interpolation interval.
#
  import numpy as np

  if ( nd == 1 ):
    c = 1.0
    xmin = xd[0]
    xmax = xd[0]
    return c, xmin, xmax

  xmin = min ( xd )
  xmax = max ( xd )
#
#  Map XD to [-1,+1].
#
  x = ( 2.0 * xd - xmin - xmax ) / ( xmax - xmin )
#
#  Form the Chebyshev Vandermonde matrix.
#
  a = np.outer ( np.arccos ( x ), np.arange ( 0, nd ) )
  a = np.cos ( a )
#
#  Solve for the expansion coefficients.
#
  c = np.linalg.solve ( a, yd )

  return c, xmin, xmax

def chebyshev_interp_1d ( nd, xd, yd, ni, xi ):

#*****************************************************************************80
#
## CHEBYSHEV_INTERP_1D determines and evaluates the Chebyshev interpolant.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    26 July 2017
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer ND, the number of data points.
#    ND must be at least 1.
#
#    Input, real XD(ND), the data locations.
#
#    Input, real YD(ND), the data values.
#
#    Input, integer NI, the number of interpolation points.
#
#    Input, real XI(NI), the interpolation points, which
#    must be each be in the interval [ min(XD), max(XD)].
#
#    Output, real YI(NI), the interpolated values.
#
  c, xmin, xmax = chebyshev_coef_1d ( nd, xd, yd )

  yi = chebyshev_value_1d ( nd, c, xmin, xmax, ni, xi )

  return yi

def chebyshev_interp_1d_test ( ):

#*****************************************************************************80
#
## CHEBYSHEV_INTERP_1D_TEST tests the CHEBYSHEV_INTERP_1D library.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    26 July 2017
#
#  Author:
#
#    John Burkardt
#
  import platform
  from test_interp import p00_prob_num

  print ( '' )
  print ( 'CHEBYSHEV_INTERP_1D_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Test the CHEBYSHEV_INTERP_1D library.' )

  r8vec_linspace_test ( )

  prob_num = p00_prob_num ( )
  for prob in range ( 1, prob_num + 1 ):
    chebyshev_interp_1d_test01 ( prob )
#
#  Terminate.
#
  print ( '' )
  print ( 'CHEBYSHEV_INTERP_1D_TEST:' )
  print ( '  Normal end of execution.' )
  return

def chebyshev_interp_1d_test01 ( prob ):

#*****************************************************************************80
#
## CHEBYSHEV_INTERP_1D_TEST01 tests CHEBYSHEV_VALUE_1D.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    26 July 2017
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np
  from test_interp import p00_data
  from test_interp import p00_data_num

  print ( '' )
  print ( 'CHEBYSHEV_INTERP_1D_TEST01:' )
  print ( '  Interpolate data from TEST_INTERP problem #%d.' % ( prob ) )

  nd = p00_data_num ( prob )
  print ( '  Number of data points = %d' % ( nd ) )

  xy = p00_data ( prob, 2, nd )
  
  xd = xy[0,0:nd]
  yd = xy[1,0:nd]

  r8vec2_print ( nd, xd, yd, '  X and Y data values:' )
#
#  #1:  Does interpolant match function at interpolation points?
#
  ni = nd
  xi = xd
  yi = chebyshev_interp_1d ( nd, xd, yd, ni, xi )

  int_error = np.linalg.norm ( yi - yd ) / float ( ni )

  print ( '' )
  print ( '  L2 interpolation error averaged per interpolant node = %g' % ( int_error ) )
#
#  #2: Compare estimated curve length to piecewise linear (minimal) curve length.
#  Assume data is sorted, and normalize X and Y dimensions by (XMAX-XMIN) and
#  (YMAX-YMIN).
#
  xmin = min ( xd )
  xmax = max ( xd )
  ymin = min ( yd )
  ymax = max ( yd )

  ni = 501
  xi = r8vec_linspace ( ni, xmin, xmax )
  yi = chebyshev_interp_1d ( nd, xd, yd, ni, xi )

  ld = 0.0
  for i in range ( 0, nd - 1 ):
    ld = ld + np.sqrt ( \
        ( ( xd[i+1] - xd[i] ) / ( xmax - xmin ) ) ** 2 \
      + ( ( yd[i+1] - yd[i] ) / ( ymax - ymin ) ) ** 2 )
 
  li = 0.0
  for i in range ( 0, ni - 1 ):
    li = li + np.sqrt ( \
        ( ( xi[i+1] - xi[i] ) / ( xmax - xmin ) ) ** 2 \
      + ( ( yi[i+1] - yi[i] ) / ( ymax - ymin ) ) ** 2 )

  print ( '\n' )
  print ( '  Normalized length of piecewise linear interpolant = %g' % ( ld ) )
  print ( '  Normalized length of polynomial interpolant       = %g' % ( li ) )
#
#  #3: Plot the data.
#
  plt.plot ( xd, yd, 'b-', linewidth = 3 )
  plt.plot ( xd, yd, 'k.', markersize = 20 )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'Piecewise Linear Interpolant' )
  plt.grid ( True )
  filename = 'p%02d_data.png' % ( prob )
  plt.savefig ( filename )
  print ( '' )
  print ( '  Created plot file "%s".' % ( filename ) )
  plt.show ( )
#
#  #4: Plot the piecewise linear and polynomial interpolants.
#
  ni = 101
  xmin = min ( xd )
  xmax = max ( xd )
  xi = r8vec_linspace ( ni, xmin, xmax )
  yi = chebyshev_interp_1d ( nd, xd, yd, ni, xi )
  plt.plot ( xi, yi, 'r-', linewidth = 3 )
  plt.plot ( xd, yd, 'b-', linewidth = 3 )
  plt.plot ( xd, yd, 'k.', markersize = 20 )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'Polynomial Interpolant using Lagrange' )
  plt.grid ( True )
  filename = 'p%02d_poly.png' % ( prob )
  plt.savefig ( filename )
  print ( '  Created plot file "%s".' % ( filename ) )
  plt.show ( )

  return

def chebyshev_value_1d ( nd, c, xmin, xmax, ni, xi ):

#*****************************************************************************80
#
## CHEBYSHEV_VALUE_1D evaluates a Chebyshev interpolant, given its coefficients.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    26 July 2017
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer ND, the number of data points.
#    ND must be at least 1.
#
#    Input, real C(ND), the Chebyshev coefficients.
#
#    Output, real XMIN, XMAX, the interpolation interval.
#
#    Input, integer NI, the number of interpolation points.
#
#    Input, real XI(NI), the interpolation points, which
#    must be each be in the interval [XMIN,XMAX].
#
#    Output, real YI(NI), the interpolated values.
#
  import numpy as np

  if ( nd == 1 ):
    yi = c.copy ( )
    return yi
#
#  Map XI to [-1,+1].
#
  x = ( 2.0 * xi - xmin - xmax ) / ( xmax - xmin )

  a = np.outer ( np.arccos ( x ), np.arange ( 0, nd ) )
  a = np.cos ( a )

  yi = np.dot ( a, c )

  return yi

def r8vec_linspace ( n, a, b ):

#*****************************************************************************80
#
## R8VEC_LINSPACE creates a column vector of linearly spaced values.
#
#  Discussion:
#
#    An R8VEC is a vector of R8's.
#
#    While MATLAB has the built in command 
#
#      x = linspace ( a, b, n )
#
#    that command has the distinct disadvantage of returning a ROW vector.
#
#    4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12.
#
#    In other words, the interval is divided into N-1 even subintervals,
#    and the endpoints of intervals are used as the points.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    02 January 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the number of entries in the vector.
#
#    Input, real A, B, the first and last entries.
#
#    Output, real X(N), a vector of linearly spaced data.
#
  import numpy as np

  x = np.zeros ( n )

  if ( n == 1 ):
    x[0] = ( a + b ) / 2.0
  else:
    for i in range ( 0, n ):
     x[i] = ( ( n - 1 - i ) * a \
            + (         i ) * b ) \
            / ( n - 1     )
 
  return x

def r8vec_linspace_test ( ):

#*****************************************************************************80
#
## R8VEC_LINSPACE_TEST tests R8VEC_LINSPACE.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    02 January 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'R8VEC_LINSPACE_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R8VEC_LINSPACE returns evenly spaced values between A and B.' )

  n = 5
  x_lo = 10.0
  x_hi = 20.0
  x = r8vec_linspace ( n, x_lo, x_hi )

  r8vec_print ( n, x, '  The vector:' )
#
#  Terminate.
#
  print ( '' )
  print ( 'R8VEC_LINSPACE_TEST' )
  print ( '  Normal end of execution.' )
  return

def r8vec_print ( n, a, title ):

#*****************************************************************************80
#
## R8VEC_PRINT prints an R8VEC.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    31 August 2014
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the dimension of the vector.
#
#    Input, real A(N), the vector to be printed.
#
#    Input, string TITLE, a title.
#
  print ( '' )
  print ( title )
  print ( '' )
  for i in range ( 0, n ):
    print ( '%6d:  %12g' % ( i, a[i] ) )

def r8vec_print_test ( ):

#*****************************************************************************80
#
## R8VEC_PRINT_TEST tests R8VEC_PRINT.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    29 October 2014
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'R8VEC_PRINT_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R8VEC_PRINT prints an R8VEC.' )

  n = 4
  v = np.array ( [ 123.456, 0.000005, -1.0E+06, 3.14159265 ], dtype = np.float64 )
  r8vec_print ( n, v, '  Here is an R8VEC:' )
#
#  Terminate.
#
  print ( '' )
  print ( 'R8VEC_PRINT_TEST:' )
  print ( '  Normal end of execution.' )
  return

def r8vec2_print ( n, a1, a2, title ):

#*****************************************************************************80
#
## R8VEC2_PRINT prints an R8VEC2.
#
#  Discussion:
#
#    An R8VEC2 is a dataset consisting of N pairs of real values, stored
#    as two separate vectors A1 and A2.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    27 June 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the number of components of the vector.
#
#    Input, real A1(N), A2(N), the vectors to be printed.
#
#    Input, string TITLE, a title.
#
  print ( '' )
  print ( title )
  print ( '' )
  for i in range ( 0, n ):
    print ( '  %6d:   %12g  %12g' % ( i, a1[i], a2[i] ) )

  return

def r8vec2_print_test ( ):

#*****************************************************************************80
#
## R8VEC2_PRINT_TEST tests R8VEC2_PRINT.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    27 June 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'R8VEC2_PRINT_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R8VEC2_PRINT prints a pair of R8VEC\'s.' )

  n = 6
  v = np.array ( [ 0.0, 0.20, 0.40, 0.60, 0.80, 1.0 ], dtype = np.float64 )
  w = np.array ( [ 0.0, 0.04, 0.16, 0.36, 0.64, 1.0 ], dtype = np.float64 )
  r8vec2_print ( n, v, w, '  Print a pair of R8VEC\'s:' )
#
#  Terminate.
#
  print ( '' )
  print ( 'R8VEC2_PRINT_TEST:' )
  print ( '  Normal end of execution.' )
  return

def timestamp ( ):

#*****************************************************************************80
#
## TIMESTAMP prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    06 April 2013
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    None
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return None

def timestamp_test ( ):

#*****************************************************************************80
#
## TIMESTAMP_TEST tests TIMESTAMP.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    03 December 2014
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    None
#
  import platform

  print ( '' )
  print ( 'TIMESTAMP_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  TIMESTAMP prints a timestamp of the current date and time.' )
  print ( '' )

  timestamp ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'TIMESTAMP_TEST:' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  timestamp ( ) 
  chebyshev_interp_1d_test ( )
  timestamp ( )