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

#*****************************************************************************80
#
## random_test() tests random().
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    09 August 2022
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'random_test():' )
  print ( '  Python version: ' + platform.python_version ( ) )
  print ( '  Test random().' )

  random_histogram ( )
  normal_histogram ( )
  dice_test ( )
  logistic_sample_test ( )
  normal_plot_test ( )
  numpy_random_random_test ( )
  spd_test ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'random_test():' )
  print ( '  Normal end of execution.' )

  return

def random_histogram ( ):

#*****************************************************************************80
#
## random_histogram() makes a histogram of n samples of uniform random values.
#
  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'random_histogram() makes a histogram of n samples' )
  print ( '  of uniform random values from numpy.random.random().' )

  n = 1000
  plt.clf ( )
  r = np.random.random ( n )
  plt.hist ( r, bins = 50 )
  plt.grid ( True )
  filename = 'random_histogram.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( )
  plt.close ( )

  mean_rand = np.mean ( r )
  var_rand = np.var ( r )
  print ( '' )
  print ( '  mean(r) = ', mean_rand )
  print ( '  var(r)  = ', var_rand )

  return

def normal_histogram ( ):

#*****************************************************************************80
#
## normal_histogram() makes a histogram of n samples of normal random values.
#
  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'normal_histogram() makes a histogram of n samples' )
  print ( '  of normal random values from numpy.random.normal().' )

  n = 1000
  plt.clf ( )
  r = np.random.normal ( loc = 2.0, scale = 1.5, size = n )
  plt.hist ( r, bins = 50 )
  plt.grid ( True )
  filename = 'normal_histogram.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( )
  plt.close ( )

  normal_mean = np.mean ( r )
  normal_stdev = np.std ( r )
  print ( '  normal mean  = ', normal_mean )
  print ( '  normal stdev = ', normal_stdev )

  return

def dice_test ( ):

#*****************************************************************************80
#
## dice_test() simulates the Newton-Pepys dice puzzle.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    29 September 2022
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  print ( '' )
  print ( 'dice_test():' )
  print ( '  Simulate the Newton-Pepys dice puzzle.' )
  test_num = 1000
  for total in range ( 1, 4 ):
    freq = 0
    for test in range ( 0, test_num ):
      dice = np.random.randint ( low = 1, high = 7, size = 6 * total )
      sixes = np.where ( dice == 6 )
      if ( total <= np.size ( sixes ) ):
        freq = freq + 1
    print ( '  Case', total, ' probability estimate is ', freq / test_num )

  return

def dif2_matrix ( m, n ):

#*****************************************************************************80
#
## dif2_matrix() returns the second difference matrix.
#
#  Example:
#
#    N = 5
#
#    2 -1  .  .  .
#   -1  2 -1  .  .
#    . -1  2 -1  .
#    .  . -1  2 -1
#    .  .  . -1  2
#
#  Rectangular Properties:
#
#    A is banded, with bandwidth 3.
#
#    A is tridiagonal.
#
#    Because A is tridiagonal, it has property A (bipartite).
#
#    A is integral: int ( A ) = A.
#
#    A is Toeplitz: constant along diagonals.
#
#  Square Properties:
#
#    A is symmetric: A' = A.
#
#    Because A is symmetric, it is normal.
#
#    Because A is normal, it is diagonalizable.
#
#    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
#
#    A is positive definite.
#
#    A is an M matrix.
#
#    A is weakly diagonally dominant, but not strictly diagonally dominant.
#
#    A has an LU factorization A = L * U, without pivoting.
#
#      The matrix L is lower bidiagonal with subdiagonal elements:
#
#        L(I+1,I) = -I/(I+1)
#
#      The matrix U is upper bidiagonal, with diagonal elements
#
#        U(I,I) = (I+1)/I
#
#      and superdiagonal elements which are all -1.
#
#    A has a Cholesky factorization A = L * L', with L lower bidiagonal.
#
#      L(I,I) =    sqrt ( (I+1) / I )
#      L(I,I-1) = -sqrt ( (I-1) / I )
#
#    The eigenvalues are
#
#      LAMBDA(I) = 2 + 2 * COS(I*PI/(N+1))
#                = 4 SIN^2(I*PI/(2*N+2))
#
#    The corresponding eigenvector X(I) has entries
#
#       X(I)(J) = sqrt(2/(N+1)) * sin ( I*J*PI/(N+1) ).
#
#    Simple linear systems:
#
#      x = (1,1,1,...,1,1),   A*x=(1,0,0,...,0,1)
#
#      x = (1,2,3,...,n-1,n), A*x=(0,0,0,...,0,n+1)
#
#    det ( A ) = N + 1.
#
#    The value of the determinant can be seen by induction,
#    and expanding the determinant across the first row:
#
#      det ( A(N) ) = 2 * det ( A(N-1) ) - (-1) * (-1) * det ( A(N-2) )
#                = 2 * N - (N-1)
#                = N + 1
#
#    The family of matrices is nested as a function of N.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    21 January 2015
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Robert Gregory, David Karney,
#    Example 3.18,
#    A Collection of Matrices for Testing Computational Algorithms,
#    Wiley, New York, 1969, page 45, 
#    LC: QA263.G68.
#
#    Morris Newman, John Todd,
#    Example A8,
#    The evaluation of matrix inversion programs,
#    Journal of the Society for Industrial and Applied Mathematics,
#    Volume 6, Number 4, pages 466-476, 1958.
#
#    John Todd,
#    Example A8,
#    Basic Numerical Mathematics,
#    Volume 2: Numerical Algebra,
#    Academic Press, 1977, page 1.
#
#    Joan Westlake,
#    Test Matrix A15,
#    A Handbook of Numerical Matrix Inversion and Solution of Linear Equations,
#    John Wiley, 1968.
#
#  Input:
#
#    integer M, N, the number of rows and columns of A.
#
#  Output:
#
#    real A(M,N), the matrix.
#
  import numpy as np

  a = np.zeros ( [ m, n ] )

  for i in range ( 0, m ):
    if ( 0 <= i - 1 and i - 1 < n ):
      a[i,i-1] = -1.0
    if ( i < n ):
      a[i,i] = 2.0
    if ( i + 1 < n ):
      a[i,i+1] = -1.0

  return a

