program bdmlib_prb
!
!*******************************************************************************
!
!! BDMLIB_PRB tests the BDMLIB Bayes Dirichlet Mixture routines.
!
  implicit none
!
  integer, parameter :: acid_num = 20
  integer, parameter :: comp_max = 9
  integer, parameter :: iunit = 1
!
  character acid_sym(acid_num)
  real beta(acid_num,comp_max)
  real beta_sum(comp_max)
  integer comp_label(comp_max)
  integer comp_num
  real comp_weight(comp_max)
  integer ierror
  character ( len = 30 ) mixture_file_name
!
  call timestamp ( )

  write ( *, * ) ' '
  write ( *, * ) 'BDMLIB_PRB'
  write ( *, * ) '  Test BDMLIB, Bayes Dirichlet Mixture'
  write ( *, * ) '  estimate weights of a Dirichlet mixture.'
!
!  Read information about the mixture.
!
  mixture_file_name = 'mixture.dat'

  open ( unit = iunit, file = mixture_file_name, form = 'formatted' )

  call mixture_read ( acid_num, acid_sym, beta, beta_sum, &
    comp_label, comp_max, comp_num, comp_weight, ierror, iunit )

  close ( unit = iunit )
!
!  Print the amino acid parameters.
!
  write ( *, * ) ' '
  write ( *, * ) 'Numeric key for amino acid abbreviations.'
  write ( *, * ) ' '

  call amino_print ( acid_num, acid_sym )
!
!  Print the component parameters.
!
  call comp_param_print ( acid_num, acid_sym, comp_max, comp_num, &
    beta, beta_sum, comp_weight )
!
!  Check the parameters.
!
  call dirichlet_mix_check ( comp_num, acid_num, acid_num, beta, comp_weight )
!
!  Perform the simple test of generating 10 isoleucines in a row.
!
  call test01 ( acid_num, beta, comp_num, comp_weight )
!
!  Now test a random sampling.
!
  call test02 ( acid_num, beta, comp_num, comp_weight )

  write ( *, * ) ' '
  write ( *, * ) 'BDMLIB_PRB'
  write ( *, * ) '  Normal end of execution.'

  stop
end
subroutine test01 ( acid_num, beta, comp_num, comp_weight )
!
!*******************************************************************************
!
!! TEST01 generates 10 isoleucine events in a row.
!
  implicit none
!
  integer acid_num
  integer comp_num
!
  integer acid_i
  real alpha(comp_num)
  real alpha_0(comp_num)
  real beta(acid_num,comp_num)
  integer comp_i
  real comp_weight(comp_num)
  real comp_weight_est(comp_num)
  integer event_i
  integer event_num
  real p(comp_num)
  real p_hat(comp_num)
  integer site_num
  integer x_sample(acid_num)
!
  write ( *, * ) ' '
  write ( *, * ) 'TEST01'
  write ( *, * ) '  Generate a (nonrandom) sequence of'
  write ( *, * ) '  10 isoleucine results in a row.'
  write ( *, * ) ' '
!
!  Initialize information about the Bayesian process.
!
  alpha_0(1:comp_num) = 1.0E+00

  call rvec_copy ( comp_num, alpha_0, alpha )

  p(1:comp_num) = 1.0E+00 / real ( comp_num )
  p_hat(1:comp_num) = 1.0E+00 / real ( comp_num )

  event_num = 10

  call rvec_print ( comp_num, comp_weight, 'Exact component weights:' )

  call rvec_print ( comp_num, alpha, 'Initial ALPHA:' )
!
!  Based on the current ALPHA's, compute the mean/expected value/estimate 
!  for the weights.
!
  call dirichlet_mean ( comp_num, alpha, comp_weight_est )

  call rvec_print ( comp_num, comp_weight_est, &
    'Initial estimated component weights:' )

  site_num = 1

  do event_i = 1, event_num
!
!  Observe a single isoleucine.
!
    x_sample(1:acid_num) = 0
    x_sample(8) = site_num
!
!  Update ALPHA, the estimated weight parameters, based on X.
!
    call event_process ( acid_num, alpha, beta, comp_num, p, p_hat, site_num, &
      x_sample )

    call rvec_print ( comp_num, alpha, 'Current ALPHA:' )
!
!  Based on the current ALPHA's, compute the mean/expected value/estimate 
!  for the weights.
!
    call dirichlet_mean ( comp_num, alpha, comp_weight_est )

    call rvec_print ( comp_num, comp_weight_est, &
      'Estimated component weights:' )

  end do

  return
end
subroutine test02 ( acid_num, beta, comp_num, comp_weight )
!
!*******************************************************************************
!
!! TEST02 generates random events.
!
  implicit none
!
  integer acid_num
  integer comp_num
!
  integer acid_i
  real alpha(comp_num)
  real alpha_0(comp_num)
  real beta(acid_num,comp_num)
  integer comp_i
  integer comp_sample
  real comp_weight(comp_num)
  real comp_weight_est(comp_num)
  integer event_i
  integer event_num
  real p(comp_num)
  real p_hat(comp_num)
  real p_sample(acid_num)
  integer site_num
  integer x_sample(acid_num)
!
  write ( *, * ) ' '
  write ( *, * ) 'TEST02'
  write ( *, * ) '  Generate many random events.'
  write ( *, * ) '  We should be able to approximate the'
  write ( *, * ) '  exact component weights.'
  write ( *, * ) ' '
!
!  Initialize information about the Bayesian process.
!
  alpha_0(1:comp_num) = 1.0E+00

  call rvec_copy ( comp_num, alpha_0, alpha )

  p(1:comp_num) = 1.0E+00 / real ( comp_num )
  p_hat(1:comp_num) = 1.0E+00 / real ( comp_num )

  event_num = 1000

  call rvec_print ( comp_num, comp_weight, 'Exact component weights:' )

  call rvec_print ( comp_num, alpha, 'Initial ALPHA:' )
!
!  Based on the current ALPHA's, compute the mean/expected value/estimate 
!  for the weights.
!
  call dirichlet_mean ( comp_num, alpha, comp_weight_est )

  call rvec_print ( comp_num, comp_weight_est, &
    'Initial estimated component weights:' )

  site_num = 10

  do event_i = 1, event_num
!
!  Randomly choose COMP_SAMPLE, the component PDF to sample.
!
    call discrete_sample ( comp_num, comp_weight, comp_sample )
!
!  Now generate the probabilities P_SAMPLE for a multinomial PDF by sampling
!  the COMP_SAMPLE-th Dirichlet distribution.
!
    call dirichlet_sample ( acid_num, beta(1,comp_sample), p_sample )
!
!  Now generate the event X_SAMPLE by sampling the multinomial PDF with
!  the given P_SAMPLE parameters.
!
    call multinomial_sample ( site_num, acid_num, p_sample, x_sample )
!
!  Update ALPHA, the estimated weight parameters, based on X.
!
    call event_process ( acid_num, alpha, beta, comp_num, p, p_hat, site_num, &
      x_sample )
!
!  Based on the current ALPHA's, compute the mean/expected value/estimate 
!  for the weights.
!
    call dirichlet_mean ( comp_num, alpha, comp_weight_est )

    if ( event_i <= 10 .or. mod ( event_i, event_num/10 ) == 0 ) then

      write ( *, * ) ' '
      write ( *, * ) 'Event ', event_i

      call rvec_print ( comp_num, alpha, 'Current ALPHA:' )

      call rvec_print ( comp_num, comp_weight_est, &
        'Estimated component weights:' )

    end if

  end do

  return
end