#! /usr/bin/env python3
#
def wolf_pack_expected ( pack_size ):

#*****************************************************************************80
#
## wolf_pack_expected(): expected value of the general full deck problem.
#
#  Discussion:
#
#    The expected value, given a pack of N wolves, is n*Hn, where
#    Hn is the n-th harmonic number.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    11 February 2023
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Herbert Wilf,
#    Some New Aspects of the Coupon Collector's Problem,
#    SIAM Review,
#    Volume 48, Number 3, September 2006, pages 549-565.
#
#  Input:
#
#    integer pack_size: the number of wolves in the pack.
#
#  Output:
#
#    real ev: the expected number of wolves that must be collected 
#    and releasted in order to see all of them at least once.
#
  ev = 0.0
  for i in range ( 1, pack_size + 1 ):
    ev = ev + pack_size / i

  return ev

def wolf_pack_expected_test ( ):

#*****************************************************************************80
#
## wolf_pack_expected_test() tests wolf_pack_expected().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    06 July 2022
#
#  Author:
#
#    John Burkardt
#
  print ( '' )
  print ( 'wolf_pack_expected_test():' )
  print ( '  wolf_pack_expected() reports the expected value for' )
  print ( '  the number of wolves trapped and released before seeing every one.' )
  print ( '' )
  print ( '  pack_size  Expected' )
  print ( '' )
  for pack_size in [ 100, 200, 300, 400 ]:
    ev = wolf_pack_expected ( pack_size )
    print ( '  %3d    %10.2f' % ( pack_size, ev ) )

  return

def wolf_pack_histogram ( n, pack_size ):

#*****************************************************************************80
#
## wolf_pack_histogram() histograms wolf_pack statistics.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    11 July 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the number of trials to make.
#
  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'wolf_pack_histogram():' )
  print ( '  Create a histogram of the wolf pack process.' )

  m = np.zeros ( n )
  for i in range ( 0, n ):
    deck = wolf_pack_try ( pack_size )
    m[i] = np.sum ( deck )

  plt.hist ( m, bins = 20, rwidth = 0.95, density = True )
  plt.grid ( True )
  plt.xlabel ( '<-- Number of draws -->' )
  plt.ylabel ( '<-- Probability -->' )
  plt.title ( 'PDF for number of draws to get full pack' )
  filename = 'wolf_pack_pdf.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.show ( block = False )
  plt.close ( )
 
  return

def wolf_pack_simulation_test ( ):

#*****************************************************************************80
#
## wolf_pack_simulation_test() tests wolf_pack_simulation().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    11 February 2023
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'wolf_pack_simulation_test():' )
  print ( '  Python version: ' + platform.python_version ( ) )
  print ( '  wolf_pack_simulation() simulates a process in which' )
  print ( '  a wolf is selected  from a pack, then returned.' )
  print ( '  The process is continued until all wolves have been seen' )
  print ( '  at least once.' )

  wolf_pack_try_test ( )
  wolf_pack_stats_test ( )
  wolf_pack_expected_test ( )
  n = 2000
  pack_size = 300
  wolf_pack_histogram ( n, pack_size )
  wolf_pack_variance_test ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'wolf_pack_simulation_test():' )
  print ( '  Normal end of execution.' )

  return

def wolf_pack_stats ( n, pack_size ):

#*****************************************************************************80
#
## wolf_pack_stats() returns statistics for N tries of the wolf pack problem.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    11 February 2023
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer n: the number of times the process will be carried out.
#
#  Output:
#
#    integer tries_min, real tries_mean, integer tries_max, real tries_std: 
#    the minimum, average, maximum and std for the number of wolves drawn 
#    over all the processes.
#
  import numpy as np

  m = np.zeros ( n )
  for i in range ( 0, n ):
    deck = wolf_pack_try ( pack_size )
    m[i] = np.sum ( deck )

  tries_min = np.min ( m )
  tries_max = np.max ( m )
  tries_mean = np.mean ( m )
  tries_std = np.std ( m )

  return tries_min, tries_mean, tries_max, tries_std

