#! /usr/bin/env python
#
def shepard_value_1d ( nd, xd, yd, p, ni, xi ):

#*****************************************************************************80
#
## SHEPARD_VALUE_1D evaluates a 1D Shepard interpolant.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    03 July 2015
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Donald Shepard,
#    A two-dimensional interpolation function for irregularly spaced data,
#    ACM '68: Proceedings of the 1968 23rd ACM National Conference,
#    ACM, pages 517-524, 1969.
#
#  Parameters:
#
#    Input, integer ND, the number of data points.
#
#    Input, real XD(ND), the data points.
#
#    Input, real YD(ND), the data values.
#
#    Input, real P, the power.
#
#    Input, integer NI, the number of interpolation points.
#
#    Input, real XI(NI), the interpolation points.
#
#    Output, real YI(NI), the interpolated values.
# 
  import numpy as np

  yi = np.zeros ( ni )

  w = np.zeros ( nd )

  for i in range ( 0, ni ):

    if ( p == 0.0 ):

      for j in range ( 0, nd ):
        w[j] = 1.0 / float ( nd )

    else:

      z = -1
      for j in range ( 0, nd ):
        w[j] = abs ( xi[i] - xd[j] )
        if ( w[j] == 0.0 ):
          z = j

      if ( z != -1 ):
        for j in range ( 0, nd ):
          w[j] = 0.0
        w[z] = 1.0
      else:
        for j in range ( 0, nd ):
          w[j] = 1.0 / w[j] ** p
        s = np.sum ( w )
        for j in range ( 0, nd ):
          w[j] = w[j] / s
    
    for j in range ( 0, nd ):
      yi[i] = yi[i] + w[j] * yd[j]
 
  return yi

def shepard_value_1d_test ( ):

#*****************************************************************************80
#
## SHEPARD_VALUE_1D_TEST tests SHEPARD_VALUE_1D.
#
#  Discussion:
#
#    f(x) = x^3 - 12 x^2 + 39 x - 28 = ( x - 1 ) * ( x - 4 ) * ( x - 7 )
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified: 
#
#    03 July 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from r8vec2_print import r8vec2_print

  nd = 4
  ni = 21

  xd = np.array ( [ 0.0, 2.0, 5.0, 10.0 ] )
  yd = np.array ( [ -28.0, +10.0, -8.0, +162.0 ] )
 
  print ( '' )
  print ( 'SHEPARD_VALUE_1D_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  SHEPARD_VALUE_1D evaluates a Shepard 1D interpolant.' )

  p = 2.0
  print ( '' )
  print ( '  Using power P = %g' % ( p ) )

  x_min = 0.0
  x_max = 10.0
  xi = np.linspace ( x_min, x_max, ni )

  yi = shepard_value_1d ( nd, xd, yd, p, ni, xi )

  r8vec2_print ( ni, xi, yi, '  Table of interpolant values:' )
#
#  Terminate.
#
  print ( '' )
  print ( 'SHEPARD_VALUE_1D_TEST:' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  shepard_value_1d_test ( )
  timestamp ( )