def fibonacci1_matrix ( n, f1, f2 ):

#*****************************************************************************80
#
## fibonacci1_matrix() returns the FIBONACCI1 matrix.
#
#  Example:
#
#    N = 5
#    F1 = 1, F2 = 2
#
#    1  2  3  5  8
#    2  3  5  8 13
#    3  5  8 13 21
#    5  8 13 21 34
#    8 13 21 34 55
#
#  Properties:
#
#    A is symmetric: A' = A.
#
#    If F1 and F2 are integral, then so is A.
#
#    If A is integral, then det ( A ) is integral, and
#    det ( A ) * inverse ( A ) is integral.
#
#    A is a Hankel matrix: constant along anti-diagonals.
#
#    If B is the Fibonacci iteration matrix:
#
#      B * A(F1,F2) = A(F2,F2+F1) = A(F2,F3)
#
#    and in general,
#
#      B^N * A(F1,F2) = A(F(N+1),F(N+2))
#
#    For 2 < N, the matrix is singular, because row 3 is the sum
#    of row 1 and row 2.
#
#    For 2 <= N,
#      rank ( A ) = 2
#
#    If N = 1, then
#      det ( A ) = 1
#    else if N = 2 then
#      det ( A ) = -1
#    else if 1 < N then
#      det ( A ) = 0
#
#    The family of matrices is nested as a function of N.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    27 January 2015
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the order of the matrix.
#
#    real F1, F2, the first two elements of the sequence
#    that will generate the Fibonacci sequence.
#
#  Output:
#
#    real A(N,N), the matrix.
#
  import numpy as np

  a = np.zeros ( [ n, n ] )

  a[0,0] = f1

  if ( 1 < n ):

    a[1,0] = f2
    a[0,1] = f2

    fnm2 = f1
    fnm1 = f2
    fn   = fnm1 + fnm2

    for k in range ( 2, n + n - 1 ):

      i = min ( k,     n - 1 )
      j = max ( 0, k - n + 1 )

      while ( 0 <= i and j < n ):
        a[i,j] = fn
        i = i - 1
        j = j + 1

      fnm2 = fnm1
      fnm1 = fn
      fn   = fnm1 + fnm2

  return a

def logistic_sample_test ( ):

