#! /usr/bin/env python
#
def pyramid01_monte_carlo_test ( ):

#*****************************************************************************80
#
## PYRAMID01_MONTE_CARLO_TEST estimates some integrals.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    21 January 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from monomial_value import monomial_value
  from pyramid01_integral import pyramid01_integral
  from pyramid01_sample import pyramid01_sample
  from pyramid01_volume import pyramid01_volume

  m = 3
  test_num = 10

  e_test = np.array ( [ \
    [ 0, 0, 2, 0, 0, 2, 0, 0, 2, 2 ], \
    [ 0, 0, 0, 2, 0, 0, 2, 0, 2, 0 ], \
    [ 0, 1, 0, 0, 2, 1, 1, 3, 0, 2 ] ] )

  print ( '' )
  print ( 'PYRAMID01_MONTE_CARLO_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Use PYRAMID01_SAMPLE to estimate integrals' )
  print ( '  over the interior of the unit pyramid in 3D.' )

  seed = 123456789

  print ( '' )
  print ( '         N' ),
  print ( '        1' ),
  print ( '               Z' ),
  print ( '             X^2' ),
  print ( '             Y^2' ),
  print ( '             Z^2' ),
  print ( '            X^2Z' ),
  print ( '            Y^2Z' ),
  print ( '             Z^3' ),
  print ( '          X^2Y^2' ),
  print ( '          X^2Z^2' )
  print ( '' )

  e = np.zeros ( m, dtype = np.int32 )
  result = np.zeros ( test_num )

  n = 1

  while ( n <= 65536 ):

    x, seed = pyramid01_sample ( n, seed )

    for j in range ( 0, test_num ):

      for i in range ( 0, m ):
        e[i] = e_test[i,j]

      value = monomial_value ( m, n, e, x )

      result[j] = pyramid01_volume ( ) * np.sum ( value ) / float ( n )

    print ( '  %8d' % ( n ) ),
    for i in range ( 0, 10 ):
      print ( '  %14.6g' % ( result[i] ) ),
    print ( '' )

    n = 2 * n

  print ( '' )

  for j in range ( 0, 10 ):

    for i in range ( 0, m ):
      e[i] = e_test[i,j]

    result[j] = pyramid01_integral ( e )

  print ( '     Exact' ),
  for i in range ( 0, 10 ):
    print ( '  %14.6g' % ( result[i] ) ),
  print ( '' )
#
#  Terminate.
#
  print ( '' )
  print ( 'PYRAMID01_MONTE_CARLO_TEST' )
  print ( '  Normal end of execution.' )
  return

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