subroutine amino_print ( acid_num, acid_sym )

!*****************************************************************************80
!
!! AMINO_PRINT prints the amino acid parameters.
!
!  Modified:
!
!    23 November 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ACID_NUM, the number of amino acids.
!
!    Input, character ACID_SYM(ACID_NUM), the one letter amino acid codes.
!
  implicit none

  integer acid_num

  integer acid_i
  character ( len = 27 ) acid_name
  character acid_sym(acid_num)
  character c

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  I  Amino Acid Symbol'
  write ( *, '(a)' ) ' '
  do acid_i = 1, acid_num
    c = acid_sym(acid_i)
    call ch_to_amino_name ( c, acid_name )
    write ( *, '(i3,2x,a,2x,a)' ) acid_i, acid_sym(acid_i), acid_name
  end do

  return
end
subroutine binomial_sample ( a, b, x )

!*****************************************************************************80
!
!! BINOMIAL_SAMPLE samples the Binomial PDF.
!
!  Modified:
!
!    02 February 1999
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    William Kennedy, James Gentle,
!    Algorithm BU,
!    Statistical Computing,
!    Dekker, 1980.
!
!  Parameters:
!
!    Input, integer A, the number of trials.
!    1 <= A.
!
!    Input, real B, the probability of success on one trial.
!    0.0 <= B <= 1.0.
!
!    Output, integer X, a sample of the PDF.
!
  implicit none

  integer a
  real b
  integer i
  real u
  integer x

  x = 0

  do i = 1, a

    call r4_random ( 0.0E+00, 1.0E+00, u )

    if ( u <= b ) then
      x = x + 1
    end if

  end do

  return
end
subroutine ch_cap ( c )

!*****************************************************************************80
!
!! CH_CAP capitalizes a single character.
!
!  Modified:
!
!    19 July 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, character C, the character to capitalize.
!
  implicit none

  character c
  integer itemp

  itemp = ichar ( c )

  if ( 97 <= itemp .and. itemp <= 122 ) then
    c = char ( itemp - 32 )
  end if

  return
end
function ch_eqi ( c1, c2 )

!*****************************************************************************80
!
!! CH_EQI is a case insensitive comparison of two characters for equality.
!
!  Example:
!
!    CH_EQI ( 'A', 'a' ) is .TRUE.
!
!  Modified:
!
!    14 August 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character C1, C2, the characters to compare.
!
!    Output, logical CH_EQI, the result of the comparison.
!
  implicit none

  logical ch_eqi
  character c1
  character c2
  character cc1
  character cc2

  cc1 = c1
  cc2 = c2

  call ch_cap ( cc1 )
  call ch_cap ( cc2 )

  if ( cc1 == cc2 ) then
    ch_eqi = .true.
  else
    ch_eqi = .false.
  end if

  return
end
subroutine ch_next ( line, cval, done )

!*****************************************************************************80
!
!! CH_NEXT "reads" space-separated characters from a string, one at a time.
!
!  Example:
!
!    Input:
!
!      LINE = ' A  B, C    DE  F'
!
!    Output:
!
!      'A', 'B', 'C', 'D', 'E', 'F', and then blanks.
!
!  Modified:
!
!    18 November 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) LINE, a string, presumably containing
!    characters, possibly separated by spaces or commas.
!
!    Output, character CVAL.  If DONE is FALSE, then CVAL contains the
!    "next" character read from LINE.  If DONE is TRUE, then
!    CVAL is blank.
!
!    Input/output, logical DONE.
!    On input with a fresh value of LINE, the user should set
!    DONE to TRUE.
!    On output, the routine sets DONE to FALSE if another character
!    was read, or TRUE if no more characters could be read.
!
  implicit none

  character cval
  logical done
  integer i
  character ( len = * ) line
  integer, save :: next = 1

  if ( done ) then
    next = 1
    done = .false.
  end if

  do i = next, len(line)

    if ( line(i:i) /= ' ' .and. line(i:i) /= ',' ) then
      cval = line(i:i)
      next = i + 1
      return
    end if

  end do

  done = .true.
  next = 1
  cval = ' '

  return
end
subroutine ch_to_amino_name ( c, amino_name )

!*****************************************************************************80
!
!! CH_TO_AMINO_NAME converts a character to an amino acid name.
!
!  Modified:
!
!    16 June 2000
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Carl Branden, John Tooze,
!    Introduction to Protein Structure,
!    Garland Publishing, 1991.
!
!  Parameters:
!
!    Input, character C, the one letter code for an amino acid.
!    Lower and upper case letters are treated the same.
!
!    Output, character ( len = * ) AMINO_NAME, the full name of the
!    corresponding amino acid.  The longest name is 27 characters.  
!    If the input code is not recognized, then AMINO_NAME will be 
!    set to '???'.
!
  implicit none

  integer, parameter :: n = 23

  character ( len = * ) amino_name
  character ( len = 27 ), dimension ( n ) :: amino_table = (/ &
    'Alanine                    ', &
    'Aspartic acid or Asparagine', &
    'Cysteine                   ', &
    'Aspartic acid              ', &
    'Glutamic acid              ', &
    'Phenylalanine              ', &
    'Glycine                    ', &
    'Histidine                  ', &
    'Isoleucine                 ', &
    'Lysine                     ', &
    'Leucine                    ', &
    'Methionine                 ', &
    'Asparagine                 ', &
    'Proline                    ', &
    'Glutamine                  ', &
    'Arginine                   ', &
    'Serine                     ', &
    'Threonine                  ', &
    'Valine                     ', &
    'Tryptophan                 ', &
    'Undetermined amino acid    ', &
    'Tyrosine                   ', &
    'Glutamic acid or Glutamine ' /)
  character c
  logical ch_eqi
  character, dimension ( n ) :: c_table = (/ &
    'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', &
    'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', &
    'X', 'Y', 'Z' /)
  integer i

  do i = 1, n
    if ( ch_eqi ( c, c_table(i) ) ) then
      amino_name = amino_table(i)
      return
    end if
  end do

  amino_name = '???'

  return
end
subroutine ch_to_digit ( c, digit )

!*****************************************************************************80
!
!! CH_TO_DIGIT returns the integer value of a base 10 digit.
!
!  Example:
!
!     C   DIGIT
!    ---  -----
!    '0'    0
!    '1'    1
!    ...  ...
!    '9'    9
!    ' '    0
!    'X'   -1
!
!  Modified:
!
!    04 August 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character C, the decimal digit, '0' through '9' or blank
!    are legal.
!
!    Output, integer DIGIT, the corresponding integer value.  If C was
!    'illegal', then DIGIT is -1.
!
  implicit none

  character c
  integer digit

  if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then

    digit = ichar ( c ) - 48

  else if ( c == ' ' ) then

    digit = 0

  else

    digit = -1

  end if

  return
end
subroutine comp_param_print ( acid_num, acid_sym, comp_max, comp_num, beta, &
  beta_sum, comp_weight )

!*****************************************************************************80
!
!! COMP_PARAM_PRINT prints the parameters for the mixture components.
!
!  Modified:
!
!    24 January 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ACID_NUM, the number of amino acids.
!
!    Input, character ACID_SYM(ACID_NUM), the one letter amino acid codes.
!
!    Input, integer COMP_MAX, the maximum number of Dirichlet mixture
!    components.
!
!    Input, integer COMP_NUM, the number of components in the Dirichlet mixture.
!
!    Input, real BETA(ACID_NUM,COMP_MAX); BETA(I,J) is the parameter
!    for the J-th acid in the I-th Dirichlet mixture component.
!
!    Input, real BETA_SUM(COMP_MAX), the sum of the values of
!    BETA(ACID_I,COMP_I) for a given component COMP_I.
!
!    Input, real COMP_WEIGHT(COMP_NUM), the mixture weight of each component.
!    These values should be nonnegative, and sum to 1.  They represent the
!    relative proportion of each component in the mixture.
!
  implicit none

  integer acid_num
  integer comp_max

  integer acid_i
  character acid_sym(acid_num)
  integer comp_i
  real beta(acid_num,comp_max)
  real beta_sum(comp_max)
  integer comp_num
  real comp_weight(comp_max)

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Number of components = ', comp_num
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' '
  write ( *, '(''Compon:'',20i8)' ) ( comp_i, comp_i = 1, comp_num )
  write ( *, '(''Weight:'',20f8.4)' ) comp_weight(1:comp_num)
  write ( *, '(a)' ) ' '

  do acid_i = 1, acid_num
    write ( *, '(i2,2x,a1,2x,20f8.4)' ) acid_i, acid_sym(acid_i), &
      beta(acid_i,1:comp_num)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a3,4x,20f8.4)' ) 'Sum', beta_sum(1:comp_num)

  return
