#! /usr/bin/env python
#
def genmul ( n, p, ncat ):

#*****************************************************************************80
#
## GENMUL generates a multinomial random deviate.
#
#  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:
#
#    Luc Devroye,
#    Non-Uniform Random Variate Generation,
#    Springer, 1986,
#    ISBN: 0387963057,
#    LC: QA274.D48.
#
#  Parameters:
#
#    Input, integer N, the number of events, which will be classified into 
#    one of the NCAT categories.
#
#    Input, real P(NCAT-1).  P(I) is the probability that an event
#    will be classified into category I.  Thus, each P(I) must be between
#    0.0 and 1.0.  Only the first NCAT-1 values of P must be defined since
#    P(NCAT) would be 1.0 minus the sum of the first NCAT-1 P's.
#
#    Input, integer NCAT, the number of categories.
#
#    Output, integer IX(NCAT), a random observation from
#    the multinomial distribution.  All IX(i) will be nonnegative and their
#    sum will be N.
#
  import numpy as np
  from ignbin import ignbin
  from sys import exit

  if ( n < 0 ):
    print ( '' )
    print ( 'GENMUL - Fatal error!' )
    print ( '  N < 0' )
    exit ( 'GENMUL - Fatal error!' )

  if ( ncat <= 1 ):
    print ( '' )
    print ( 'GENMUL - Fatal error!' )
    print ( '  NCAT <= 1' )
    exit ( 'GENMUL - Fatal error!' )

  for i in range ( 0, ncat - 1 ):

    if ( p[i] < 0.0 ):
      print ( ' ' )
      print ( 'GENMUL - Fatal error!' )
      print ( '  Some P(i) < 0.' )
      exit ( 'GENMUL - Fatal error!' )

    if ( 1.0 < p[i] ):
      print ( ' ' )
      print ( 'GENMUL - Fatal error!' )
      print ( '  Some 1 < P(i).' )
      exit ( 'GENMUL - Fatal error!' )

  ptot = np.sum ( p[0:ncat-1] )

  if ( 0.99999 < ptot ):
    print ( ' ' )
    print ( 'GENMUL - Fatal error!' )
    print ( '  1 < Sum of P().' )
    exit ( 'GENMUL - Fatal error!' )
#
#  Initialize variables.
#
  ntot = n
  ptot = 1.0
  ix = np.zeros ( ncat )
#
#  Generate the observation.
#
  for icat in range ( 0, ncat - 1 ):
    prob = p[icat] / ptot
    ix[icat] = ignbin ( ntot, prob )
    ntot = ntot - ix[icat]
    if ( ntot <= 0 ):
      return ix
    ptot = ptot - p[icat]

  ix[ncat-1] = ntot

  return ix

def genmul_test ( phrase ):

#*****************************************************************************80
#
## GENMUL_TEST tests GENMUL.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    03 September 2018
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from initialize import initialize
  from phrtsd import phrtsd
  from set_initial_seed import set_initial_seed

  print ( '' )
  print ( 'GENMUL_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  GENMUL generates a multinomial random deviate.' )
#
#  Initialize the generators.
#
  initialize ( )
#
#  Set the seeds based on the phrase.
#
  seed1, seed2 = phrtsd ( phrase )
#
#  Initialize all generators.
#
  set_initial_seed ( seed1, seed2 )

  n = 100
  p = np.array ( [ 0.35, 0.15, 0.05, 0.30 ] )
  ncat = 5

  print ( '' )
  print ( '  Try 10 successive trials:' )
  print ( '' )

  for test in range ( 0, 10 ):
    ix = genmul ( n, p, ncat )
    print ( '    ' ),
    for i in range ( 0, ncat ):
      print ( '%2d' % ( ix[i] ) ),
    print ( '' )
#
#  Terminate.
#
  print ( '' )
  print ( 'GENMUL_TEST:' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  phrase = 'Randomize'
  genmul_test ( phrase )
  timestamp ( )