#! /usr/bin/env python
#
def triangle01_monte_carlo ( xy_num, triangle01_integrand, seed ):

#*****************************************************************************80
#
## TRIANGLE01_MONTE_CARLO applies the Monte Carlo rule to integrate a function.
#
#  Discussion:
#
#    The function f(x,y) is to be integrated over the unit triangle.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    04 November 2016
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer XY_NUM, the number of sample points.
#
#    Input, external TRIANGLE01_INTEGRAND, the integrand routine.
#
#    Input/output, integer SEED, a seed for the random
#    number generator.
#
#    Output, real RESULT(F_NUM), the approximate integrals.
#
  import numpy as np
  from triangle01_sample import triangle01_sample

  area = 0.5

  xy, seed = triangle01_sample ( xy_num, seed )

  fxy = triangle01_integrand ( xy )
 
  result = area * np.sum ( fxy ) / float ( xy_num )

  return result, seed

def triangle01_monte_carlo_test ( ):

#*****************************************************************************80
#
## TRIANGLE01_MONTE_CARLO_TEST estimates integrals over the unit triangle.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    04 November 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  from triangle01_monomial_integral import triangle01_monomial_integral

  global e

  m = 2

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

  print ( '' )
  print ( 'TRIANGLE01_MONTE_CARLO_TEST01' )
  print ( '  TRIANGLE01_MONTE_CARLO estimates an integral over' )
  print ( '  the unit triangle using the Monte Carlo method.' )

  seed = 123456789

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

  n = 1
  e = np.zeros ( m, dtype = np.int32 )

  while ( n <= 65536 ):

    print ( '  %8d' % ( n ) ),

    for j in range ( 0, 7 ):

      e[0:m] = e_test[j,0:m]

      result, seed = triangle01_monte_carlo ( n, triangle01_integrand, seed )

      print ( '  %14.6g' % ( result ) ),

    print ( '' )

    n = 2 * n

  print ( '' )
  print ( '     Exact' ),

  for j in range ( 0, 7 ):

    e[0:m] = e_test[j,0:m]

    result = triangle01_monomial_integral ( e[0], e[1] )

    print ( '  %14.6g' % ( result ) ),

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

def triangle01_integrand ( xy ):

#*****************************************************************************80
#
## TRIANGLE01_INTEGRAND is a sample integrand function.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    04 November 2016
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, real XY(XY_NUM), the number of evaluation points.
#
#    Output, real FXY(XY_NUM), the function values.
#
  from monomial_value import monomial_value

  m = 2
  n = xy.shape[1]

  fxy = monomial_value ( m, n, e, xy )

  return fxy

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