#! /usr/bin/env python
#
def pyramid01_monomial_integral ( expon ):

#*****************************************************************************80
#
## PYRAMID01_MONOMIAL_INTEGRAL: monomial integral in a unit pyramid.
#
#  Discussion:
#
#    This routine returns the integral of
#
#      product ( 1 <= I <= 3 ) X(I)^EXPON(I)
#
#    over the unit pyramid.
#
#    The unit pyramid is defined as:
#
#    - ( 1 - Z ) <= X <= 1 - Z
#    - ( 1 - Z ) <= Y <= 1 - Z
#              0 <= Z <= 1.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    22 June 2015
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Arthur Stroud,
#    Approximate Calculation of Multiple Integrals,
#    Prentice Hall, 1971,
#    ISBN: 0130438936,
#    LC: QA311.S85.
#
#  Parameters:
#
#    Input, integer EXPON(3), the exponents.
#
#    Output, real VALUE, the integral of the monomial.
#
  from r8_choose import r8_choose
  from r8_mop import r8_mop

  if ( ( ( expon[0] % 2 ) == 0 ) and ( ( expon[1] % 2 ) == 0 ) ):

    i_hi = 2 + expon[0] + expon[1]

    value = 0.0
    for i in range ( 0, i_hi + 1 ):
      value = value + r8_mop ( i ) * r8_choose ( i_hi, i ) / float ( i + expon[2] + 1 )

    value = value * 2.0 / float ( expon[0] + 1 ) * 2.0 / float ( expon[1] + 1 )

  else:

    value = 0.0

  return value

def pyramid01_monomial_integral_test ( ):

#*****************************************************************************80
#
## PYRAMID01_MONOMIAL_INTEGRAL_TEST tests PYRAMID01_MONOMIAL_INTEGRAL.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    22 June 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from monomial_value import monomial_value
  from pyramid01_sample import pyramid01_sample
  from pyramid01_volume import pyramid01_volume

  m = 3
  n = 500000
  e_max = 6

  print ( '' )
  print ( 'PYRAMID01_MONOMIAL_INTEGRAL_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  PYRAMID01_MONOMIAL_INTEGRAL returns the integral of a monomial' )
  print ( '  over the unit pyramid in 3D.' )
  print ( '  Compare to a Monte Carlo estimate.' )
#
#  Get sample points.
#
  seed = 123456789
  x, seed = pyramid01_sample ( n, seed )

  print ( '' )
  print ( '  Number of sample points used is %d' % ( n ) )
  print ( '' )
  print ( '   E1  E2  E3     MC-Estimate      Exact           Error' )
  print ( '' )
#
#  Check all monomials, with only even dependence on X or Y, 
#  up to total degree E_MAX.
#
  e = np.zeros ( 3, dtype = np.int32 )

  for e3 in range ( 0, e_max + 1 ):
    e[2] = e3
    for e2 in range ( 0, e_max - e3 + 1, 2 ):
      e[1] = e2
      for e1 in range ( 0 , e_max - e3 - e2 + 1, 2 ):
        e[0] = e1

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

        q = pyramid01_volume ( ) * np.sum ( value ) / float ( n )
        exact = pyramid01_monomial_integral ( e )
        error = abs ( q - exact )

        print ( '  %2d  %2d  %2d  %14.6g  %14.6g  %10.2g' \
          % ( e[0], e[1], e[2], q, exact, error ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'PYRAMID01_MONOMIAL_INTEGRAL_TEST:' )
  print ( '  Normal end of execution.' )
  return

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