#! /usr/bin/env python
#
def triangle_monte_carlo ( t, n, triangle_integrand, seed ):

#*****************************************************************************80
#
## TRIANGLE_MONTE_CARLO applies the Monte Carlo rule to integrate a function.
#
#  Discussion:
#
#    The function f(x,y) is to be integrated over a triangle T.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    12 February 2010
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, real T(2,3), the triangle vertices.
#
#    Input, integer N, the number of sample points.
#
#    Input, external TRIANGLE_INTEGRAND, the integrand routine.
#
#    Input/output, integer SEED, a seed for the random
#    number generator.
#
#    Output, real RESULT, the approximate integral.
#
  import numpy as np
  from reference_to_physical_t3 import reference_to_physical_t3
  from triangle_area import triangle_area
  from triangle01_sample import triangle01_sample

  area = triangle_area ( t )

  p, seed = triangle01_sample ( n, seed )

  p2 = reference_to_physical_t3 ( t, n, p )

  fp = triangle_integrand ( p2 )
  
  result = area * np.sum ( fp[:] ) / float ( n )

  return result, seed

def triangle_monte_carlo_test ( ):

#*****************************************************************************80
#
## TRIANGLE_MONTE_CARLO_TEST estimates integrals over a general triangle.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    18 July 2018
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  from r8mat_transpose_print import r8mat_transpose_print

  global e

  m = 2

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

  print ( '' )
  print ( 'TRIANGLE_MONTE_CARLO_TEST' )
  print ( '  TRIANGLE_MONTE_CARLO estimates an integral over' )
  print ( '  a general triangle using the Monte Carlo method.' )

  t = np.array ( [ \
    [ 2.0, 3.0, 0.0 ], \
    [ 0.0, 4.0, 3.0 ] ] )

  r8mat_transpose_print ( 2, 3, t, '  Triangle vertices:' )

  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 = triangle_monte_carlo ( t, n, triangle_integrand, seed )

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

    print ( '' )

    n = 2 * n
#
#  Terminate.
#
  print ( '' )
  print ( 'TRIANGLE_MONTE_CARLO_TEST:' )
  print ( '  Normal end of execution.' )
  return

def triangle_integrand ( xy ):

#*****************************************************************************80
#
## TRIANGLE_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 ( )
  triangle_monte_carlo_test ( )
  timestamp ( )