#! /usr/bin/env python3
#
def walker_verify ( n, x, y, a ):

#****************************************************************************80
#
#  Purpose:
#
#    WALKER_VERIFY verifies a Walker Sampler structure.
#
#  Discussion:
#
#    This test applies the sampling algorithms to a Zipfian distribution.
#
#  Modified:
#
#    20 February 2016
#
#  Author:
#
#    Original C version by Warren Smith.
#    Python version by John Burkardt.
#
#  Parameters:
#
#    Input, unsigned int N, indicates the size of X.
#
#    Input, double X[N+2], contains in X[1] through X[N] the
#    probabilities of outcomes 1 through N.
#
#    Input, double Y[N+2], the Walker threshold vector.
#
#    Input, unsigned int A[N+2], the Walker index vector.
#
  import numpy as np

  z = np.zeros ( n + 2, dtype = np.float64 )
#
#  Reverse the scaling.
#
  for i in range ( 0, n + 2 ):
    z[i] = y[i] / float ( n )
#
#  Add back the adjustments.
#
  for i in range ( 1, n + 1 ):
    z[a[i]] = z[a[i]] + ( 1.0 - y[i] ) / float ( n )
#
#  Check for discrepancies between Z and X.
#
  v = 0.0
  for i in range ( 1, n + 1 ):
    v = v + abs ( z[i] - x[i] )

  return v

def walker_verify_test ( ):

#*****************************************************************************80
#
## WALKER_VERIFY_TEST tests WALKER_VERIFY.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    20 February 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from r8vec_print import r8vec_print
  from walker_build import walker_build

  print ( '' )
  print ( 'WALKER_VERIFY_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  WALKER_VERIFY verifies the Walker sampler data vectors Y and A,' )
  print ( '  for a given probability vector X.' )

  n = 9
  x = np.zeros ( n + 2, dtype = np.float64 )
  for i in range ( 1, n + 1 ):
    x[i] = np.log ( 1.0 + 1.0 / float ( i ) ) / np.log ( float ( n + 1 ) )
  r8vec_print ( n + 2, x, '  Benford PDF (ignore first and last entries):' )

  y, a = walker_build ( n, x )

  print ( '' )
  print ( '   I    A[I]    Y[i] (ignore first and last entries)' )
  print ( '' )
  for i in range ( 0, n + 2 ):
    print ( '  %2d  %2d  %10.4g' % ( i, a[i], y[i] ) )

  v = walker_verify ( n, x, y, a )

  print ( '' )
  print ( '  The verification sum = %g' % ( v ) )
  print ( '  It should be very close to zero.' )
#
#  Terminate.
#
  print ( '' )
  print ( 'WALKER_VERIFY_TEST' )
  print ( '  Normal end of execution.' )
  return

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