#! /usr/bin/env python
#
def genbet ( aa, bb ):

#*****************************************************************************80
#
## GENBET generates a beta random deviate.
#
#  Discussion:
#
#    This procedure returns a single random deviate from the beta distribution
#    with parameters A and B.  The density is
#
#      x^(a-1) * (1-x)^(b-1) / Beta(a,b) for 0 < x < 1
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    03 September 2018
#
#  Author:
#
#    Original FORTRAN77 version by Barry Brown, James Lovato.
#    Python version by John Burkardt.
#
#  Reference:
#
#    Russell Cheng,
#    Generating Beta Variates with Nonintegral Shape Parameters,
#    Communications of the ACM,
#    Volume 21, Number 4, April 1978, pages 317-322.
#
#  Parameters:
#
#    Input, real AA, the first parameter of the beta distribution.
#    0.0 < AA.
#
#    Input, real BB, the second parameter of the beta distribution.
#    0.0 < BB.
#
#    Output, real VALUE, a beta random variate.
#
  import numpy as np
  from r4_exp import r4_exp
  from r4_uni_01 import r4_uni_01
  from sys import exit

  if ( aa <= 0.0 ):
    fprintf ( 1, '' )
    fprintf ( 1, 'GENBET - Fatal error!' )
    fprintf ( 1, '  AA <= 0.0' )
    exit ( 'GENBET - Fatal error!' )

  if ( bb <= 0.0 ):
    fprintf ( 1, '' )
    fprintf ( 1, 'GENBET - Fatal error!' )
    fprintf ( 1, '  BB <= 0.0' )
    exit ( 'GENBET - Fatal error!' )
#
#  Algorithm BB
#
  if ( 1.0 < aa and 1.0 < bb ):

    a = min ( aa, bb )
    b = max ( aa, bb )
    alpha = a + b
    beta = np.sqrt ( ( alpha - 2.0 ) / ( 2.0 * a * b - alpha ) )
    gamma = a + 1.0 / beta

    while ( True ):

      u1 = r4_uni_01 ( )
      u2 = r4_uni_01 ( )
      v = beta * np.log ( u1 / ( 1.0 - u1 ) )
#
#  EXP ( V ) replaced by R4_EXP ( V ), to truncate at large magnitude inputs.
#
      w = a * r4_exp ( v )

      z = u1 ** 2 * u2
      r = gamma * v - np.log ( 4.0 )
      s = a + r - w

      if ( 5.0 * z <= s + 1.0 + np.log ( 5.0 ) ):
        break

      t = np.log ( z )
      if ( t <= s ):
        break

      if ( t <= ( r + alpha * np.log ( alpha / ( b + w ) ) ) ):
        break
#
#  Algorithm BC
#
  else:

    a = max ( aa, bb )
    b = min ( aa, bb )
    alpha = a + b
    beta = 1.0 / b
    delta = 1.0 + a - b

    k1 = delta * ( 1.0 / 72.0 + b / 24.0 ) / ( a / b - 7.0 / 9.0 )
    k2 = 0.25 + ( 0.5 + 0.25 / delta ) * b

    while ( True ):

      u1 = r4_uni_01 ( )
      u2 = r4_uni_01 ( )

      if ( u1 < 0.5 ):

        y = u1 * u2
        z = u1 * y

        if ( k1 <= 0.25 * u2 + z - y ):
          continue

      else:

        z = u1 ** 2 * u2

        if ( z <= 0.25 ):

          v = beta * np.log ( u1 / ( 1.0 - u1 ) )
          w = a * np.exp ( v )

          if ( aa == a ):
            value = w / ( b + w )
          else:
            value = b / ( b + w )

          return value

        if ( k2 < z ):
          continue

      v = beta * np.log ( u1 / ( 1.0 - u1 ) )
      w = a * np.exp ( v )

      if ( np.log ( z ) <= alpha * ( np.log ( alpha / ( b + w ) ) + v ) - np.log ( 4.0 ) ):
        break

  if ( aa == a ):
    value = w / ( b + w )
  else:
    value = b / ( b + w )

  return value

def genbet_test ( phrase ):

#*****************************************************************************80
#
## GENBET_TEST tests GENBET, which generates Beta deviates.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    03 September 2019
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  from genunf import genunf
  from initialize import initialize
  from phrtsd import phrtsd
  from set_initial_seed import set_initial_seed
  from stats import stats
  from trstat import trstat

  n = 1000

  print ( '' )
  print ( 'RANLIB_TEST_GENBET' )
  print ( '  GENBET generates Beta deviates.' )
#
#  Initialize the generators.
#
  initialize ( )
#
#  Set the seeds based on the phrase.
#
  seed1, seed2 = phrtsd ( phrase )
#
#  Initialize all generators.
#
  set_initial_seed ( seed1, seed2 )
#
#  Select the parameters at random within a given range.
#
  low = 1.0
  high = 10.0
  a = genunf ( low, high )

  low = 1.0
  high = 10.0
  b = genunf ( low, high )

  print ( '' )
  print ( '  N = %d' % ( n ) )
  print ( '' )
  print ( '  Parameters:' )
  print ( '' )
  print ( '  A = %g' % ( a ) )
  print ( '  B = %g' % ( b ) )
#
#  Generate N samples.
#
  array = np.zeros ( n )
  for i in range ( 0, n ):
    array[i] = genbet ( a, b )
#
#  Compute statistics on the samples.
#
  av, var, xmin, xmax = stats ( array, n )
#
#  Request expected value of statistics for this distribution.
#
  pdf = 'bet'
  param = np.array ( [ a, b ] )

  avtr, vartr = trstat ( pdf, param )

  print ( '' )
  print ( '  Sample data range:          %14.6g  %14.6g' % ( xmin, xmax ) )
  print ( '  Sample mean, variance:      %14.6g  %14.6g' % ( av,   var ) )
  print ( '  Distribution mean, variance %14.6g  %14.6g' % ( avtr, vartr ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'GENBET_TEST' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  phrase = 'Randomizer'
  genbet_test ( phrase )
  timestamp ( )