#*****************************************************************************80
#
## logistic_sample_test() demonstrates sampling on logistic distribution.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    11 August 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'logistic_sample_test():' )
  print ( '  Use the inverse CDF of the logistic distribution' )
  print ( '  to produce 1000 sample points.' )

  n = 10000
  c = np.random.random ( n )

  a = 0.0
  b = 1.0
  x = logistic_cdf_inv ( c, a, b )

  plt.clf ( )
  plt.hist ( x, bins = 25, density = True )
  x2 = np.linspace ( -8.0, 8.0, 101 )
  pdf = logistic_pdf ( x2, a, b )
  plt.plot ( x2, pdf, 'r-', linewidth = 3 )
  plt.grid ( True )
  plt.xlabel ( '<-- X -->' )
  plt.ylabel ( '<-- Probability -->' )
  plt.title ( 'Estimated and exact logistic PDF' )
  filename = 'logistic_sample.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.close ( )

  return

def logistic_cdf ( x, a, b ):

#*****************************************************************************80
#
## logistic_cdf() evaluates the Logistic CDF.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    04 April 2016
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real X, the argument of the CDF.
#
#    real A, B, the parameters of the PDF.
#    0.0 < B.
#
#  Output:
#
#    real CDF, the value of the CDF.
#
  import numpy as np

  cdf = 1.0 / ( 1.0 + np.exp ( ( a - x ) / b ) )

  return cdf

def logistic_cdf_inv ( cdf, a, b ):

#*****************************************************************************80
#
## logistic_cdf_inv() inverts the Logistic CDF.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    04 April 2016
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real CDF, the value of the CDF.
#    0.0 <= CDF <= 1.0.
#
#    real A, B, the parameters of the PDF.
#    0.0 < B.
#
#  Output:
#
#    real X, the corresponding argument.
#
  import numpy as np

  x = a - b * np.log ( ( 1.0 - cdf ) / cdf )

  return x

def logistic_pdf ( x, a, b ):

#*****************************************************************************80
#
## logistic_pdf() evaluates the Logistic PDF.
#
#  Discussion:
#
#    PDF(X)(A,B) = EXP ( ( A - X ) / B ) /
#      ( B * ( 1 + EXP ( ( A - X ) / B ) )^2 )
#
#    The Logistic PDF is also known as the Sech-Squared PDF.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    04 April 2016
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real X, the argument of the PDF.
#
#    real A, B, the parameters of the PDF.
#    0.0 < B.
#
#  Output:
#
#    real PDF, the value of the PDF.
#
  import numpy as np

  pdf = np.exp ( ( a - x ) / b ) / ( b * ( 1.0 + np.exp ( ( a - x ) / b ) ) ** 2 )

  return pdf

def normal_plot_test ( ):

#*****************************************************************************80
#
## normal_plot_test() plots the normal distribution.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    26 September 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'normal_plot_test():' )
  print ( '  Plot the normal PDF for several parameter sets.' )

  x = np.linspace ( -4.0, 4.0, 101 )
  y1 = normal_pdf ( x, 0.0, 1.0 )
  y2 = normal_pdf ( x, 1.0, 1.0 )
  y3 = normal_pdf ( x, 0.0, 2.0 )
  z = np.zeros ( 101 )

  plt.clf ( )
  plt.plot ( x, y1, 'r-', linewidth = 3 )
  plt.plot ( x, y2, 'g-', linewidth = 3 )
  plt.plot ( x, y3, 'b-', linewidth = 3 )
  plt.plot ( x, z, 'k--', linewidth = 2 )
  plt.grid ( True )
  plt.xlabel ( '<-- X -->' )
  plt.ylabel ( '<-- Probability -->' )
  plt.title ( 'Three normal distributions' )
  plt.legend ( ['loc=0,scale=1', 'loc=1,scale=1', 'loc=0,scale=2'] )
  filename = 'normal_plot.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.close ( )

  return

def normal_pdf ( x, mu, sigma ):

#*****************************************************************************80
#
## normal_pdf() evaluates the Normal PDF.
#
#  Discussion:
#
#    The normal PDF is also known as the Gaussian PDF.
#
#  Formula:
#
#    PDF(X;MU,SIGMA) = 
#      EXP ( - 0.5 * ( ( X - MU ) / SIGMA )^2 ) / SQRT ( 2 * PI * SIGMA^2 )
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    21 March 2016
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real X(), the argument of the PDF.
#
#    real MU, SIGMA, the mean and standard deviation.
#    SIGMA should be nonzero.
#
#  Output:
#
#    real PDF(), the value of the PDF.
#
  import numpy as np

  pdf = np.exp ( - 0.5 * ( ( x - mu ) / sigma ) ** 2 ) \
    / np.sqrt ( 2.0 * np.pi * sigma ** 2 )

  return pdf

def numpy_random_random_test ( ):

