#! /usr/bin/env python3
#
def python11 ( ):

  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'python11():' )
  print ( '  Exercises for random numbers, simulation, testing.' )
#
#  Seed the default random number generator.
#
  for seed in [ 123456789, 987654321, 123456789 ]:
    print ( '' )
    print ( '  Initialize random number generator with seed = ', seed )
    rng = np.random.default_rng ( seed )
    r = rng.integers ( low = 0, high = 100, size = 10 )
    print ( '  Five random integers 0 <= i < 100:' )
    print ( r )
#
#  Return a scalar, vector, and matrix of random integers.
#
  print ( '' )
  r = rng.integers ( low = 0, high = 100 )
  print ( '  r = rng.integers ( low = 0, high = 100 )' )
  print ( r )
  r = rng.integers ( low = 0, high = 100, size = 1 )
  print ( '  r = rng.integers ( low = 0, high = 100, size = 1 )' )
  print ( r )
  r = rng.integers ( low = 0, high = 100, size = 3 )
  print ( '  r = rng.integers ( low = 0, high = 100, size = 3 )' )
  print ( r )
  r = rng.integers ( low = 0, high = 100, size = ( 3, 2 ) )
  print ( '  r = rng.integers ( low = 0, high = 100, size = ( 3, 2 ) )' )
  print ( r )
#
#  Return a scalar, vector, and matrix of random floats in [0,1].
#
  print ( '' )
  print ( '  Return random floats in [0,1]' )
  print ( '' )
  r = rng.random ( )
  print ( '  r = rng.random ( )' )
  print ( r )
  r = rng.random ( 1 )
  print ( '  rng.random ( 1 )' )
  print ( r )
  r = rng.random ( 3 )
  print ( '  rng.random ( 3 )' )
  print ( r )
  r = rng.random ( ( 3, 2 ) )
  print ( '  rng.random ( ( 3, 2 ) )' )
  print ( r )
#
#  Uniform statistics
#
  print ( '' )
  print ( '  statistics for random():' )
  print ( '' )
  x = rng.random ( 10000 )
  x_mean = np.mean ( x )
  x_std = np.std ( x )
  x_min = np.min ( x )
  x_max = np.max ( x )
  print ( '  x = rng._random ( 10000 )' )
  print ( '  x[0:5] = ', x[0:5] )
  print ( '  mean ( x ) = ', x_mean, ' Distribution mean = 0' )
  print ( '  std ( x ) = ', x_std, ' Distribution std = ', 1.0 / np.sqrt ( 12.0 ) )
  print ( '  min ( x ) = ', x_min, ' Distribution min = 0.0' )
  print ( '  min ( x ) = ', x_max, ' Distribution max = 1.0' )
#
#  Return random floats in [a,b].
#
  a = np.pi
  b = 10.0
  print ( '' )
  print ( '  Return random floats in [', a, ',', b, ']' )
  print ( '' )

  r = a + ( b - a ) * rng.random ( )
  print ( '  r = a + ( b - a ) * rng.random ( )' )
  print ( r )
  r = a + ( b - a ) * rng.random ( 1 )
  print ( '  r = a + ( b - a ) * rng.random ( 1 )' )
  print ( r )
  r = a + ( b - a ) * rng.random ( 3 )
  print ( '  r = a + ( b - a ) * rng.random ( 3 )' )
  print ( r )
  r = a + ( b - a ) * rng.random ( ( 3, 2 ) )
  print ( '  r = a + ( b - a ) * rng.random ( ( 3, 2 ) )' )
  print ( r )
#
#  Return a scalar, vector, and matrix of random standard normals.
#
  print ( '' )
  r = rng.standard_normal ( )
  print ( '  r = rng.standard_normal ( )' )
  print ( r )
  r = rng.standard_normal ( 1 )
  print ( '  rng.standard_normal ( 1 )' )
  print ( r )
  r = rng.standard_normal ( 3 )
  print ( '  rng.standard_normal ( 3 )' )
  print ( r )
  r = rng.standard_normal ( ( 3, 2 ) )
  print ( '  rng.standard_normal ( ( 3, 2 ) )' )
  print ( r )