end
subroutine dirichlet_mean ( n, a, mean )

!*****************************************************************************80
!
!! DIRICHLET_MEAN returns the means of the Dirichlet PDF.
!
!  Modified:
!
!    23 November 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of components.
!
!    Input, real A(N), the probabilities for each component.
!    Each A(I) should be nonnegative, and at least one should be positive.
!
!    Output, real MEAN(N), the means of the PDF.
!
  implicit none

  integer n

  real a(n)
  real mean(n)

  if ( any ( a(1:n) < 0.0 ) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DIRICHLET_MEAN - Fatal error!'
    write ( *, '(a)' ) '  At least one entry of A is negative!'
    stop
  end if

  if ( all ( a(1:n) == 0.0 ) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DIRICHLET_MEAN - Fatal error!'
    write ( *, '(a)' ) '  All entries of A are zero!'
    stop
  end if

  mean(1:n) = a(1:n)

  call r4vec_unit_sum ( n, mean )

  return
end
subroutine dirichlet_mix_check ( comp_num, elem_max, elem_num, a, comp_weight )

!*****************************************************************************80
!
!! DIRICHLET_MIX_CHECK checks the parameters of a Dirichlet mixture PDF.
!
!  Modified:
!
!    13 December 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer COMP_NUM, the number of components in the Dirichlet
!    mixture density, that is, the number of distinct Dirichlet PDF's
!    that are mixed together.
!
!    Input, integer ELEM_MAX, the leading dimension of A, which must
!    be at least ELEM_NUM.
!
!    Input, integer ELEM_NUM, the number of elements of an observation.
!
!    Input, real A(ELEM_MAX,COMP_NUM), the probabilities for element ELEM_NUM
!    in component COMP_NUM.
!    Each A(I,J) should be greater than or equal to 0.0.
!
!    Input, integer COMP_WEIGHT(COMP_NUM), the mixture weights of the densities.
!    These do not need to be normalized.  The weight of a given component is
!    the relative probability that that component will be used to generate
!    the sample.
!
  implicit none

  integer comp_num
  integer elem_max
  integer elem_num

  real a(elem_max,comp_num)
  integer comp_i
  real comp_weight(comp_num)
  integer elem_i
  logical positive

  do comp_i = 1, comp_num

    do elem_i = 1, elem_num
      if ( a(elem_i,comp_i) < 0.0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'DIRICHLET_MIX_CHECK - Fatal error!'
        write ( *, '(a)' ) '  A(ELEM,COMP) < 0.'
        write ( *, '(a,i6)' ) '  COMP = ', comp_i
        write ( *, '(a,i6)' ) '  ELEM = ', elem_i
        write ( *, '(a,g14.6)' ) '  A(COMP,ELEM) = ', a(elem_i,comp_i)
        stop
      end if
    end do

  end do

  positive = .false.

  do comp_i = 1, comp_num

    if ( comp_weight(comp_i) < 0.0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DIRICHLET_MIX_CHECK - Fatal error!'
      write ( *, '(a)' ) '  COMP_WEIGHT(COMP) < 0.'
      write ( *, '(a,i6)' ) '  COMP = ', comp_i
      write ( *, '(a,g14.6)' ) '  COMP_WEIGHT(COMP) = ', comp_weight(comp_i)
      stop
    else if ( comp_weight(comp_i) > 0.0 ) then
      positive = .true.
    end if

  end do

  if ( .not. positive ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DIRICHLET_MIX_CHECK - Fatal error!'
    write ( *, '(a)' ) '  All component weights are zero.'
    stop
  end if

  return
end
subroutine dirichlet_multinomial_pdf ( x, a, b, c, pdf )

!*****************************************************************************80
!
!! DIRICHLET_MULTINOMIAL_PDF evaluates a Dirichlet Multinomial PDF.
!
!  Formula:
!
!    PDF(X)(A,B,C) = Comb(A,B,X) * ( Gamma(C_Sum) / Gamma(C_Sum+A) )
!      Product ( 1 <= I <= B ) Gamma(C(I)+X(I)) / Gamma(C(I))
!
!    where:
!
!      Comb(A,B,X) is the multinomial coefficient C( A; X(1), X(2), ..., X(B) ),
!      C_Sum = Sum ( 1 <= I <= B ) C(I)
!
!  Modified:
!
!    17 December 1999
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Kenneth Lange,
!    Mathematical and Statistical Methods for Genetic Analysis,
!    Springer, 1997, page 45.
!
!  Parameters:
!
!    Input, integer X(B); X(I) counts the number of occurrences of
!    outcome I, out of the total of A trials.
!
!    Input, integer A, the total number of trials.
!
!    Input, integer B, the number of different possible outcomes on
!    one trial.
!
!    Input, integer C(B); C(I) is the Dirichlet parameter associated
!    with outcome I.
!
!    Output, real PDF, the value of the Dirichlet multinomial PDF.
!
  implicit none

  integer b

  integer a
  real c(b)
  real c_sum
  real gamma_log
  integer i
  real pdf
  real pdf_log
  integer x(b)

  c_sum = sum ( c(1:b) )

  pdf_log = - gamma_log ( c_sum + real ( a ) ) + gamma_log ( c_sum ) &
    + gamma_log ( real ( a + 1 ) )

  do i = 1, b

    pdf_log = pdf_log + gamma_log ( c(i) + real ( x(i) ) ) &
      - gamma_log ( c(i) ) - gamma_log ( real ( x(i) + 1 ) )
  end do

  pdf = exp ( pdf_log )

  return
end
subroutine dirichlet_sample ( n, a, x )

!*****************************************************************************80
!
!! DIRICHLET_SAMPLE samples the Dirichlet PDF.
!
!  Modified:
!
!    23 November 1999
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Jerry Banks, editor,
!    Handbook of Simulation,
!    Engineering and Management Press Books, 1998, page 169.
!
!  Parameters:
!
!    Input, integer N, the number of components.
!
!    Input, real A(N), the probabilities for each component.
!    Each A(I) should be nonnegative, and at least one should be
!    positive.
!
!    Output, real X(N), a sample of the PDF.  The entries of X should
!    sum to 1.
!
  implicit none

  integer n

  real a(n)
  real, parameter :: a2 = 0.0E+00
  real, parameter :: b2 = 1.0E+00
  real c2
  integer fail
  integer i
  real x(n)

  fail = 0

  do

    do i = 1, n
      c2 = a(i)
      call gamma_sample ( a2, b2, c2, x(i) )
    end do
!
!  Rescale the vector to have unit sum.
!
    if ( sum ( x(1:n) ) /= 0.0 ) then
      exit
    end if

    fail = fail + 1

    if ( fail == 1 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DIRICHLET_SAMPLE - Warning!'
      write ( *, '(a)' ) '  The Gamma samples sum to 0.'
      write ( *, '(a)' ) '  We will retry the calculation.'
      call r4vec_print ( n, a, '  The probability vector A:' )
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DIRICHLET_SAMPLE - Fatal error!'
      write ( *, '(a)' ) '  The Gamma samples sum to 0 on second try.'
      write ( *, '(a)' ) '  We abort execution.'
      stop
    end if

  end do

  call r4vec_unit_sum ( n, x )

  return
end
subroutine discrete_cdf_inv ( cdf, a, b, x )

!*****************************************************************************80
!
!! DISCRETE_CDF_INV inverts the Discrete CDF.
!
!  Modified:
!
!    05 December 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real CDF, the value of the CDF.
!    0.0 <= CDF <= 1.0.
!
!    Input, integer A, the number of probabilities assigned.
!
!    Input, real B(A), the relative probabilities of outcomes 1 through A.
!    Each entry must be nonnegative.
!
!    Output, integer X, the corresponding argument for which
!    CDF(X-1) < CDF <= CDF(X)
!
  implicit none

  integer a

  real b(a)
  real b_sum
  real cdf
  real cum
  integer j
  integer x

  if ( cdf < 0.0 .or. cdf > 1.0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DISCRETE_CDF_INV - Fatal error!'
    write ( *, '(a)' ) '  CDF < 0 or 1 < CDF.'
    stop
  end if

  b_sum = sum ( b(1:a) )

  cum = 0.0E+00

  do j = 1, a

    cum = cum + b(j) / b_sum

    if ( cdf <= cum ) then
      x = j
      return
    end if

  end do

  x = a

  return
end
subroutine discrete_sample ( a, b, x )

!*****************************************************************************80
!
!! DISCRETE_SAMPLE samples the Discrete PDF.
!
!  Modified:
!
!    05 December 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer A, the number of probabilities assigned.
!
!    Input, real B(A), the relative probabilities of outcomes 1 through A.
!    Each entry must be nonnegative.
!
!    Output, integer X, a sample of the PDF.
!
  implicit none

  integer a

  real b(a)
  real b_sum
  real cdf
  integer x

  b_sum = sum ( b(1:a) )

  call r4_random ( 0.0E+00, 1.0E+00, cdf )

  call discrete_cdf_inv ( cdf, a, b, x )

  return
end
subroutine event_process ( acid_num, alpha, beta, comp_num, p, p_hat, &
  site_num, x_sample )

!*****************************************************************************80
!
!! EVENT_PROCESS updates the mixture weight distribution parameters.
!
!  Discussion:
!
!    This routine updates the values of ALPHA.  It does this by
!    considering the results of the most recent event.  If we knew
!    which component PDF had generated the event, then we would 
!    simply add 1 to the ALPHA for that component.  Instead, we
!    use Bayesian analysis to estimate the proportion of the event
!    that is to be added to each ALPHA.
!
!  Modified:
!
!    18 December 1999
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    BS Everitt, DJ Hand,
!    Finite Mixture Distributions,
!    Chapman and Hall, 1981.
!
!    AFM Smith, UE Makov,
!    A Quasi-Bayes Sequential Procedure for Mixtures,
!    Journal of the Royal Statistical Society,
!    Volume 40, Number 1, B, 1978, pages 106-112.
!
!  Parameters:
!
!    Input, integer ACID_NUM, the number of amino acids.
!
!    Input/output, real ALPHA(COMP_NUM), the Dirichlet parameters for the 
!    weights.
!
!    Input, real BETA(ACID_NUM,COMP_MAX); BETA(I,J) is the multinomial Dirichlet
!    parameter for the J-th acid in the I-th Dirichlet mixture component.
!
!    Input, integer COMP_NUM, the number of components in the Dirichlet mixture.
!
!    Input/output, real P(COMP_NUM); P(I) is the Bayesian posterior probability 
!    of component I, given the observation of the most recent event, which is 
!    proportional to the probability of the event under the component I PDF,
!    times the prior probability of component I.
!
!    Input/output, real P_HAT(COMP_NUM), the prior probablities of the
!    components.
!
!    Input, integer SITE_NUM, the number of sites observed for this event.
!    This value might change from call to call, although in the demonstration
!    I'm keeping it fixed.
!
!    Input, integer X_SAMPLE(ACID_NUM), the "current event", namely, the count
!    vector for the number of occurrences of each acid out of the total
!    of SITE_NUM sites analyzed.  This is the evidence used to update the
!    "theory" for the value of ALPHA.
!
  implicit none

  integer acid_num
  integer comp_num

  real alpha(comp_num)
  real alpha_sum
  real beta(acid_num,comp_num)
  integer comp_i
  real comp_pdf
  real p(comp_num)
  real p_hat(comp_num)
  integer site_num
  integer x_sample(acid_num)
!
!  Sum the parameters.
!
  alpha_sum = sum ( alpha(1:comp_num) )
!
!  Update P_HAT.
!
  do comp_i = 1, comp_num
    p_hat(comp_i) = ( ( alpha_sum - 1.0E+00 ) * p_hat(comp_i) + p(comp_i) ) &
      / alpha_sum
  end do
!
!  Generate the new P's.
!  P(COMP_I) = the Bayesian posterior probability of component I,
!  given the observation of event EVENT_I, which is proportional
!  to the probability of event EVENT_I in the component I PDF,
!  times the prior probability of component I.
!
  do comp_i = 1, comp_num
!
!  I think this is where I need the infamous Dirichlet-multinomial PDF.
!
    call dirichlet_multinomial_pdf ( x_sample, site_num, acid_num, &
      beta(1,comp_i), comp_pdf )

    p(comp_i) = comp_pdf * p_hat(comp_i)

  end do
!
!  Normalize the P's.
!
  if ( sum ( p(1:comp_num) ) == 0.0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'EVENT_PROCESS - Fatal error!'
    write ( *, '(a)' ) '  The P''s sum to 0.'
    stop
  end if

  call r4vec_unit_sum ( comp_num, p )
!
!  Update the alpha's by adding adding appropriate portions of
!  the most recent event to each component's parameter.
!
  do comp_i = 1, comp_num
    alpha(comp_i) = alpha(comp_i) + p(comp_i)
  end do

  return
end
subroutine exponential_01_cdf_inv ( cdf, x )

!*****************************************************************************80
!
!! EXPONENTIAL_01_CDF_INV inverts the Exponential 01 CDF.
!
!  Modified:
!
!    09 December 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real CDF, the value of the CDF.
!    0.0 <= CDF <= 1.0.
!
!    Output, real X, the corresponding argument.
!
  implicit none

  real cdf
  real x

  if ( cdf < 0.0E+00 .or. cdf > 1.0E+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'EXPONENTIAL_01_CDF_INV - Fatal error!'
    write ( *, '(a)' ) '  CDF < 0 or 1 < CDF.'
    stop
  end if

  x = - log ( 1.0E+00 - cdf )

  return
end
subroutine exponential_01_sample ( x )

!*****************************************************************************80
!
!! EXPONENTIAL_01_SAMPLE samples the Exponential PDF with parameter 1.
!
!  Modified:
!
!    16 June 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real X, a sample of the PDF.
!
  implicit none

  real cdf
  real x

  call r4_random ( 0.0E+00, 1.0E+00, cdf )

  call exponential_01_cdf_inv ( cdf, x )

  return
end
subroutine exponential_cdf_inv ( cdf, a, b, x )

!*****************************************************************************80
!
!! EXPONENTIAL_CDF_INV inverts the Exponential CDF.
!
!  Modified:
!
!    01 May 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real CDF, the value of the CDF.
!    0.0 <= CDF <= 1.0.
!
!    Input, real A, B, the parameters of the PDF.
!    0.0 < B.
!
!    Output, real X, the corresponding argument.
!
  implicit none

  real a
  real b
  real cdf
  real x

  if ( cdf < 0.0E+00 .or. cdf > 1.0E+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'EXPONENTIAL_CDF_INV - Fatal error!'
    write ( *, '(a)' ) '  CDF < 0 or 1 < CDF.'
    stop
  end if

  x = a - b * log ( 1.0E+00 - cdf )

  return
end
subroutine exponential_sample ( a, b, x )

!*****************************************************************************80
!
!! EXPONENTIAL_SAMPLE samples the Exponential PDF.
!
!  Modified:
!
!    01 May 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real A, B, the parameters of the PDF.
!    0.0 < B.
!
!    Output, real X, a sample of the PDF.
!
  implicit none

  real a
  real b
  real cdf
  real x

  call r4_random ( 0.0E+00, 1.0E+00, cdf )

  call exponential_cdf_inv ( cdf, a, b, x )

  return
end
function gamma_log ( x )

!*****************************************************************************80
!
!! GAMMA_LOG calculates the natural logarithm of GAMMA ( X ) for positive X.
!
!  Discussion:
!
!    The program uses rational functions that theoretically approximate
!    LOG(GAMMA(X)) to at least 18 significant decimal digits.  
!
!  Modified:
!
!    16 June 1999
!
!  Author:
!
!    William Cody, Laura Stoltz
!    FORTRAN90 version by John Burkardt
!
!  References:
!
!    W. J. Cody and K. E. Hillstrom,
!    Chebyshev Approximations for the Natural Logarithm of the Gamma Function,
!    Math. Comp.
!    Volume 21, 1967, pages 198-203.
!
!    K. E. Hillstrom,
!    ANL/AMD Program ANLC366S, DGAMMA/DLGAMA,
!    May 1969.
!
!    Hart, Et. Al.,
!    Computer Approximations,
!    Wiley and sons, New York, 1968.
!
!  Parameters:
!
!    Input, real X, the argument of the Gamma function.  X must be positive.
!
!    Output, real GAMMA_LOG, the logarithm of the Gamma function of X.
!
!  Explanation of machine-dependent constants
!
!  BETA   - radix for the floating-point representation.
!
!  MAXEXP - the smallest positive power of BETA that overflows.
!
!  XBIG   - largest argument for which LN(GAMMA(X)) is representable
!           in the machine, i.e., the solution to the equation
!             LN(GAMMA(XBIG)) = BETA**MAXEXP.
!
!  FRTBIG - Rough estimate of the fourth root of XBIG
!
!
!  Approximate values for some important machines are:
!
!                            BETA      MAXEXP         XBIG
!
!  CRAY-1        (S.P.)        2        8191       9.62E+2461
!  Cyber 180/855
!    under NOS   (S.P.)        2        1070       1.72E+319
!  IEEE (IBM/XT,
!    SUN, etc.)  (S.P.)        2         128       4.08E+36
!  IEEE (IBM/XT,
!    SUN, etc.)  (D.P.)        2        1024       2.55D+305
!  IBM 3033      (D.P.)       16          63       4.29D+73
!  VAX D-Format  (D.P.)        2         127       2.05D+36
!  VAX G-Format  (D.P.)        2        1023       1.28D+305
!
!
!                           FRTBIG
!
!  CRAY-1        (S.P.)   3.13E+615
!  Cyber 180/855
!    under NOS   (S.P.)   6.44E+79
!  IEEE (IBM/XT,
!    SUN, etc.)  (S.P.)   1.42E+9
!  IEEE (IBM/XT,
!    SUN, etc.)  (D.P.)   2.25D+76
!  IBM 3033      (D.P.)   2.56D+18
!  VAX D-Format  (D.P.)   1.20D+9
!  VAX G-Format  (D.P.)   1.89D+76
!
  implicit none

  real, parameter, dimension ( 7 ) :: c = (/ &
    -1.910444077728E-03, &
     8.4171387781295E-04, &
    -5.952379913043012E-04, &
     7.93650793500350248E-04, &
    -2.777777777777681622553E-03, &
     8.333333333333333331554247E-02, &
     5.7083835261e-03 /)
  real corr
  real, parameter :: d1 = - 5.772156649015328605195174E-01
  real, parameter :: d2 =   4.227843350984671393993777E-01
  real, parameter :: d4 =   1.791759469228055000094023E+00
  integer i
  real, parameter :: frtbig = 1.42E+09
  real gamma_log
  real, parameter, dimension ( 8 ) :: p1 = (/ &
    4.945235359296727046734888E+00, &
    2.018112620856775083915565E+02, &
    2.290838373831346393026739E+03, &
    1.131967205903380828685045E+04, &
    2.855724635671635335736389E+04, &
    3.848496228443793359990269E+04, &
    2.637748787624195437963534E+04, &
    7.225813979700288197698961E+03 /)
  real, parameter, dimension ( 8 ) :: p2 = (/ &
    4.974607845568932035012064E+00, &
    5.424138599891070494101986E+02, &
    1.550693864978364947665077E+04, &
    1.847932904445632425417223E+05, &
    1.088204769468828767498470E+06, &
    3.338152967987029735917223E+06, &
    5.106661678927352456275255E+06, &
    3.074109054850539556250927E+06 /)
  real, parameter, dimension ( 8 ) :: p4 = (/ &
    1.474502166059939948905062E+04, &
    2.426813369486704502836312E+06, &
    1.214755574045093227939592E+08, &
    2.663432449630976949898078E+09, &
    2.940378956634553899906876E+10, &
    1.702665737765398868392998E+11, &
    4.926125793377430887588120E+11, &
    5.606251856223951465078242E+11 /)
  real, parameter :: pnt68 = 0.6796875E+00
  real, parameter, dimension ( 8 ) :: q1 = (/ &
    6.748212550303777196073036E+01, &
    1.113332393857199323513008E+03, &
    7.738757056935398733233834E+03, &
    2.763987074403340708898585E+04, &
    5.499310206226157329794414E+04, &
    6.161122180066002127833352E+04, &
    3.635127591501940507276287E+04, &
    8.785536302431013170870835E+03 /)
  real, parameter, dimension ( 8 ) :: q2 = (/ &
    1.830328399370592604055942E+02, &
    7.765049321445005871323047E+03, &
    1.331903827966074194402448E+05, &
    1.136705821321969608938755E+06, &
    5.267964117437946917577538E+06, &
    1.346701454311101692290052E+07, &
    1.782736530353274213975932E+07, &
    9.533095591844353613395747E+06 /)
  real, parameter, dimension ( 8 ) :: q4 = (/ &
    2.690530175870899333379843E+03, &
    6.393885654300092398984238E+05, &
    4.135599930241388052042842E+07, &
    1.120872109616147941376570E+09, &
    1.488613728678813811542398E+10, &
    1.016803586272438228077304E+11, &
    3.417476345507377132798597E+11, &
    4.463158187419713286462081E+11 /)
  real res
  real, parameter :: sqrtpi = 0.9189385332046727417803297E+00
  real x
  real, parameter :: xbig = 4.08E+36
  real xden
  real xm1
  real xm2
  real xm4
  real xnum
  real xsq
!
!  Return immediately if the argument is out of range.
!
  if ( x <= 0.0E+00 .or. x > xbig ) then
    gamma_log = huge ( gamma_log )
    return
  end if

  if ( x <= epsilon ( x ) ) then

    res = - log ( x )

  else if ( x <= 1.5E+00 ) then

    if ( x < pnt68 ) then
      corr = - log ( x )
      xm1 = x
    else
      corr = 0.0E+00
      xm1 = ( x - 0.5E+00 ) - 0.5E+00
    end if

    if ( x <= 0.5E+00 .or. x >= pnt68 ) then

      xden = 1.0E+00
      xnum = 0.0E+00

      do i = 1, 8
        xnum = xnum * xm1 + p1(i)
        xden = xden * xm1 + q1(i)
      end do

      res = corr + ( xm1 * ( d1 + xm1 * ( xnum / xden ) ) )

    else

      xm2 = ( x - 0.5E+00 ) - 0.5E+00
      xden = 1.0E+00
      xnum = 0.0E+00
      do i = 1, 8
        xnum = xnum * xm2 + p2(i)
        xden = xden * xm2 + q2(i)
      end do

      res = corr + xm2 * ( d2 + xm2 * ( xnum / xden ) )

    end if

  else if ( x <= 4.0E+00 ) then

    xm2 = x - 2.0E+00
    xden = 1.0E+00
    xnum = 0.0E+00
    do i = 1, 8
      xnum = xnum * xm2 + p2(i)
      xden = xden * xm2 + q2(i)
    end do

    res = xm2 * ( d2 + xm2 * ( xnum / xden ) )

  else if ( x <= 12.0E+00 ) then

    xm4 = x - 4.0E+00
    xden = - 1.0E+00
    xnum = 0.0E+00
    do i = 1, 8
      xnum = xnum * xm4 + p4(i)
      xden = xden * xm4 + q4(i)
    end do

    res = d4 + xm4 * ( xnum / xden )

  else

    res = 0.0E+00

    if ( x <= frtbig ) then

      res = c(7)
      xsq = x * x

      do i = 1, 6
        res = res / xsq + c(i)
      end do

    end if

    res = res / x
    corr = log ( x )
    res = res + sqrtpi - 0.5E+00 * corr
    res = res + x * ( corr - 1.0E+00 )

  end if

  gamma_log = res

  return
end
subroutine gamma_sample ( a, b, c, x )

!*****************************************************************************80
!
!! GAMMA_SAMPLE samples the Gamma PDF.
!
!  Modified:
!
!    15 September 2000
!
!  Reference:
!
!    JH Ahrens, U Dieter,
!    Generating Gamma Variates by a Modified Rejection Technique,
!    Communications of the ACM,
!    Volume 25, Number 1, January 1982, pages 47 - 54.
!
!    JH Ahrens, U Dieter,
!    Computer Methods for Sampling from Gamma, Beta, Poisson and
!    Binomial Distributions.
!    Computing, Volume 12, 1974, pages 223 - 246.
!
!    JH Ahrens, KD Kohrt, U Dieter,
!    Algorithm 599,
!    ACM Transactions on Mathematical Software,
!    Volume 9, Number 2, June 1983, pages 255-257.
!
!  Parameters:
!
!    Input, real A, B, C, the parameters of the PDF.
!    0.0E+00 < B,
!    0.0E+00 < C.
!
!    Output, real X, a sample of the PDF.
!
  implicit none

  real a
  real, parameter :: a1 =   0.3333333E+00
  real, parameter :: a2 = - 0.2500030E+00
  real, parameter :: a3 =   0.2000062E+00
  real, parameter :: a4 = - 0.1662921E+00
  real, parameter :: a5 =   0.1423657E+00
  real, parameter :: a6 = - 0.1367177E+00
  real, parameter :: a7 =   0.1233795E+00
  real b
  real bcoef
  real c
  real co
  real d
  real e
  real, parameter :: e1 = 1.0E+00
  real, parameter :: e2 = 0.4999897E+00
  real, parameter :: e3 = 0.1668290E+00
  real, parameter :: e4 = 0.0407753E+00
  real, parameter :: e5 = 0.0102930E+00
  real, parameter :: euler = 2.71828182845904E+00
  real p
  real q
  real q0
  real, parameter :: q1 =   0.04166669E+00
  real, parameter :: q2 =   0.02083148E+00
  real, parameter :: q3 =   0.00801191E+00
  real, parameter :: q4 =   0.00144121E+00
  real, parameter :: q5 = - 0.00007388E+00
  real, parameter :: q6 =   0.00024511E+00
  real, parameter :: q7 =   0.00024240E+00
  real r
  real s
  real si
  real s2
  real t
  real u
  real v
  real w
  real x
!
!  Allow C = 0.
!
  if ( c == 0.0E+00 ) then
    x = a
    return
  end if
!
!  C < 1.
!
  if ( c < 1.0E+00 ) then

    do

      call r4_random ( 0.0E+00, 1.0E+00, u )
      t = 1.0E+00 + c / euler
      p = u * t

      call exponential_01_sample ( s )

      if ( p < 1.0E+00 ) then
        x = exp ( log ( p ) / c )
        if ( x <= s ) then
          exit
        end if
      else
        x = - log ( ( t - p ) / c )
        if ( ( 1.0E+00 - c ) * log ( x ) <= s ) then
          exit
        end if
      end if

    end do

    x = a + b * x
    return
!
!  1 <= C.
!
  else

    s2 = c - 0.5E+00
    s = sqrt ( c - 0.5E+00 )
    d = sqrt ( 32.0E+00 ) - 12.0E+00 * sqrt ( c - 0.5E+00 )

    call normal_01_sample ( t )
    x = ( sqrt ( c - 0.5E+00 ) + 0.5E+00 * t )**2

    if ( t >= 0.0E+00 ) then
      x = a + b * x
      return
    end if

    call r4_random ( 0.0E+00, 1.0E+00, u )

    if ( d * u <= t**3 ) then
      x = a + b * x
      return
    end if

    r = 1.0E+00 / c

    q0 = ( ( ( ( ( ( &
           q7   * r &
         + q6 ) * r &
         + q5 ) * r &
         + q4 ) * r &
         + q3 ) * r &
         + q2 ) * r &
         + q1 ) * r

    if ( c <= 3.686E+00 ) then
      bcoef = 0.463E+00 + s - 0.178E+00 * s2
      si = 1.235E+00
      co = 0.195E+00 / s - 0.079E+00 + 0.016 * s
    else if ( c <= 13.022E+00 ) then
      bcoef = 1.654E+00 + 0.0076E+00 * s2
      si = 1.68E+00 / s + 0.275E+00
      co = 0.062E+00 / s + 0.024E+00
    else
      bcoef = 1.77E+00
      si = 0.75E+00
      co = 0.1515E+00 / s
    end if

    if ( sqrt ( c - 0.5E+00 ) + 0.5E+00 * t > 0.0E+00 ) then

      v = 0.5E+00 * t / s

      if ( abs ( v ) > 0.25E+00 ) then
        q = q0 - s * t + 0.25E+00 * t**2 + 2.0E+00 * s2 * log ( 1.0E+00 + v )
      else
        q = q0 + 0.5E+00 * t**2 * ( ( ( ( ( ( &
               a7   * v &
             + a6 ) * v &
             + a5 ) * v &
             + a4 ) * v &
             + a3 ) * v &
             + a2 ) * v &
             + a1 ) * v
      end if

      if ( log ( 1.0E+00 - u ) <= q ) then
        x = a + b * x
        return
      end if

    end if

    do

      call exponential_01_sample ( e )

      call r4_random ( 0.0E+00, 1.0E+00, u )
      u = 2.0E+00 * u - 1.0E+00
      t = bcoef + sign ( si * e, u )

      if ( t >= - 0.7187449E+00 ) then

        v = 0.5E+00 * t / s

        if ( abs ( v ) > 0.25E+00 ) then
          q = q0 - s * t + 0.25 * t**2 + 2.0E+00 * s2 * log ( 1.0E+00 + v )
        else
          q = q0 + 0.5E+00 * t**2 * ( ( ( ( ( ( &
               a7   * v &
             + a6 ) * v &
             + a5 ) * v &
             + a4 ) * v &
             + a3 ) * v &
             + a2 ) * v &
             + a1 ) * v
        end if

        if ( q > 0.0E+00 ) then

          if ( q > 0.5E+00 ) then
            w = exp ( q ) - 1.0E+00
          else
            w = ( ( ( ( &
                    e5   * q &
                  + e4 ) * q &
                  + e3 ) * q &
                  + e2 ) * q &
                  + e1 ) * q
          end if

          if ( co * abs ( u ) <= w * exp ( e - 0.5E+00 * t**2 ) ) then
            x = a + b * ( s + 0.5E+00 * t )**2
            return
          end if

        end if

      end if

    end do

  end if

  return
end
subroutine i4_next ( line, ival, done )

!*****************************************************************************80
!
!! I4_NEXT "reads" integers from a string, one at a time.
!
!  Modified:
!
!    14 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) LINE, a string, presumably containing
!    integers.  These may be separated by spaces or commas.
!
!    Output, integer IVAL.  If DONE is FALSE, then IVAL contains the
!    "next" integer read from LINE.  If DONE is TRUE, then
!    IVAL is zero.
!
!    Input/output, logical DONE.
!    On input with a fresh value of LINE, the user should set
!    DONE to TRUE.
!    On output, the routine sets DONE to FALSE if another integer
!    was read, or TRUE if no more integers could be read.
!
  implicit none

  logical done
  integer ierror
  integer ival
  integer lchar
  character ( len = * ) line
  integer, save :: next = 1

  ival = 0

  if ( done ) then
    next = 1
    done = .false.
  end if

  if ( next > len(line) ) then
    done = .true.
    return
  end if

  call s_to_i4 ( line(next:), ival, ierror, lchar )

  if ( ierror /= 0 .or. lchar == 0 ) then
    done = .true.
    next = 1
  else
    done = .false.
    next = next + lchar
  end if

  return
end
subroutine mixture_read ( acid_num, acid_sym, beta, beta_sum, comp_label, &
  comp_max, comp_num, comp_weight, ierror, iunit )

!*****************************************************************************80
!
!! MIXTURE_READ reads the Dirichlet mixture parameters from a file.
!
!  Discussion:
!
!    The data in the file is delimited by keywords.
!
!    The first lines (not necessarily in order!) may include
!
!      ClassName = string
!      NumDistr = N           the number of components in the mixture.
!      Alphabet = string
!      Order = A C D E ...    the order of the amino acids.
!      AlphaChar = 20
!      NumDistr = 9           the number of distributions
!      EndClassName = string
!
!    For each component, there are four lines:
!
!      Number= N              the component number, starting with 0
!      Mixture= N             the mixture weight, out of a total of 1.0
!      Alpha=  |A| A1 A2 ...  the parameter sum, and individual parameters
!      Comment=               a comment, which describes the frequencies.
!
!    In the comment, the symbol "><" indicates the mean background frequency;
!    residues to the left of that symbol occur more frequently
!    than background, residues to the right less frequently.  Commas separate
!    residues differing in frequency by a factor of 2.
!
!    For example, the comment
!      S A T , C G P >< N V M , Q H R I K F L D W , E Y
!    indicates that for this component, the frequency of
!    proline is just above the mean, and serine, alanine and
!    threonine are twice as frequent in this component than they
!    are on average.  By contrast, tyrosine and glutamic acid are
!    between 4 and 8 times less likely in this component than on
!    average.
!
!  Modified:
!
!    21 January 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ACID_NUM, the number of amino acids.
!
!    Output, character ACID_SYM(ACID_NUM), the one letter amino acid codes.
!
!    Output, real BETA(ACID_NUM,COMP_MAX); BETA(I,J) is the parameter
!    for the J-th acid in the I-th Dirichlet mixture component.
!
!    Output, real BETA_SUM(COMP_MAX), the sum of the values of
!    BETA(ACID_I,COMP_I) for a given component COMP_I.
!
!    Output, integer COMP_LABEL(COMP_NUM), the label of each component.
!    Normally, component I has label I.
!
!    Input, integer COMP_MAX, the maximum number of Dirichlet mixture
!    components.
!
!    Output, integer COMP_NUM, the number of components in the Dirichlet
!    mixture.
!
!    Output, real COMP_WEIGHT(COMP_NUM), the mixture weight of each component.
!    These values should be nonnegative, and sum to 1.  They represent the
!    relative proportion of each component in the mixture.
!
!    Output, integer IERROR, error indicator.
!    0: no error occurred; nonzero: an error occurred.
!
!    Input, integer IUNIT, the FORTRAN unit from which the data is to
!    be read.
!
  implicit none

  integer acid_num
  integer comp_max

  integer acid_i
  character acid_sym(acid_num)
  real beta(acid_num,comp_max)
  real beta_sum(comp_max)
  integer comp_i
  integer comp_label(comp_max)
  integer comp_num
  real comp_weight(comp_max)
  logical done
  integer iequal
  integer ierror
  integer iunit
  integer ngoofy
  integer nrec
  logical s_begin
  character ( len = 500 ) string

  ierror = 0
  comp_i = 0
  comp_num = 0
  nrec = 0
  ngoofy = 0

10    continue

  read ( iunit, '(a)', end = 20 ) string
  nrec = nrec + 1
!
!  Ignore blank lines.
!
  if ( string == ' ' ) then
!
!  Ignore the CLASSNAME field.
!
  else if ( s_begin ( string, 'CLASSNAME' ) ) then
!
!  Ignore the ENDCLASSNAME field.
!
  else if ( s_begin ( string, 'ENDCLASSNAME' ) ) then
!
!  Ignore the NAME field.
!
  else if ( s_begin ( string, 'NAME' ) ) then
!
!  Ignore the ALPHABET field.
!
  else if ( s_begin ( string, 'ALPHABET' ) ) then
!
!  Read the ORDER field, since it tells us how to interpret the ALPHA's.
!
  else if ( s_begin ( string, 'ORDER' ) ) then

    iequal = index ( string, '=' )
    done = .true.
    do acid_i = 1, acid_num
      call ch_next ( string(iequal+1:), acid_sym(acid_i), done )
    end do
!
!  Ignore the ALPHACHAR field.
!
  else if ( s_begin ( string, 'ALPHACHAR' ) ) then
!
!  Read the NUMDISTR field.
!
  else if ( s_begin ( string, 'NUMDISTR' ) ) then

    iequal = index ( string, '=' )
    done = .true.
    call i4_next ( string(iequal+1:), comp_num, done )

    if ( comp_num < 1 ) then
      ierror = 1
      return
    else if ( comp_num > comp_max ) then
      ierror = 2
      return
    end if
!
!  Read the NUMBER field.
!
  else if ( s_begin ( string, 'NUMBER' ) ) then

    comp_i = comp_i + 1

    if ( comp_i > comp_num ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'MIXTURE_READ - Fatal error!'
      write ( *, '(a,i6)' ) '  Number of components = ', comp_i
      write ( *, '(a,i6)' ) '  exceeding reported value of ', comp_num
      stop
    end if

    iequal = index ( string, '=' )
    done = .true.
    call i4_next ( string(iequal+1:), comp_label(comp_i), done )
!
!  Read the MIXTURE field.
!
  else if ( s_begin ( string, 'MIXTURE' ) ) then

    iequal = index ( string, '=' )
    done = .true.
    call r4_next ( string(iequal+1:), comp_weight(comp_i), done )
!
!  Read the ALPHA field.
!
  else if ( s_begin ( string, 'ALPHA' ) ) then

    iequal = index ( string, '=' )
    done = .true.
    call r4_next ( string(iequal+1:), beta_sum(comp_i), done )
    do acid_i = 1, acid_num
      call r4_next ( string(iequal+1:), beta(acid_i,comp_i), done )
    end do
!
!  Ignore the COMMENT field.
!
  else if ( s_begin ( string, 'COMMENT' ) ) then
!
!  Unexpected field:
!
  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'MIXTURE_READ - Warning!'
    write ( *, '(a)' ) '  Goofy record: '
    write ( *, '(a)' ) trim ( string )

    ngoofy = ngoofy + 1

  end if

  go to 10

20    continue

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'MIXTURE_READ - Note:'
  write ( *, '(a,i6)' ) '  Number of records read was ', nrec
  write ( *, '(a,i6)' ) '  Number of goofy records was ', ngoofy

  return
end
subroutine multinomial_sample ( a, b, c, x )

!*****************************************************************************80
!
!! MULTINOMIAL_SAMPLE samples the Multinomial PDF.
!
!  Modified:
!
!    14 February 1999
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Luc Devroye,
!    Non-Uniform Random Variate Generation,
!    Springer-Verlag, New York, 1986, page 559.
!
!  Parameters:
!
!    Input, integer A, the total number of trials.
!    0 <= A.
!
!    Input, integer B, the number of outcomes possible on one trial.
!    1 <= B.
!
!    Input, real C(B).  C(I) is the probability of outcome I on
!    any trial.
!    0.0 <= C(I) <= 1.0,
!    SUM ( 1 <= I <= B) C(I) = 1.0.
!
!    Output, integer X(B); X(I) is the number of
!    occurrences of event I during the N trials.
!
  implicit none

  integer b

  integer a
  real c(b)
  integer i
  integer ifactor
  integer ntot
  real prob
  real sum2
  integer x(b)

  ntot = a

  sum2 = 1.0E+00

  x(1:b) = 0

  do ifactor = 1, b - 1

    prob = c(ifactor) / sum2
!
!  Generate a binomial random deviate for NTOT trials with
!  single trial success probability PROB.
!
    call binomial_sample ( ntot, prob, x(ifactor) )

    ntot = ntot - x(ifactor)
    if ( ntot <= 0 ) then
      return
    end if

    sum2 = sum2 - c(ifactor)

  end do
!
!  The last factor gets what's left.
!
  x(b) = ntot

  return
end
subroutine normal_01_sample ( x )

!*****************************************************************************80
!
!! NORMAL_01_SAMPLE samples the Normal 01 PDF.
!
!  Discussion:
!
!    The standard normal distribution has mean 0 and standard
!    deviation 1.
!
!    The Box-Muller method is used.
!
!  Modified:
!
!    10 February 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real X, a sample of the PDF.
!
  implicit none

  integer, save :: iset = 0
  real, parameter :: PI = &
    3.14159265358979323846264338327950288419716939937510E+00
  real v1
  real v2
  real x
  real, save :: xsave = 0.0E+00

  if ( iset == 0 ) then

    call r4_random ( 0.0E+00, 1.0E+00, v1 )

    if ( v1 <= 0.0E+00 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'NORMAL_01_SAMPLE - Fatal error!'
      write ( *, '(a)' ) '  V1 <= 0.'
      write ( *, '(a,g14.6)' ) '  V1 = ', v1
      stop
    end if

    call r4_random ( 0.0E+00, 1.0E+00, v2 )

    if ( v2 <= 0.0E+00 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'NORMAL_01_SAMPLE - Fatal error!'
      write ( *, '(a)' ) '  V2 <= 0.'
      write ( *, '(a,g14.6)' ) '  V2 = ', v2
      stop
    end if

    x = sqrt ( - 2.0E+00 * log ( v1 ) ) * cos ( 2.0E+00 * PI * v2 )

    xsave = sqrt ( - 2.0E+00 * log ( v1 ) ) * sin ( 2.0E+00 * PI * v2 )

    iset = 1

  else

    x = xsave
    iset = 0

  end if

  return
end
subroutine r4_next ( line, rval, done )

!*****************************************************************************80
!
!! R4_NEXT "reads" real numbers from a string, one at a time.
!
!  Modified:
!
!    14 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) LINE, a string, presumably containing real
!    numbers.  These may be separated by spaces or commas.
!
!    Output, real RVAL.  If DONE is FALSE, then RVAL contains the
!    "next" real value read from LINE.  If DONE is TRUE, then
!    RVAL is zero.
!
!    Input/output, logical DONE.
!    On input with a fresh value of LINE, the user should set
!    DONE to TRUE.
!    On output, the routine sets DONE to FALSE if another real
!    value was read, or TRUE if no more reals could be read.
!
  implicit none

  logical done
  integer ierror
  integer lchar
  character ( len = * ) line
  integer, save :: next = 1
  real rval

  rval = 0.0E+00

  if ( done ) then
    next = 1
    done = .false.
  end if

  if ( next > len(line) ) then
    done = .true.
    return
  end if

  call s_to_r4 ( line(next:), rval, ierror, lchar )

  if ( ierror /= 0 .or. lchar == 0 ) then
    done = .true.
    next = 1
  else
    done = .false.
    next = next + lchar
  end if

  return
end
subroutine r4_random ( rlo, rhi, r )

!*****************************************************************************80
!
!! R4_RANDOM returns a random real in a given range.
!
!  Modified:
!
!    01 December 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real RLO, RHI, the minimum and maximum values.
!
!    Output, real R, the randomly chosen value.
!
  implicit none

  logical, save :: seed = .false.
  real r
  real rhi
  real rlo
  real t

  if ( .not. seed ) then
    call random_seed ( )
    seed = .true.
  end if
!
!  Pick a random number in (0,1).
!
  call random_number ( harvest = t )
!
!  Set R.
!
  r = ( 1.0 - t ) * rlo + t * rhi

  return
end
subroutine r4vec_print ( n, a, title )

!*****************************************************************************80
!
!! R4VEC_PRINT prints a real vector, with an optional title.
!
!  Modified:
!
!    16 December 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of components of the vector.
!
!    Input, real A(N), the vector to be printed.
!
!    Input, character ( len = * ) TITLE, a title to be printed first.
!    TITLE may be blank.
!
  implicit none

  integer n

  real a(n)
  integer i
  character ( len = * ) title

  if ( title /= ' ' ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) title
  end if

  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(i6,g14.6)' ) i, a(i)
  end do

  return
end
subroutine r4vec_unit_sum ( n, a )

!*****************************************************************************80
!
!! R4VEC_UNIT_SUM normalizes a real vector to have unit sum.
!
!  Modified:
!
!    28 February 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input/output, real A(N), the vector to be normalized.  On output,
!    the entries of A should have unit sum.  However, if the input vector
!    has zero sum, the routine halts.
!
  implicit none

  integer n

  real a(n)
  real a_sum

  a_sum = sum ( a(1:n) )

  if ( a_sum == 0.0E+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'RVEC_UNIT_SUM - Fatal error!'
    write ( *, '(a)' ) '  The vector entries sum to 0.'
    stop
  end if

  a(1:n) = a(1:n) / a_sum

  return
end
function s_begin ( s1, s2 )

!*****************************************************************************80
!
!! S_BEGIN is TRUE if one string matches the beginning of the other.
!
!  Examples:
!
!     S1              S2      S_BEGIN
!
!    'Bob'          'BOB'     TRUE
!    '  B  o b '    ' bo b'   TRUE
!    'Bob'          'Bobby'   TRUE
!    'Bobo'         'Bobb'    FALSE
!    ' '            'Bob'     FALSE    (Do not allow a blank to match
!                                       anything but another blank string.)
!
!  Modified:
!
!    20 January 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S1, S2, the strings to be compared.
!
!    Output, logical S_BEGIN, .TRUE. if the strings match up to
!    the end of the shorter string, ignoring case,
!    FALSE otherwise.
!
  implicit none

  logical ch_eqi
  integer i1
  integer i2
  integer len1
  integer len2
  logical s_begin
  character ( len = * ) s1
  character ( len = * ) s2

  len1 = len_trim ( s1 )
  len2 = len_trim ( s2 )
!
!  If either string is blank, then both must be blank to match.
!  Otherwise, a blank string matches anything, which is not
!  what most people want.
!
  if ( len1 == 0 .or. len2 == 0 ) then

    if ( len1 == 0 .and. len2 == 0 ) then
      s_begin = .true.
    else
      s_begin = .false.
    end if

    return

  end if

  i1 = 0
  i2 = 0
!
!  Find the next nonblank in S1.
!
10    continue

  i1 = i1 + 1

  if ( i1 > len1 ) then
    s_begin = .true.
    return
  end if

  if ( s1(i1:i1) == ' ' ) then
    go to 10
  end if
!
!  Find the next nonblank in S2.
!
20    continue

  i2 = i2 + 1

  if ( i2 > len2 ) then
    s_begin = .true.
    return
  end if

  if ( s2(i2:i2) == ' ' ) then
    go to 20
  end if
!
!  If the characters match, get the next pair.
!
  if ( ch_eqi ( s1(i1:i1), s2(i2:i2) ) ) then
    go to 10
  end if

  s_begin = .false.

  return
end
subroutine s_to_i4 ( string, ival, ierror, last )

!*****************************************************************************80
!
!! S_TO_I4 reads an integer value from a string.
!
!  Modified:
!
!    28 August 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) STRING, a string to be examined.
!
!    Output, integer IVAL, the integer value read from the string.
!    If STRING is blank, then IVAL will be returned 0.
!
!    Output, integer IERROR, an error flag.
!    0, no error.
!    1, an error occurred.
!
!    Output, integer LAST, the last character in STRING that was
!    part of the representation of IVAL.
!
  implicit none

  character c
  integer i
  integer ierror
  integer isgn
  integer istate
  integer ival
  integer last
  integer lens
  character ( len = * ) string

  ierror = 0
  istate = 0

  isgn = 1
  ival = 0

  lens = len ( string )

  i = 0

10    continue

  i = i + 1

  c = string(i:i)

  if ( istate == 0 ) then

    if ( c == ' ' ) then

    else if ( c == '-' ) then
      istate = 1
      isgn = -1
    else if ( c == '+' ) then
      istate = 1
      isgn = + 1
    else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
      istate = 2
      ival = ichar ( c ) - ichar ( '0' )
    else
      ierror = 1
      return
    end if

  else if ( istate == 1 ) then

    if ( c == ' ' ) then

    else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
      istate = 2
      ival = ichar ( c ) - ichar ( '0' )
    else
      ierror = 1
      return
    end if

  else if ( istate == 2 ) then

    if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
      ival = 10 * ival + ichar ( c ) - ichar ( '0' )
    else
      istate = 3
    end if

  end if
!
!  Continue or exit?
!
  if ( istate == 3 ) then
    ival = isgn * ival
    last = i - 1
    return
  else if ( i >= lens ) then
    if ( istate == 2 ) then
      ival = isgn * ival
      last = lens
    else
      ierror = 1
      last = 0
    end if
    return
  end if

  go to 10
end
subroutine s_to_r4 ( string, rval, ierror, lchar )

!*****************************************************************************80
!
!! S_TO_R4 reads a real number from a string.
!
!  Discussion:
!
!    This routine will read as many characters as possible until it reaches
!    the end of the string, or encounters a character which cannot be
!    part of the real number.
!
!    Legal input is:
!
!       1 blanks,
!       2 '+' or '-' sign,
!       2.5 spaces
!       3 integer part,
!       4 decimal point,
!       5 fraction part,
!       6 'E' or 'e' or 'D' or 'd', exponent marker,
!       7 exponent sign,
!       8 exponent integer part,
!       9 exponent decimal point,
!      10 exponent fraction part,
!      11 blanks,
!      12 final comma or semicolon.
!
!    with most quantities optional.
!
!  Examples:
!
!    STRING            RVAL
!
!    '1'               1.0
!    '     1   '       1.0
!    '1A'              1.0
!    '12,34,56'        12.0
!    '  34 7'          34.0
!    '-1E2ABCD'        -100.0
!    '-1X2ABCD'        -1.0
!    ' 2E-1'           0.2
!    '23.45'           23.45
!    '-4.2E+2'         -420.0
!    '17d2'            1700.0
!    '-14e-2'         -0.14
!    'e2'              100.0
!    '-12.73e-9.23'   -12.73 * 10.0**(-9.23)
!
!  Modified:
!
!    14 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) STRING, the string containing the
!    data to be read.  Reading will begin at position 1 and
!    terminate at the end of the string, or when no more
!    characters can be read to form a legal real.  Blanks,
!    commas, or other nonnumeric data will, in particular,
!    cause the conversion to halt.
!
!    Output, real RVAL, the real value that was read from the string.
!
!    Output, integer IERROR, error flag.
!
!    0, no errors occurred.
!
!    1, 2, 6 or 7, the input number was garbled.  The
!    value of IERROR is the last type of input successfully
!    read.  For instance, 1 means initial blanks, 2 means
!    a plus or minus sign, and so on.
!
!    Output, integer LCHAR, the number of characters read from
!    STRING to form the number, including any terminating
!    characters such as a trailing comma or blanks.
!
  implicit none

  logical ch_eqi
  character chrtmp
  integer ierror
  integer ihave
  integer isgn
  integer iterm
  integer jbot
  integer jsgn
  integer jtop
  integer lchar
  integer nchar
  integer ndig
  real rbot
  real rexp
  real rtop
  real rval
  character ( len = * ) string
  character, parameter :: TAB = char ( 9 )

  nchar = len ( string )
  ierror = 0
  rval = 0.0E+00
  lchar = - 1
  isgn = 1
  rtop = 0.0E+00
  rbot = 1.0E+00
  jsgn = 1
  jtop = 0
  jbot = 1
  ihave = 1
  iterm = 0

10    continue

  lchar = lchar + 1
  chrtmp = string(lchar+1:lchar+1)
!
!  Blank or TAB character.
!
  if ( chrtmp == ' ' .or. chrtmp == TAB ) then
!
!  20 November 1993
!
!  I would like to allow input like "+ 2", where there is a space
!  between the plus and the number.  So I am going to comment out
!  this line, because I think that's all that's keeping me from
!  doing this.
!
!       if ( ihave == 2 .or.
!    &     ihave == 6 .or.
!    &     ihave == 7 ) then

    if ( ihave == 2 ) then

    else if ( ihave == 6 .or. ihave == 7 ) then
      iterm = 1
    else if ( ihave > 1 ) then
      ihave = 11
    end if
!
!  Comma.
!
  else if ( chrtmp == ',' .or. chrtmp == ';' ) then

    if ( ihave /= 1 ) then
      iterm = 1
      ihave = 12
      lchar = lchar + 1
    end if
!
!  Minus sign.
!
  else if ( chrtmp == '-' ) then

    if ( ihave == 1 ) then
      ihave = 2
      isgn = - 1
    else if ( ihave == 6 ) then
      ihave = 7
      jsgn = - 1
    else
      iterm = 1
    end if
!
!  Plus sign.
!
  else if ( chrtmp == '+' ) then

    if ( ihave == 1 ) then
      ihave = 2
    else if ( ihave == 6 ) then
      ihave = 7
    else
      iterm = 1
    end if
!
!  Decimal point.
!
  else if ( chrtmp == '.' ) then

    if ( ihave < 4 ) then
      ihave = 4
    else if ( ihave >= 6 .and. ihave <= 8 ) then
      ihave = 9
    else
      iterm = 1
    end if
!
!  Exponent marker.
!
  else if ( ch_eqi ( chrtmp, 'E' ) .or. ch_eqi ( chrtmp, 'D' ) ) then

    if ( ihave < 6 ) then
      ihave = 6
    else
      iterm = 1
    end if
!
!  Digit.
!
  else if ( ihave < 11 .and. lge ( chrtmp, '0' ) .and. lle ( chrtmp, '9' ) ) then

    if ( ihave <= 2 ) then
      ihave = 3
    else if ( ihave == 4 ) then
      ihave = 5
    else if ( ihave == 6 .or. ihave == 7 ) then
      ihave = 8
    else if ( ihave == 9 ) then
      ihave = 10
    end if

    call ch_to_digit ( chrtmp, ndig )

    if ( ihave == 3 ) then
      rtop = 10.0E+00 * rtop + real ( ndig )
    else if ( ihave == 5 ) then
      rtop = 10.0E+00 * rtop + real ( ndig )
      rbot = 10.0E+00 * rbot
    else if ( ihave == 8 ) then
      jtop = 10 * jtop + ndig
    else if ( ihave == 10 ) then
      jtop = 10 * jtop + ndig
      jbot = 10 * jbot
    end if
!
!  Anything else is regarded as a terminator.
!
  else
    iterm = 1
  end if
!
!  If we haven't seen a terminator, and we haven't examined the
!  entire string, go get the next character.
!
  if ( iterm /= 1 .and. lchar+1 < nchar ) then
    go to 10
  end if
!
!  If we haven't seen a terminator, and we have examined the
!  entire string, then we're done, and LCHAR is equal to NCHAR.
!
  if ( iterm /= 1 .and. lchar+1 == nchar ) then
    lchar = nchar
  end if
!
!  Number seems to have terminated.  Have we got a legal number?
!  Not if we terminated in states 1, 2, 6 or 7!
!
  if ( ihave == 1 .or. ihave == 2 .or. ihave == 6 .or. ihave == 7 ) then

    ierror = ihave

    return
  end if
!
!  Number seems OK.  Form it.
!
  if ( jtop == 0 ) then
    rexp = 1.0E+00
  else

    if ( jbot == 1 ) then
      rexp = 10.0E+00**( jsgn * jtop )
    else
      rexp = jsgn * jtop
      rexp = rexp / jbot
      rexp = 10.0E+00**rexp
    end if

  end if

  rval = isgn * rexp * rtop / rbot

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone

  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end