#*****************************************************************************80
#
## logistic_pdf() evaluates the Logistic PDF.
#
#  Modified:
#
#    23 September 2022
#
  import numpy as np

  print ( '' )
  print ( 'numpy_random_random_test():' )
  print ( '  Demonstrate some features of the numpy random number generators.' )

  x = np.random.random ( )
  print ( '  x = np.random.random ( ): ' )
  print ( x )

  s = np.random.random ( 1 )
  print ( '  s = np.random.random ( size = 1 ): ' )
  print ( s )

  v = np.random.random ( 2 )
  print ( '  s = np.random.random ( size = 2 ): ' )
  print ( s )

  a = np.random.random ( [ 3, 4 ] )
  print ( '  a = np.random.random ( [ 3, 4 ] ): ' )
  print ( a )

  print ( '' )
  print ( '  np.random.seed() can reset the seed:' )
  print ( '' )

  seed_value = 123456789
  np.random.seed ( seed_value )
  print ( '  np.random.seed ( seed_value ) ' )
  x = np.random.random ( )
  print ( '  x = np.random.random ( ): ' )
  print ( x )
  x = np.random.random ( )
  print ( '  x = np.random.random ( ): ' )
  print ( x )
  np.random.seed ( seed_value )
  print ( '  np.random.seed ( seed_value ) ' )
  x = np.random.random ( )
  print ( '  x = np.random.random ( ): ' )
  print ( x )

  print ( '' )
  print ( '  x =  np.random.uniform(low=a,high=b,size=?) returns values in [a,b]' )
  print ( '' )
  a = 10
  b = 15
  print ( '  Interval is [', a, ',', b, ']' )
  x = np.random.uniform ( low = a, high = b )
  print ( '  x = np.random.uniform ( low = a, high = b ): ' )
  print ( x )
  v = np.random.uniform (  low = a, high = b, size = 5 )
  print ( '  v = np.random.uniform (  low = a, high = b, size = 5 ): ' )
  print ( v )

  print ( '' )
  print ( '  x = np.random.standard_normal(size=?)' )
  print ( '  returns normal random values with mean 0, standard deviation 1.' )
  print ( '' )
  v = np.random.standard_normal ( size = [ 3, 4 ] )
  print ( '  v = np.random.standard_normal ( size = [ 3, 4 ] ): ' )
  print ( v )

  print ( '' )
  print ( '  x = np.random.normal(loc=mu,scale=sigma,size=?)' )
  print ( '  returns normal random values with mean mu, standard deviation sigma.' )
  print ( '' )
  v = np.random.normal ( loc = 10, scale = 0.5, size = [ 3, 4 ] )
  print ( '  v = np.random.normal ( loc = 10, scale = 0.5, size = [ 3, 4 ] ): ' )
  print ( v )

  return

def spd_test ( ):

#*****************************************************************************80
#
## spd_test() uses random sampling to test if a matrix is SPD.
#
#  Discussion:
#
#    Check that x'Ax > 0 for all x.  Test is bad if we only use "random()"
#    because all entries are nonnegative.  Test is better if we use
#    standard_normal().
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    23 September 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'spd_test():' )
  print ( '  Matrix A is SPD if x\'Ax > 0 for all nonzero vectors x' )

  m = 5
  n = 5
  A1 = dif2_matrix ( m, n )
  A2 = fibonacci1_matrix ( n, 1.0, 2.0 )

  for A, name in [ [ A1, 'dif2' ], [ A2, 'fibonacci' ] ]:

    fail = 0
    for test in range ( 0, 100 ):
      x = np.random.random ( size = n )
      xtAx = np.dot ( x.T, np.matmul ( A, x ) )
      if ( xtAx <= 0.0 ):
        fail = fail + 1
    if ( 0 == fail ):
      print ( '  ' + name + ' matrix passes SPD test #1' )
    else:
      print ( '  ' + name + ' matrix fails SPD test #1' )

    fail = 0
    for test in range ( 0, 100 ):
      x = np.random.standard_normal ( size = n )
      xtAx = np.dot ( x.T, np.matmul ( A, x ) )
      if ( xtAx <= 0.0 ):
        fail = fail + 1
    if ( 0 == fail ):
      print ( '  ' + name + ' matrix passes SPD test #2' )
    else:
      print ( '  ' + name + ' matrix fails SPD test #2' )

  return

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    21 August 2019
#
#  Author:
#
#    John Burkardt
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return

if ( __name__ == '__main__' ):
  timestamp ( )
  random_test ( )
  timestamp ( )