#
#  Standard normal statistics
#
  print ( '' )
  print ( '  standard_normal:' )
  print ( '' )
  x = rng.standard_normal ( 10000 )
  x_mean = np.mean ( x )
  x_std = np.std ( x )
  x_min = np.min ( x )
  x_max = np.max ( x )
  print ( '  x = rng._standard_normal ( 10000 )' )
  print ( '  x[0:5] = ', x[0:5] )
  print ( '  mean ( x ) = ', x_mean, ' Distribution mean = 0.0' )
  print ( '  std ( x ) = ', x_std, ' Distribution std = 1.0' )
  print ( '  min ( x ) = ', x_min, ' Distribution min = -inf' )
  print ( '  max ( x ) = ', x_max, ' Distribution max = +inf' )
#
#  Normal statistics
#
  print ( '' )
  print ( '  normal:' )
  print ( '' )
  mu = 10.0
  sigma = 2.0
  x = rng.normal ( mu, sigma, 10000 )
  x_mean = np.mean ( x )
  x_std = np.std ( x )
  x_min = np.min ( x )
  x_max = np.max ( x )
  print ( '  x = rng.normal ( ', mu, ',', sigma, ', 1000 )' )
  print ( '  x[0:5] = ', x[0:5] )
  print ( '  mean ( x ) = ', x_mean, ' Distribution mean = ', mu )
  print ( '  std ( x ) = ', x_std, ' Distribution std = ', sigma )
  print ( '  min ( x ) = ', x_min, ' Distribution min = -inf' )
  print ( '  max ( x ) = ', x_max, ' Distribution max = +inf' )
#
#  Randomly permute a vector.
#
  print ( '' )
  print ( '  Randomly permute a vector.' )
  print ( '' )
  a = np.arange ( 10 )
  print ( '  a = np.arange ( 10 )' )
  print ( a )
  b = rng.permutation ( a )
  print ( '  b = np.permutation ( a )' )
  print ( b )
#
#  Randomly permute a matrix.
#
  print ( '' )
  print ( '  Randomly permute a matrix.' )
  print ( '' )
  c = np.array ( [ \
    [  0,  1,  2,  3,  4 ], \
    [ 10, 11, 12, 13, 14 ], \
    [ 20, 21, 22, 23, 24 ] ] )
  print ( '  c:' )
  print ( c )
  d = rng.permutation ( c )
  print ( '  d = rng.permutation ( c )' )
  print ( d )
  print ( '  e = np.transpose ( d )' )
  e = np.transpose ( d )
  print ( e )
  print ( '  f = rng.permutation ( e )' )
  f = rng.permutation ( e )
  print ( f )
  print ( '  g = rng.transpose ( e )' )
  g = np.transpose ( f )
  print ( g )

  defg = np.transpose ( rng.permutation ( np.transpose ( rng.permutation ( c ) ) ) )
  print ( '  defg = np.transpose ( rng.permutation ( np.transpose ( rng.permutation ( c ) ) ) )' )
  print ( defg )
#
#  Randomly select 3 items from 12, using permutation.
#
  print ( '' )
  print ( '  Randomly select 3 items from 12, using permutation.' )
  print ( '' )
  months = np.array ( [ 'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec' ] )
  i = np.arange ( 10 )
  j = rng.permutation ( i )
  three_months = months[j[0:3]]
  print ( three_months )
#
#  Randomly select 5 items from a list of 52.  (No repeats)
#
  print ( '' )
  print ( '  Randomly select 5 cards from 52, no repeats' )
  print ( '' )
  deck = np.arange ( 52 )
  print ( '  deck = np.arange ( 52 )' )
  print ( deck )
  shuffled_deck = rng.permutation ( deck )
  print ( '  shuffled_deck = rng.permutation ( deck )' )
  print ( shuffled_deck )
  hand = shuffled_deck[0:5]
  print ( '  hand = shuffled_deck[0:5]' )
  print ( hand )
#
#  Randomly select 5 cards from a deck of 52.  (Could have repeats)
#
  print ( '' )
  print ( '  Randomly select 5 cards from 52, can have repeats' )
  print ( '' )
  deck = np.arange ( 52 )
  print ( '  deck = np.arange ( 52 )' )
  print ( deck )
  b = rng.choice ( deck, size = 5 )
  print ( '  b = rng.choice ( deck, size = 5 )' )
  print ( b )