def wolf_pack_stats_test ( ):

#*****************************************************************************80
#
## wolf_pack_stats_test() tests wolf_pack_stats().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    06 July 2022
#
#  Author:
#
#    John Burkardt
#
  print ( '' )
  print ( 'wolf_pack_stats_test():' )
  print ( '  wolf_pack_stats() reports statistics for N instances' )
  print ( '  of the wolf pack process.' )
  print ( '' )
  print ( '       N  Min    Mean  Max    Std' )
  print ( '' )

  pack_size = 300
  for n in [ 10, 20, 50, 100, 1000, 10000 ]:
    tmin, tmean, tmax, tstd = wolf_pack_stats ( n, pack_size )
    print ( '   %5d  %3d  %6.2f  %2d  %6.2f' % ( n, tmin, tmean, tmax, tstd ) )

  return

def wolf_pack_try ( pack_size ):

#*****************************************************************************80
#
## wolf_pack_try() returns all the cards drawn for a full deck problem.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    11 February 2023
#
#  Author:
#
#    John Burkardt
#
#  Output:
#
#    integer count(pack_size): the number of times each wolf was drawn.
#
  from numpy.random import default_rng
  import numpy as np

  rng = default_rng()

  count = np.zeros ( pack_size )
  while ( True ):
    card = rng.integers ( low = 0, high = pack_size )
    count[card] = count[card] + 1
    if ( np.all ( 0 < count ) ):
      break

  return count

def wolf_pack_try_test ( ):

#*****************************************************************************80
#
## wolf_pack_try_test() tests wolf_pack_try().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    11 February 2023
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  print ( '' )
  print ( 'wolf_pack_try_test():' )
  print ( '  wolf_pack_try() repeatedly draws a random wolf' )
  print ( '  from a pack until all cards have been seen' )
  print ( '  at least once.' )
  print ( '' )
  print ( '  Try  Min    Mean Max  Total' )
  print ( '' )

  pack_size = 300
  for i in range ( 0, 10 ):
    cards = wolf_pack_try ( pack_size )
    print ( '   %2d  %3d  %6.2f  %2d   %4d' % \
      ( i, np.min ( cards ), np.mean ( cards ), np.max ( cards ), np.sum ( cards ) ) )

  return

def wolf_pack_variance ( pack_size ):

#*****************************************************************************80
#
## wolf_pack_variance() returns the variance of the wolf pack problem.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    06 July 2022
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Herbert Wilf,
#    Some New Aspects of the Coupon Collector's Problem,
#    SIAM Review,
#    Volume 48, Number 3, September 2006, pages 549-565.
#
#  Input:
#
#    integer pack_size: the number of cards in the deck.
#
#  Output:
#
#    real variance: the variance in the number of wolves that must be 
#    collected in order to acquire at least one of every kind.
#
  variance = 0.0
  for i in range ( 2, pack_size + 1 ):
    variance = variance + ( i - 1 ) / ( pack_size - i + 1 )**2

  variance = pack_size * variance

  return variance

def wolf_pack_variance_test ( ):

#*****************************************************************************80
#
## wolf_pack_variance_test() tests wolf_pack_variance().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    06 July 2022
#
#  Author:
#
#    John Burkardt
#
  print ( '' )
  print ( 'wolf_pack_variance_test():' )
  print ( '  wolf_pack_variance() reports the variance for' )
  print ( '  the number of wolves drawn in the full pack process.' )
  print ( '' )
  print ( '    N     Variance' )
  print ( '' )
  for pack_size in [ 100, 200, 300, 400 ]:
    variance = wolf_pack_variance ( pack_size )
    print ( '  %3d  %10.2f' % ( pack_size, variance ) )

  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 ( )
  wolf_pack_simulation_test ( )
  timestamp ( )