#
#  Generate a random win/loss/tie record.
#
  print ( '' )
  print ( '  Randomly select W/L/T values for a game season' )
  print ( '' )
  winlosstie = [ 'win', 'loss', 'tie' ]
  print ( "  winlosstie = [ 'win', 'loss', 'tie' ]" )
  print ( winlosstie )
  record = rng.choice ( winlosstie, size = 10 )
  print ( '  record = rng.choice ( winlosstie, size = 10 )' )
  print ( record )
  record = rng.choice ( winlosstie, size = 10 )
  print ( '  record = rng.choice ( winlosstie, size = 10 )' )
  print ( record )
#
#  Generate a binary matrix.
#
  print ( '' )
  print ( '  Generate a binary matrix.' )
  print ( '' )
  binary = [ 0, 1 ]
  print ( '  binary = [ 0, 1 ]' )
  print ( binary )
  binary_matrix = rng.choice ( binary, size = ( 4, 5 ) )
  print ( '  binary_matrix = rng.choice ( binary, size = ( 4, 5 ) )' )
  print ( binary_matrix )
#
#  Integral estimation.
#
  print ( '' )
  print ( '  Compute estimate Q for integral(a<=x<=b) f(x) dx' )
  print ( '' )
  print ( '  N    Q    Error' )
  print ( '' )
  a = 1.0
  b = 10.0
  exact = 1.0 / 2.0 * np.log ( ( 2.0 * b + 3.0 ) / ( 2.0 * a + 3.0 ) )
  for n in [ 10, 100, 1000 ]:
    x = a + ( b - a ) * rng.random ( n )
    fx = 1.0 / ( 2.0 * x + 3.0 )
    q = ( b - a ) * sum ( fx ) / n
    err = abs ( q - exact )
    print ( n, q, err )
#
#  Area estimation.
#
  print ( '' )
  print ( '  Compute estimate Q for area of ellipse' )
  print ( '' )
  print ( '  N    Q    Error' )
  print ( '' )
  exact = 2.0 * np.pi
  inside = 0
  xlo = -2.0
  xhi = 2.0
  ylo = -1.0
  yhi = 1.0
  for n in [ 10, 100, 1000 ]:
    x = xlo + ( xhi - xlo ) * rng.random ( n )
    y = ylo + ( yhi - ylo ) * rng.random ( n )
    fx = 0.25 * x**2 + y**2
    inside = sum ( fx <= 1.0 )
    q = ( xhi - xlo ) * ( yhi - ylo ) * inside / n
    err = abs ( q - exact )
    print ( n, q, err )
#
#  Expected value estimation.
#
  print ( '' )
  print ( '  Estimate average distance between two points in unit square.')
  print ( '' )
  print ( '  N    Estimate  Error' )
  print ( '' )
  exact = ( np.sqrt ( 2.0 ) + 2.0 + 5.0 * np.log ( 1.0 + np.sqrt ( 2.0 ) ) ) / 15.0
  for n in [ 10, 100, 1000, 10000 ]:
    x1 = rng.random ( [ n, 2 ] )
    x2 = rng.random ( [ n, 2 ] )
    dx = x1 - x2
    dist = np.sqrt ( dx[:,0]**2 + dx[:,1]**2 )
    dist_ave = np.mean ( dist )
    err = np.abs ( exact - dist_ave )
    print ( n, dist_ave, err )
#
#  L2 norm <= Frobenius norm?
#
  print ( '' )
  print ( '  ||A||2 <= ||A||F always?' )
  print ( '' )
  test_num = 1000
  true_num = 0
  for test in range ( 0, test_num ):
    A = rng.standard_normal ( ( 3, 3 ) )
    A_norm2 = np.linalg.norm ( A, 2 )
    A_normf = np.linalg.norm ( A, 'fro' )
    if ( A_norm2 <= A_normf ):
      true_num = true_num + 1

  print ( '  Statement was true ', true_num, ' times out of ', test_num )

  return

if ( __name__ == "__main__" ):
  python11 ( )