# include <cmath>
# include <cstdlib>
# include <iostream>
# include <fstream>
# include <iomanip>
# include <ctime>
# include <cstring>

using namespace std;

# include "dream1.hpp"
# include "rnglib.hpp"

//****************************************************************************80

void chain_init ( int chain_num, double fit[], int gen_num, int mnor_ind[], 
  int mnor_num, int par_num, double parm[], int prior_den_par_num[], 
  double prior_den_par_val[], int prior_den_type[], 
  double sample_likelihood ( int par_num, double zp[] ),
  int uni_ind[], int uni_num, double z[] )

//****************************************************************************80
//
//  Purpose:
//
//    CHAIN_INIT starts Markov chains from a prior distribution.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    13 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, int CHAIN_NUM, the total number of chains.
//
//    Output, double FIT[CHAIN_NUM*GEN_NUM], the likelihood of
//    each sample.
//
//    Input, int GEN_NUM, the total number of generations.
//
//    Input, int MNOR_IND[MNOR_NUM], indices of parameters
//    associated with multivariate normal distributions.
//
//    Input, int MNOR_NUM, the number of parameters with
//    a multivariate normal distribution.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input, double PARM[MNOR_NUM*(MNOR_NUM+3)/2+1], parameters 
//    for a multivariate normal distribution.
//
//    Input, int PRIOR_DEN_PAR_NUM[PAR_NUM], the number of
//    variables associated with the prior density for each parameter.
//
//    Input, double PRIOR_DEN_PAR_VAL[2*PAR_NUM], parameter values
//    associated with prior densities.
//
//    Input, int PRIOR_DEN_TYPE[PAR_NUM], the type
//    of the prior density.
//
//    Input, double SAMPLE_LIKELIHOOD ( int PAR_NUM, double ZP[] ),
//    a user-defined function to evaluate the log likelihood function of
//    a sample vector ZP of dimension PAR_NUM.
//
//    Input, int UNI_IND[UNI_NUM], indices of parameters
//    associated with univariate distributions.
//
//    Input, int UNI_NUM, the number of parameters with
//    a univariate distribution.
//
//    Output, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain 
//    sample data.
//
{
  int c;
  int p;
  double *zp;

  for ( c = 0; c < chain_num; c++ )
  {
    zp = prior_sample ( mnor_ind, mnor_num, par_num, parm, 
      prior_den_par_num, prior_den_par_val, prior_den_type, uni_ind, 
      uni_num );

    for ( p = 0; p < par_num; p++ )
    {
      z[p+c*par_num+0*par_num*chain_num] = zp[p];
    }

    fit[c+0*chain_num] = sample_likelihood ( par_num, zp );

    delete [] zp;
  }
  return;
}
//****************************************************************************80

void chain_init_print ( int chain_num, double fit[], int gen_num, int par_num, 
  bool restart_flag, string restart_read_filename, double z[] )

//****************************************************************************80
//
//  Purpose:
//
//    CHAIN_INIT_PRINT prints the initial values for Markov chains.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    03 May 2013
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int CHAIN_NUM, the total number of chains.
//
//    Input, double FIT[CHAIN_NUM*GEN_NUM], the likelihood of
//    each sample.
//
//    Input, int GEN_NUM, the total number of generations.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input, bool RESTART_FLAG, is TRUE if the chains are to be initialized
//    from a restart file.
//
//    Input, string RESTART_READ_FILENAME, the name of the 
//    restart file.
//
//    Input, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain 
//    sample data.
//
{
  int i;
  int j;

  cout << "\n";
  cout << "CHAIN_INIT_PRINT\n";
  cout << "  Display initial values of Markov chains.\n";

  if ( restart_flag )
  {
    cout << "  Initialization from restart file \"" 
         << restart_read_filename << "\"\n";
  }
  else
  {
    cout << "  Initialization by sampling prior density.\n";
  }
  for ( j = 0; j < chain_num; j++ )
  {
    cout << "\n";
    cout << "  Chain " << j << "\n";
    cout << "  Fitness " << fit[j+0*chain_num] << "\n";
    for ( i = 0; i < par_num; i++ )
    {
      cout << "  " << setw(14) << z[i+j*par_num+0*par_num*chain_num];
      if ( ( i + 1 ) % 5 == 0 || i == par_num - 1 )
      {
        cout << "\n";
      }
    }
  }
  return;
}
//****************************************************************************80

void chain_outliers ( int chain_num, int gen_index, int gen_num, int par_num,
  double fit[], double z[] )

//****************************************************************************80
//
//  Purpose:
//
//    CHAIN_OUTLIERS identifies and modifies outlier chains during burn-in.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Jasper Vrugt, CJF ter Braak, CGH Diks, Bruce Robinson, James Hyman, 
//    Dave Higdon,
//    Accelerating Markov Chain Monte Carlo Simulation by Differential 
//    Evolution with Self-Adaptive Randomized Subspace Sampling,
//    International Journal of Nonlinear Sciences and Numerical Simulation,
//    Volume 10, Number 3, March 2009, pages 271-288.
//
//  Parameters:
//
//    Input, int CHAIN_NUM, the total number of chains.
//
//    Input, int GEN_INDEX, the index of the current generation.
//    2 <= GEN_INDEX <= GEN_NUM.
//
//    Input, int GEN_NUM, the total number of generations.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input/output, double FIT[CHAIN_NUM*GEN_NUM], the likelihood of
//    each sample.
//
//    Input/output, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov
//    chain sample data.
//
{
  double *avg;
  double avg_max;
  double *avg_sorted;
  int best;
  int i;
  int ind1;
  int ind3;
  int j;
  int klo;
  int knum;
  int k;
  int outlier_num;
  double q1;
  double q3;
  double qr;
  double t;

  klo = ( ( gen_index + 1 ) / 2 ) - 1;
  knum = gen_index + 1 - klo;

  avg = new double[chain_num];

  for ( j = 0; j < chain_num; j++ )
  {
    t = 0.0;
    for ( k = klo; k <= gen_index; k++ )
    {
      t = t + fit[j+k*chain_num];
    }
    avg[j] = t  / ( double ) ( knum );
  }
//
//  Set BEST to be the index of the chain with maximum average.
//
  best = 0;
  avg_max = avg[0];
  for ( j = 1; j < chain_num; j++ )
  {
    if ( avg_max < avg[j] )
    {
      best = j;
      avg_max = avg[j];
    }
  }
//
//  Determine the indices of the chains having averages 1/4 "above" 
//  and "below" the average.
//
  avg_sorted = r8vec_copy_new ( chain_num, avg );

  r8vec_sort_heap_a ( chain_num, avg_sorted );

  ind1 = r8_round_i4 ( 0.25 * ( double ) ( chain_num ) );
  ind3 = r8_round_i4 ( 0.75 * ( double ) ( chain_num ) );

  q1 = avg_sorted[ind1];
  q3 = avg_sorted[ind3];
  qr = q3 - q1;

  delete [] avg_sorted;
//
//  Identify outlier chains, and replace their later samples
//  with values from the "best" chain.
//
  outlier_num = 0;
  for ( j = 0; j < chain_num; j++ )
  {
    if ( avg[j] < q1 - 2.0 * qr )
    {
      outlier_num = outlier_num + 1;
      for ( i = 0; i < par_num; i++ )
      {
        z[i+j*par_num+gen_index*par_num*chain_num] = 
          z[i+best*par_num+gen_index*par_num*chain_num];
      }
      for ( k = klo; k <= gen_index; k++ )
      {
        fit[j+k*chain_num]  = fit[best+k*chain_num];
      }
    }
  }
//
//  List the outlier chains.
//
  if ( 0 < outlier_num )
  {
    cout << "\n";
    cout << "CHAIN_OUTLIERS:\n";
    cout << "  At iteration " << gen_index 
         << " found " << outlier_num << " outlier chains,\n";
    cout << "  whose indices appear below, and for which samples\n";
    cout << "  from the chain with the largest log likelihood function,\n";
    cout << "  index number " << best << " will be substituted.\n";

    for ( j = 0; j < chain_num; j++ )
    {
      if ( avg[j] < q1 - 2.0 * qr )
      {
        cout << "  " << j << "\n";
      }
    }
  }

  delete [] avg;

  return;
}
//****************************************************************************80

void chain_write ( string chain_filename, int chain_num, double fit[], 
  int gen_num, int par_num, double z[] )

//****************************************************************************80
//
//  Purpose:
//
//    CHAIN_WRITE writes samples of each chain to separate files.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version John Burkardt.
//
//  Parameters:
//
//    Input, string CHAIN_FILENAME, the base name for
//    the files.  It should include a string of 0's that will be replaced
//    by the indices of each chain.
//
//    Input, int CHAIN_NUM, the total number of chains.
//
//    Input, double FIT[CHAIN_NUM*GEN_NUM], the likelihood of
//    each sample.
//
//    Input, int GEN_NUM, the total number of generations.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain 
//    sample data.
//
{
  ofstream chain;
  string chain_filename2;
  int i;
  int j;
  int k;
//
//  Make a temporary copy of the filename template, which we can alter.
//
  chain_filename2 = chain_filename;
//
//  Write parameter samples of all chains.
//
  cout << "\n";
  cout << "CHAIN_WRITE:\n";

  for ( j = 0; j < chain_num; j++ )
  {
    chain.open ( chain_filename2.c_str ( ) );

    if ( ! chain )
    {
      cerr << "\n";
      cerr << "CHAIN_WRITE - Fatal error!\n";
      cerr << "  Could not open file \"" << chain_filename2 << "\".\n";
      exit ( 1 );
    }

    chain <<
      "DREAM.CPP:Parameters_and_log_likelihood_for_chain_#" << j << "\n";

    for ( k = 0; k < gen_num; k++ )
    {
      chain << "  " << k
            << "  " << fit[j+k*chain_num];
      for ( i = 0; i < par_num; i++ )
      {
        chain << "  " << z[i+j*par_num+k*par_num*chain_num];
      }
      chain << "\n";
    }

    chain.close ( );

    cout << "  Created file \"" << chain_filename2 << "\".\n";

    filename_inc ( &chain_filename2 );
  }
  return;
}
//****************************************************************************80

void cr_dis_update ( int chain_index, int chain_num, double cr_dis[], 
  int cr_index, int cr_num, int cr_ups[], int gen_index, int gen_num, 
  int par_num, double z[] )

//****************************************************************************80
//
//  Purpose:
//
//    CR_DIS_UPDATE updates the CR distance.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    29 August 2018
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, int CHAIN_INDEX, the index of the chain.
//    0 <= CHAIN_INDEX < CHAIN_NUM.
//
//    Input, int CHAIN_NUM, the total number of chains.
//
//    Input/output, double CR_DIS[CR_NUM], the CR distances.
//
//    Input, int CR_INDEX, the index of the CR.
//    0 <= CR_INDEX < CR_NUM.
//
//    Input, int CR_NUM, the total number of CR values.
//
//    Input/output, int CR_UPS[CR_NUM], the number of updates 
//    for each CR.
//
//    Input, int GEN_INDEX, the index of the generation.
//    0 <= GEN_INDEX < GEN_NUM.
//
//    Input, int GEN_NUM, the total number of generations.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain 
//    sample data.
//
{
  int i;
  int i1;
  int i2;
  double *std;
//
//  Compute the standard deviations.
//
  std = std_compute ( chain_num, gen_index, gen_num, par_num, z );
//
//  Increment the update count.
//
  cr_ups[cr_index] = cr_ups[cr_index] + 1;
//
//  Update the CR distance.
//
  for ( i = 0; i < par_num; i++ )
  {
    i1 = i + chain_index * par_num +   gen_index       * par_num * chain_num;
    i2 = i + chain_index * par_num + ( gen_index - 1 ) * par_num * chain_num;
    cr_dis[cr_index] = cr_dis[cr_index] + pow ( ( z[i2] - z[i1] ) / std[i], 2 );
  }

  delete [] std;

  return;
}
//****************************************************************************80

int cr_index_choose ( int cr_num, double cr_prob[] )

//****************************************************************************80
//
//  Purpose:
//
//    CR_INDEX_CHOOSE chooses a CR index.
//
//  Discussion:
//
//    Index I is chosen with probability CR_PROB(I).
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, int CR_NUM, the total number of CR values.
//
//    Input, double CR_PROB[CR_NUM], the probability of each CR.
//
//    Output, int CR_INDEX_CHOOSE, the index of the CR.
//    0 <= CR_INDEX_CHOOSE < CR_NUM.
//
{
  int cr_index;
  int i;
  int n;
  int *tmp_index;

  if ( cr_num == 1 )
  {
    cr_index = 0;
  }
  else
  { 
    n = 1;
    tmp_index = i4vec_multinomial_sample ( n, cr_prob, cr_num );

    for ( i = 0; i < cr_num; i++ )
    {
      if ( tmp_index[i] == 1 )
      {
        cr_index = i;
        break;
      }
    }
    delete [] tmp_index;
  }
  return cr_index;
}
//****************************************************************************80

void cr_init ( double cr[], double cr_dis[], int cr_num, double cr_prob[], 
  int cr_ups[] )

//****************************************************************************80
//
//  Purpose:
//
//    CR_INIT initializes the crossover probability values.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Output, double CR[CR_NUM], the CR values.
//
//    Output, double CR_DIS[CR_NUM], the CR distances.
//
//    Input, int CR_NUM, the total number of CR values.
//
//    Output, double CR_PROB[CR_NUM], the probability of each CR.
//
//    Output, int CR_UPS[CR_NUM], the number of updates
//    for each CR.
//
{
  int i;

  for ( i = 0; i < cr_num; i++ )
  {
    cr[i] = ( double ) ( i + 1 ) / ( double ) ( cr_num );
    cr_dis[i] = 1.0;
    cr_prob[i] = 1.0 / ( double ) ( cr_num );
    cr_ups[i] = 1;
  }
  return;
}
//****************************************************************************80

void cr_prob_update ( double cr_dis[], int cr_num, double cr_prob[], 
  int cr_ups[] )

//****************************************************************************80
//
//  Purpose:
//
//    CR_PROB_UPDATE updates the CR probabilities.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double CR_DIS[CR_NUM], the CR distances.
//
//    Input, int CR_NUM, the total number of CR values.
//
//    Output, double CR_PROB[CR_NUM], the updated CR probabilities.
//
//    Input, int CR_UPS[CR_NUM], the number of updates 
//    for each CR.
//
{
  double cr_prob_sum;
  int i;
 
  for ( i = 0; i < cr_num - 1; i++ )
  {
    cr_prob[i] = cr_dis[i] / ( double ) cr_ups[i];
  }

  cr_prob_sum = r8vec_sum ( cr_num, cr_prob );

  for ( i = 0; i < cr_num - 1; i++ )
  {
    cr_prob[i] = cr_prob[i] / cr_prob_sum;
  }

  return;
}
//****************************************************************************80

double *diff_compute ( int chain_num, int gen_index, int gen_num, 
  int jump_dim[], int jump_num, int pair_num, int par_num, int r[], 
  double z[] ) 

//****************************************************************************80
//
//  Purpose:
//
//    DIFF_COMPUTE computes the differential evolution.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Jasper Vrugt, CJF ter Braak, CGH Diks, Bruce Robinson, James Hyman, 
//    Dave Higdon,
//    Accelerating Markov Chain Monte Carlo Simulation by Differential 
//    Evolution with Self-Adaptive Randomized Subspace Sampling,
//    International Journal of Nonlinear Sciences and Numerical Simulation,
//    Volume 10, Number 3, March 2009, pages 271-288.
//
//  Parameters:
//
//    Input, int CHAIN_NUM, the total number of chains.
//
//    Input, int GEN_INDEX, the index of the current generation.
//    1 <= GEN_INDEX <= GEN_NUM.
//
//    Input, int GEN_NUM, the total number of generations.
//
//    Input, int JUMP_DIM[JUMP_NUM], the dimensions in which
//    a jump is to be made.
//
//    Input, int JUMP_NUM, the number of dimensions in which
//    a jump will be made.  0 <= JUMP_NUM <= PAR_NUM.
//
//    Input, int PAIR_NUM, the number of pairs of 
//    crossover chains.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input, int R[2*PAIR_NUM], pairs of chains used
//    to compute differences.
//
//    Input, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain 
//    sample data.
//
//    Output, double DIFF_COMPUTE[PAR_NUM], the vector of pair differences.
//
{
  double *diff;
  int i1;
  int i2;
  int j;
  int k;
  int pair;
  int r1;
  int r2;
//
//  Produce the difference of the pairs used for population evolution.
//
  diff = r8vec_zero_new ( par_num );

  for ( pair = 0; pair < pair_num; pair++ )
  {
    r1 = r[0+pair*2];
    r2 = r[1+pair*2];
    for ( j = 0; j < jump_num; j++ )
    {
      k = jump_dim[j];
      i1 = k+r1*par_num+(gen_index-1)*par_num*chain_num;
      i2 = k+r2*par_num+(gen_index-1)*par_num*chain_num;
      diff[k] = diff[k] + ( z[i1] - z[i2] );
    }
  }

  return diff;
}
//****************************************************************************80

void dream_algm ( int chain_num, double cov_mat[], int cr_num, double fit[], 
  int gen_num, double gr[], int &gr_conv, int &gr_count, int gr_num, 
  double gr_threshold, double jumprate_table[], int jumpstep, double limits[], 
  double logdet, int mnor_ind[], double mnor_mean[], int mnor_num, 
  int pair_num, int par_num, int printstep, int prior_den_par_num[], 
  double prior_den_par_val[], int prior_den_type[], 
  double sample_likelihood ( int par_num, double zp[] ),
  int uni_ind[], int uni_num, double z[] )
    
//****************************************************************************80
//
//  Purpose:
//
//    DREAM_ALGM gets a candidate parameter sample.
//
// Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    13 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Jasper Vrugt, CJF ter Braak, CGH Diks, Bruce Robinson, James Hyman, 
//    Dave Higdon,
//    Accelerating Markov Chain Monte Carlo Simulation by Differential 
//    Evolution with Self-Adaptive Randomized Subspace Sampling,
//    International Journal of Nonlinear Sciences and Numerical Simulation,
//    Volume 10, Number 3, March 2009, pages 271-288.
//
//  Parameters:
//
//    Input, int CHAIN_NUM, the total number of chains.
//
//    Input, double COV_MAT[MNOR_NUM*MNOR_NUM], the covariance
//    matrix for a multivariate normal density.
//
//    Input, int CR_NUM, the total number of CR values.
//
//    Input, double FIT[CHAIN_NUM*GEN_NUM], the likelihood of
//    each sample.
//
//    Input, int GEN_NUM, the total number of generations.
//
//    Input, double GR[PAR_NUM*GR_NUM], 
//    the Gelman-Rubin R statistic.
//
//    Input/output, int &GR_CONV, the Gelman-Rubin convergence flag.
//
//    Input/output, int &GR_COUNT, counts the number of generations
//    at which the Gelman-Rubin statistic has been computed.
//
//    Input, int GR_NUM, the number of times the Gelman-Rubin
//    statistic may be computed.
//
//    Input, double GR_THRESHOLD, the convergence tolerance for
//    the Gelman-Rubin statistic.
//
//    Input, double JUMPRATE_TABLE[PAR_NUM], the jumprate table.
//
//    Input, int JUMPSTEP, forces a "long jump" every
//    JUMPSTEP generations.
//
//    Input, double LIMITS[2*PAR_NUM], lower and upper bounds
//    for each parameter.
//
//    Input, double LOGDET, the logarithm of the determinant
//    of the covariance matrix.
//
//    Input, int MNOR_IND[MNOR_NUM], the indices of
//    parameters with multivariate density.
//
//    Input, int MNOR_MEAN[MNOR_NUM], the means associated
//    with the multivariate density.
//
//    Input, int MNOR_NUM, the number of parameters with
//    a multivariate normal distribution.
//
//    Input, int PAIR_NUM, the number of pairs of 
//    crossover chains.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input, int PRINTSTEP, the interval between generations on 
//    which the Gelman-Rubin statistic will be computed and written to a file.
//
//    Input, int PRIOR_DEN_PAR_NUM[PAR_NUM], the number of
//    variables associated with the prior density for each parameter.
//
//    Input, double PRIOR_DEN_PAR_VAL[2*PAR_NUM], parameter values
//    associated with prior densities.
//
//    Input, int PRIOR_DEN_TYPE[PAR_NUM], the type
//    of the prior density.
//
//    Input, double SAMPLE_LIKELIHOOD ( int PAR_NUM, double ZP[] ),
//    a user-defined function to evaluate the log likelihood function of
//    a sample vector ZP of dimension PAR_NUM.
//
//    Input, int UNI_IND[UNI_NUM], indices of parameters 
//    associated with univariate distributions.
//
//    Input, int UNI_NUM, the number of parameters with
//    a univariate distribution.
//
//    Output, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain 
//    sample data.
//
//  Local parameters:
//
//    Local, int CHAIN_INDEX, the index of the current chain.
//    1 <= CHAIN_INDEX <= CHAIN_NUM.
//
//    Local, double CR[CR_NUM], the CR values.
//
//    Local, double CR_DIS[CR_NUM], the CR distances.
//
//    Local, int CR_INDEX, the index of the selected CR value.
//    1 <= CR_INDEX <= CR_NUM.
//
//    Local, double CR_PROB[CR_NUM], the probability of each CR.
//
//    Local, double CR_UPS[CR_NUM], the number of updates for each CR.
//
//    Local, int GEN_INDEX, the index of the current generation.
//    1 <= GEN_INDEX <= GEN_NUM.
//
//    Local, double ZP[PAR_NUM], a candidate sample.
//
//    Local, int ZP_ACCEPT, the number of candidates accepted.
//
//    Local, double ZP_ACCEPT_RATE, the rate at which generated
//    candidates were accepted.
//
//    Local, int ZP_COUNT, the number of candidates generated.
//
//    Local, double ZP_RATIO, the Metropolis ratio for a candidate.
//
{
  int chain_index;
  double *cr;
  double *cr_dis;
  int cr_index;
  double *cr_prob;
  int *cr_ups;
  int gen_index;
  int i;
  double pd1;
  double pd2;
  double r;
  double *zp;
  int zp_accept;
  double zp_accept_rate;
  int zp_count;
  double zp_fit;
  double *zp_old;
  double zp_old_fit;
  double zp_ratio;

  zp_old = new double[par_num];
  zp_count = 0;
  zp_accept = 0;
//
//  Initialize the CR values.
//
  cr = new double[cr_num];
  cr_dis = new double[cr_num];
  cr_prob = new double[cr_num];
  cr_ups = new int[cr_num];

  cr_init ( cr, cr_dis, cr_num, cr_prob, cr_ups );

  for ( gen_index = 1; gen_index < gen_num; gen_index++ )
  {
    for ( chain_index = 0; chain_index < chain_num; chain_index++ )
    {
//
//  Choose CR_INDEX, the index of a CR.
//
      cr_index = cr_index_choose ( cr_num, cr_prob );
//
//  Generate a sample candidate ZP.
//
      zp = sample_candidate ( chain_index, chain_num, cr, cr_index, cr_num, 
        gen_index, gen_num, jumprate_table, jumpstep, limits, pair_num, 
        par_num, z );

      zp_count = zp_count + 1;
//
//  Compute the log likelihood function for ZP.
//
      zp_fit = sample_likelihood ( par_num, zp );

      for ( i = 0; i < par_num; i++ )
      {
        zp_old[i] = z[i+chain_index*par_num+(gen_index-1)*par_num*chain_num];
      }
      zp_old_fit = fit[chain_index+(gen_index-1)*chain_num];
//
//  Compute the Metropolis ratio for ZP.
//
      pd1 = prior_density ( cov_mat, logdet, mnor_ind, mnor_mean, mnor_num, 
        par_num, prior_den_par_num, prior_den_par_val, prior_den_type, 
        uni_ind, uni_num, zp );

      pd2 = prior_density ( cov_mat, logdet, mnor_ind, mnor_mean, mnor_num, 
        par_num, prior_den_par_num, prior_den_par_val, prior_den_type, 
        uni_ind, uni_num, 
        z+0+chain_index*par_num+(gen_index-1)*par_num*chain_num );

      zp_ratio = exp ( 
        ( zp_fit     + log ( pd1 ) ) - 
        ( zp_old_fit + log ( pd2 ) ) );

      zp_ratio = r8_min ( zp_ratio, 1.0 );
//
//  Accept the candidate, or copy the value from the previous generation.
//
      r = r8_uniform_01_sample ( );

      if ( r <= zp_ratio )
      {
        for ( i = 0; i < par_num; i++ )
        {
          z[i+chain_index*par_num+gen_index*par_num*chain_num] = zp[i];
        }
        zp_accept = zp_accept + 1;
        fit[chain_index+gen_index*chain_num] = zp_fit;
      }
      else
      {
        for ( i = 0; i < par_num; i++ )
        {
          z[i+chain_index*par_num+gen_index*par_num*chain_num] = zp_old[i];
        }
        fit[chain_index+gen_index*chain_num] = zp_old_fit; 
      }
//
//  Update the CR distance.
//
      if ( ! gr_conv )
      {
        if ( 1 < cr_num )
        {
          cr_dis_update ( chain_index, chain_num, cr_dis, cr_index, 
            cr_num, cr_ups, gen_index, gen_num, par_num, z );
        }
      }

      delete [] zp;
    }
//
//  Update the multinomial distribution of CR.
//
    if ( ! gr_conv )
    {
      if ( 1 < cr_num )
      {
        if ( ( gen_index + 1 ) % 10 == 0 )
        {
          cr_prob_update ( cr_dis, cr_num, cr_prob, cr_ups );
        }
      }
    }
//
//  Every PRINTSTEP interval,
//  * compute the Gelman Rubin R statistic for this generation,
//    and determine if convergence has occurred.
//
    if ( ( gen_index + 1 ) % printstep == 0 )
    {
      gr_compute ( chain_num, gen_index, gen_num, gr, gr_conv, gr_count, 
        gr_num, gr_threshold, par_num, z );
    }
//
//  Check for outlier chains.
//
    if ( ! gr_conv )
    {
      if ( ( gen_index + 1 ) % 10 == 0 )
      {
        chain_outliers ( chain_num, gen_index, gen_num, par_num, fit, z );
      }
    }
  }    
//
//  Compute the acceptance rate.
//
  zp_accept_rate = ( double ) ( zp_accept ) / ( double ) ( zp_count );

  cout << "\n";
  cout << "  The acceptance rate is " << zp_accept_rate << "\n";

  delete [] cr;
  delete [] cr_dis;
  delete [] cr_prob;
  delete [] cr_ups;
  delete [] zp_old;

  return;
}
//****************************************************************************80

void filename_inc ( string *filename )

//****************************************************************************80
//
//  Purpose:
//
//    FILENAME_INC increments a partially numeric file name.
//
//  Discussion:
//
//    It is assumed that the digits in the name, whether scattered or
//    connected, represent a number that is to be increased by 1 on
//    each call.  If this number is all 9's on input, the output number
//    is all 0's.  Non-numeric letters of the name are unaffected.
//
//    If the name is empty, then the routine stops.
//
//    If the name contains no digits, the empty string is returned.
//
//  Example:
//
//      Input            Output
//      -----            ------
//      "a7to11.txt"     "a7to12.txt"  (typical case.  Last digit incremented)
//      "a7to99.txt"     "a8to00.txt"  (last digit incremented, with carry.)
//      "a9to99.txt"     "a0to00.txt"  (wrap around)
//      "cat.txt"        " "           (no digits to increment)
//      " "              STOP!         (error)
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 November 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input/output, string *FILENAME, the filename to be incremented.
//
{
  char c;
  int change;
  int i;
  int lens;

  lens = (*filename).length ( );

  if ( lens <= 0 )
  {
    cerr << "\n";
    cerr << "FILENAME_INC - Fatal error!\n";
    cerr << "  The input string is empty.\n";
    exit ( 1 );
  }

  change = 0;

  for ( i = lens - 1; 0 <= i; i-- )
  {
    c = (*filename)[i];

    if ( '0' <= c && c <= '9' )
    {
      change = change + 1;

      if ( c == '9' )
      {
        c = '0';
        (*filename)[i] = c;
      }
      else
      {
        c = c + 1;
        (*filename)[i] = c;
        return;
      }
    }
  }
//
//  No digits were found.  Return blank.
//
  if ( change == 0 )
  {
    for ( i = lens - 1; 0 <= i; i-- )
    {
      (*filename)[i] = ' ';
    }
  }

  return;
}
//****************************************************************************80

void gr_compute ( int chain_num, int gen_index, int gen_num, double gr[], 
  int &gr_conv, int &gr_count, int gr_num, double gr_threshold, int par_num, 
  double z[] )

//****************************************************************************80
//
//  Purpose:
//
//    GR_COMPUTE computes the Gelman Rubin statistics R used to check
//    convergence.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Jasper Vrugt, CJF ter Braak, CGH Diks, Bruce Robinson, James Hyman, 
//    Dave Higdon,
//    Accelerating Markov Chain Monte Carlo Simulation by Differential 
//    Evolution with Self-Adaptive Randomized Subspace Sampling,
//    International Journal of Nonlinear Sciences and Numerical Simulation,
//    Volume 10, Number 3, March 2009, pages 271-288.
//
//  Parameters:
//
//    Input, int CHAIN_NUM, the total number of chains.
//
//    Input, int GEN_INDEX, the index of the current generation.
//    0 < GEN_INDEX < GEN_NUM.
//
//    Input, int GEN_NUM, the total number of generations.
//
//    Output, double GR[PAR_NUM*GR_NUM], the Gelman-Rubin R statistic.
//
//    Output, int &GR_CONV, the Gelman-Rubin convergence flag.
//
//    Input/output, int &GR_COUNT, counts the number of 
//    generations at which the Gelman-Rubin statistic has been computed.
//
//    Input, int GR_NUM, the number of times the Gelman-Rubin
//    statistic may be computed.
//
//    Input, double GR_THRESHOLD, the convergence tolerance for the
//    Gelman-Rubin statistic.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain 
//    sample data.
//
{
  double b_var;
  int chain_index;
  int ind0;
  int k;
  double mean_all;
  double *mean_chain;
  int par_index;
  double rnd0;
  double s;
  double s_sum;
  double var;
  double w_var;

  ind0 = ( ( gen_index + 1 ) / 2 ) - 1;
  rnd0 = ( double ) ( ind0 + 1 );

  mean_chain = new double[chain_num];

  for ( par_index = 0; par_index < par_num; par_index++ )
  {
    for ( chain_index = 0; chain_index < chain_num; chain_index++ )
    {
      mean_chain[chain_index] = 0.0;
      for ( k = ind0; k <= gen_index; k++ )
      {
        mean_chain[chain_index] = mean_chain[chain_index] 
          + z[par_index+chain_index*par_num+k*par_num*chain_num];
      }
      mean_chain[chain_index] = mean_chain[chain_index] / rnd0;
    }

    mean_all = r8vec_sum ( chain_num, mean_chain ) / ( double ) chain_num;

    b_var = 0.0;
    for ( chain_index = 0; chain_index < chain_num; chain_index++ )
    {
      b_var = b_var + pow ( mean_chain[chain_index] - mean_all, 2 );
    }
    b_var = rnd0 * b_var / ( double ) ( chain_num - 1 );

    s_sum = 0.0;
    for ( chain_index = 0; chain_index < chain_num; chain_index++ )
    {
      s = 0.0;
      for ( k = ind0; k <= gen_index; k++ )
      {
        s = s + pow ( z[par_index+chain_index*par_num+k*par_num*chain_num] 
          - mean_chain[chain_index], 2 );
      }
      s_sum = s_sum + s;
    }
    s_sum = s_sum / ( rnd0 - 1.0 );

    w_var = s_sum / ( double ) ( chain_num );

    var = ( ( rnd0 - 1.0 ) * w_var + b_var ) / rnd0;

    gr[par_index+gr_count*par_num] = sqrt ( var / w_var );
  }
//
//  Set the convergence flag.
//
  gr_conv = 1;

  for ( par_index = 0; par_index < par_num; par_index++ )
  {
    if ( gr_threshold < gr[par_index+gr_count*par_num] )
    {
      gr_conv = 0;
      break;
    }
  }

  if ( gr_conv ) 
  {
    cout << "\n";
    cout << "GR_COMPUTE:\n";
    cout << "  GR convergence at iteration: " << gen_index << "\n";
  }

  delete [] mean_chain;

  gr_count = gr_count + 1;

  return;
}
//****************************************************************************80

void gr_init ( double gr[], int &gr_conv, int &gr_count, int gr_num, 
  int par_num )

//****************************************************************************80
//
//  Purpose:
//
//    GR_INIT initializes Gelman-Rubin variables.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Output, double GR[PAR_NUM*GR_NUM], the Gelman-Rubin statistic.
//
//    Output, int &GR_CONV, the convergence flag.
//
//    Output, int &GR_COUNT, counts the number of generations
//    at which the Gelman-Rubin statistic has been computed.
//
//    Input, int GR_NUM, the number of times the Gelman-Rubin
//    statistic may be computed.
//
//    Input, int PAR_NUM, the number of parameters.
//
{
  int i;
  int j;

  for ( j = 0; j < gr_num; j++ )
  {
    for ( i = 0; i < par_num; i++ )
    {
      gr[i+j*par_num] = 0.0;
    }
  }
  gr_conv = 0;
  gr_count = 0;

  return;
}
//****************************************************************************80

void gr_write ( double gr[], string gr_filename, int gr_num, int par_num, 
  int printstep )

//****************************************************************************80
//
//  Purpose:
//
//    GR_WRITE writes Gelman-Rubin R statistics into a file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double GR[PAR_NUM*GR_NUM], the Gelman-Rubin R statistic.
//
//    Input, string GR_FILENAME, the Gelman-Rubin filename.
//
//    Input, int GR_NUM, the number of times the Gelman-Rubin
//    statistic may be computed.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input, int PRINTSTEP, the interval between generations on 
//    which the Gelman-Rubin statistic will be computed and written to a file.
//
{
  ofstream gr_unit;
  int i;
  int j;

  gr_unit.open ( gr_filename.c_str ( ) );

  if ( ! gr_unit )
  {
    cerr << "\n";
    cerr << "GR_WRITE - Fatal error!\n";
    cerr << "  Could not open the file \"" << gr_filename << "\"\n";
    exit ( 1 );
  }

  gr_unit <<
    "DREAM.CPP:Monitored_parameter_interchains_Gelman_Rubin_R_statistic\n";

  for ( j = 0; j < gr_num; j++ )
  {
    gr_unit << printstep * ( j + 1 ) - 1;
    for ( i = 0; i < par_num; i++ )
    {
      gr_unit << "  " << gr[i+j*par_num];
    }
    gr_unit << "\n";
  }

  gr_unit.close ( );

  cout << "\n";
  cout << "GR_WRITE:\n";
  cout << "  Created the file \"" << gr_filename << "\".\n";

  return;
}
//****************************************************************************80

int i4_binomial_sample ( int n, double pp )

//****************************************************************************80
//
//  Purpose:
//
//    I4_BINOMIAL_SAMPLE generates a binomial random deviate.
//
//  Discussion:
//
//    This procedure generates a single random deviate from a binomial
//    distribution whose number of trials is N and whose
//    probability of an event in each trial is P.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN77 version by Barry Brown, James Lovato.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Voratas Kachitvichyanukul, Bruce Schmeiser,
//    Binomial Random Variate Generation,
//    Communications of the ACM,
//    Volume 31, Number 2, February 1988, pages 216-222.
//
//  Parameters:
//
//    Input, int N, the number of binomial trials, from which a
//    random deviate will be generated.
//    0 < N.
//
//    Input, double PP, the probability of an event in each trial of
//    the binomial distribution from which a random deviate is to be generated.
//    0.0 < PP < 1.0.
//
//    Output, int I4_BINOMIAL_SAMPLE, a random deviate from the
//    distribution.
//
{
  double al;
  double alv;
  double amaxp;
  double c;
  double f;
  double f1;
  double f2;
  double ffm;
  double fm;
  double g;
  int i;
  int ix;
  int ix1;
  int k;
  int m;
  int mp;
  double p;
  double p1;
  double p2;
  double p3;
  double p4;
  double q;
  double qn;
  double r;
  double t;
  double u;
  double v;
  int value;
  double w;
  double w2;
  double x;
  double x1;
  double x2;
  double xl;
  double xll;
  double xlr;
  double xm;
  double xnp;
  double xnpq;
  double xr;
  double ynorm;
  double z;
  double z2;

  if ( pp <= 0.0 || 1.0 <= pp )
  {
    cerr << "\n";
    cerr << "I4_BINOMIAL_SAMPLE - Fatal error!\n";
    cerr << "  PP is out of range.\n";
    exit ( 1 );
  }

  p = r8_min ( pp, 1.0 - pp );
  q = 1.0 - p;
  xnp = ( double ) ( n ) * p;

  if ( xnp < 30.0 )
  {
    qn = pow ( q, n );
    r = p / q;
    g = r * ( double ) ( n + 1 );

    for ( ; ; )
    {
      ix = 0;
      f = qn;
      u = r8_uniform_01_sample ( );

      for ( ; ; )
      {
        if ( u < f )
        {
          if ( 0.5 < pp )
          {
            ix = n - ix;
          }
          value = ix;
          return value;
        }

        if ( 110 < ix )
        {
          break;
        }
        u = u - f;
        ix = ix + 1;
        f = f * ( g / ( double ) ( ix ) - r );
      }
    }
  }
  ffm = xnp + p;
  m = ffm;
  fm = m;
  xnpq = xnp * q;
  p1 = ( int ) ( 2.195 * sqrt ( xnpq ) - 4.6 * q ) + 0.5;
  xm = fm + 0.5;
  xl = xm - p1;
  xr = xm + p1;
  c = 0.134 + 20.5 / ( 15.3 + fm );
  al = ( ffm - xl ) / ( ffm - xl * p );
  xll = al * ( 1.0 + 0.5 * al );
  al = ( xr - ffm ) / ( xr * q );
  xlr = al * ( 1.0 + 0.5 * al );
  p2 = p1 * ( 1.0 + c + c );
  p3 = p2 + c / xll;
  p4 = p3 + c / xlr;
//
//  Generate a variate.
//
  for ( ; ; )
  {
    u = r8_uniform_01_sample ( ) * p4;
    v = r8_uniform_01_sample ( );
//
//  Triangle
//
    if ( u < p1 )
    {
      ix = xm - p1 * v + u;
      if ( 0.5 < pp ) 
      {
        ix = n - ix;
      }
      value = ix;
      return value;
    }
//
//  Parallelogram
//
    if ( u <= p2 )
    {
      x = xl + ( u - p1 ) / c;
      v = v * c + 1.0 - fabs ( xm - x ) / p1;

      if ( v <= 0.0 || 1.0 < v )
      {
        continue;
      }
      ix = x;
    }
    else if ( u <= p3 )
    {
      ix = xl + log ( v ) / xll;
      if ( ix < 0 )
      {
        continue;
      }
      v = v * ( u - p2 ) * xll;
    }
    else
    {
      ix = xr - log ( v ) / xlr;
      if ( n < ix )
      {
        continue;
      }
      v = v * ( u - p3 ) * xlr;
    }
    k = abs ( ix - m );

    if ( k <= 20 || xnpq / 2.0 - 1.0 <= k )
    {
      f = 1.0;
      r = p / q;
      g = ( n + 1 ) * r;

      if ( m < ix )
      {
        mp = m + 1;
        for ( i = mp; i <= ix; i++ )
        {
          f = f * ( g / i - r );
        }
      }
      else if ( ix < m )
      {
        ix1 = ix + 1;
        for ( i = ix1; i <= m; i++ )
        {
          f = f / ( g / i - r );
        }
      }

      if ( v <= f )
      {
        if ( 0.5 < pp )
        {
          ix = n - ix;
        }
        value = ix;
        return value;
      }
    }
    else
    {
      amaxp = ( k / xnpq ) * ( ( k * ( k / 3.0 
        + 0.625 ) + 0.1666666666666 ) / xnpq + 0.5 );
      ynorm = - ( double ) ( k * k ) / ( 2.0 * xnpq );
      alv = log ( v );

      if ( alv < ynorm - amaxp )
      {
        if ( 0.5 < pp )
        {
          ix = n - ix;
        }
        value = ix;
        return value;
      }

      if ( ynorm + amaxp < alv )
      {
        continue;
      }

      x1 = ( double ) ( ix + 1 );
      f1 = fm + 1.0;
      z = ( double ) ( n + 1 ) - fm;
      w = ( double ) ( n - ix + 1 );
      z2 = z * z;
      x2 = x1 * x1;
      f2 = f1 * f1;
      w2 = w * w;

      t = xm * log ( f1 / x1 ) + ( n - m + 0.5 ) * log ( z / w ) 
        + ( double ) ( ix - m ) * log ( w * p / ( x1 * q )) 
        + ( 13860.0 - ( 462.0 - ( 132.0 - ( 99.0 - 140.0 
        / f2 ) / f2 ) / f2 ) / f2 ) / f1 / 166320.0 
        + ( 13860.0 - ( 462.0 - ( 132.0 - ( 99.0 - 140.0 
        / z2 ) / z2 ) / z2 ) / z2 ) / z / 166320.0 
        + ( 13860.0 - ( 462.0 - ( 132.0 - ( 99.0 - 140.0 
        / x2 ) / x2 ) / x2 ) / x2 ) / x1 / 166320.0 
        + ( 13860.0 - ( 462.0 - ( 132.0 - ( 99.0 - 140.0 
        / w2 ) / w2 ) / w2 ) / w2 ) / w / 166320.0;

      if ( alv <= t )
      {
        if ( 0.5 < pp )
        {
          ix = n - ix;
        }
        value = ix;
        return value;
      }
    }
  }
  return value;
}
//****************************************************************************80

int i4_max ( int i1, int i2 )

//****************************************************************************80
//
//  Purpose:
//
//    I4_MAX returns the maximum of two I4's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    13 October 1998
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int I1, I2, are two integers to be compared.
//
//    Output, int I4_MAX, the larger of I1 and I2.
//
{
  int value;

  if ( i2 < i1 )
  {
    value = i1;
  }
  else
  {
    value = i2;
  }
  return value;
}
//****************************************************************************80

int i4_min ( int i1, int i2 )

//****************************************************************************80
//
//  Purpose:
//
//    I4_MIN returns the minimum of two I4's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    13 October 1998
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int I1, I2, two integers to be compared.
//
//    Output, int I4_MIN, the smaller of I1 and I2.
//
{
  int value;

  if ( i1 < i2 )
  {
    value = i1;
  }
  else
  {
    value = i2;
  }
  return value;
}
//****************************************************************************80

void i4mat_print ( int m, int n, int a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    I4MAT_PRINT prints an I4MAT.
//
//  Discussion:
//
//    An I4MAT is an MxN array of I4's, stored by (I,J) -> [I+J*M].
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    10 September 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, the number of rows in A.
//
//    Input, int N, the number of columns in A.
//
//    Input, int A[M*N], the M by N matrix.
//
//    Input, string TITLE, a title.
//
{
  i4mat_print_some ( m, n, a, 1, 1, m, n, title );

  return;
}
//****************************************************************************80

void i4mat_print_some ( int m, int n, int a[], int ilo, int jlo, int ihi,
  int jhi, string title )

//****************************************************************************80
//
//  Purpose:
//
//    I4MAT_PRINT_SOME prints some of an I4MAT.
//
//  Discussion:
//
//    An I4MAT is an MxN array of I4's, stored by (I,J) -> [I+J*M].
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    20 August 2010
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, the number of rows of the matrix.
//    M must be positive.
//
//    Input, int N, the number of columns of the matrix.
//    N must be positive.
//
//    Input, int A[M*N], the matrix.
//
//    Input, int ILO, JLO, IHI, JHI, designate the first row and
//    column, and the last row and column to be printed.
//
//    Input, string TITLE, a title.
//
{
# define INCX 10

  int i;
  int i2hi;
  int i2lo;
  int j;
  int j2hi;
  int j2lo;

  cout << "\n";
  cout << title << "\n";

  if ( m <= 0 || n <= 0 )
  {
    cout << "\n";
    cout << "  (None)\n";
    return;
  }
//
//  Print the columns of the matrix, in strips of INCX.
//
  for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
  {
    j2hi = j2lo + INCX - 1;
    j2hi = i4_min ( j2hi, n );
    j2hi = i4_min ( j2hi, jhi );

    cout << "\n";
//
//  For each column J in the current range...
//
//  Write the header.
//
    cout << "  Col:";
    for ( j = j2lo; j <= j2hi; j++ )
    {
      cout << "  " << setw(6) << j - 1;
    }
    cout << "\n";
    cout << "  Row\n";
    cout << "\n";
//
//  Determine the range of the rows in this strip.
//
    i2lo = i4_max ( ilo, 1 );
    i2hi = i4_min ( ihi, m );

    for ( i = i2lo; i <= i2hi; i++ )
    {
//
//  Print out (up to INCX) entries in row I, that lie in the current strip.
//
      cout << setw(5) << i - 1 << ":";
      for ( j = j2lo; j <= j2hi; j++ )
      {
        cout << "  " << setw(6) << a[i-1+(j-1)*m];
      }
      cout << "\n";
    }
  }

  return;
# undef INCX
}
//****************************************************************************80

int *i4vec_multinomial_sample ( int n, double p[], int ncat )

//****************************************************************************80
//
//  Purpose:
//
//    I4VEC_MULTINOMIAL_SAMPLE generates a multinomial random deviate.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN77 version by Barry Brown, James Lovato.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Luc Devroye,
//    Non-Uniform Random Variate Generation,
//    Springer, 1986,
//    ISBN: 0387963057,
//    LC: QA274.D48.
//
//  Parameters:
//
//    Input, int N, the number of events, which will be
//    classified into one of the NCAT categories.
//
//    Input, double 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, int NCAT, the number of categories.
//
//    Output, int I4VEC_MULTINOMIAL_SAMPLE[NCAT], a random observation from 
//    the multinomial distribution.  All IX(i) will be nonnegative and their 
//    sum will be N.
//
{
  int i;
  int icat;
  int *ix;
  int ntot;
  double prob;
  double ptot;

  if ( n < 0 )
  {
    cerr << "\n";
    cerr << "I4VEC_MULTINOMIAL_SAMPLE - Fatal error!\n";
    cerr << "  N < 0\n";
    exit ( 1 );
  }

  if ( ncat <= 1 )
  {
    cerr << "\n";
    cerr << "I4VEC_MULTINOMIAL_SAMPLE - Fatal error!\n";
    cerr << "  NCAT <= 1\n";
    exit ( 1 );
  }

  for ( i = 0; i < ncat - 1; i++ )
  {
    if ( p[i] < 0.0 )
    {
      cerr << "\n";
      cerr << "I4VEC_MULTINOMIAL_SAMPLE - Fatal error!\n";
      cerr << "  Some P(i) < 0.\n";
      exit ( 1 );
    }

    if ( 1.0 < p[i] )
    {
      cerr << "\n";
      cerr << "I4VEC_MULTINOMIAL_SAMPLE - Fatal error!\n";
      cerr << "  Some 1 < P(i).\n";
      exit ( 1 );
    }
  }

  ptot = 0.0;
  for ( i = 0; i < ncat - 1; i++ )
  {
    ptot = ptot + p[i];
  }

  if ( 0.99999 < ptot ) 
  {
    cerr << "\n";
    cerr << "I4VEC_MULTINOMIAL_SAMPLE - Fatal error!\n";
    cerr << "  1.0 < Sum of P().\n";
    exit ( 1 );
  }
//
//  Initialize variables.
//
  ntot = n;
  ptot = 1.0;

  ix = new int[ncat];
  for ( i = 0; i < ncat; i++ )
  {
    ix[i] = 0;
  }
//
//  Generate the observation.
//
  for ( icat = 0; icat < ncat - 1; icat++ )
  {
    prob = p[icat] / ptot;
    ix[icat] = i4_binomial_sample ( ntot, prob );
    ntot = ntot - ix[icat];
    if ( ntot <= 0 )
    {
      return ix;
    }
    ptot = ptot - p[icat];
  }

  ix[ncat-1] = ntot;

  return ix;
}
//****************************************************************************80

void i4vec_transpose_print ( int n, int a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    I4VEC_TRANSPOSE_PRINT prints an I4VEC "transposed".
//
//  Discussion:
//
//    An I4VEC is a vector of I4's.
//
//  Example:
//
//    A = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 }
//    TITLE = "My vector:  "
//
//    My vector:      1    2    3    4    5
//                    6    7    8    9   10
//                   11
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    03 July 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of components of the vector.
//
//    Input, int A[N], the vector to be printed.
//
//    Input, string TITLE, a title.
//
{
  int i;
  int ihi;
  int ilo;
  int title_len;

  title_len = title.length ( );

  for ( ilo = 1; ilo <= n; ilo = ilo + 5 )
  {
    ihi = i4_min ( ilo + 5 - 1, n );
    if ( ilo == 1 )
    {
      cout << title;
    }
    else
    {
      for ( i = 1; i <= title_len; i++ )
      {
        cout << " ";
      }
    }
    for ( i = ilo; i <= ihi; i++ )
    {
      cout << setw(12) << a[i-1];
    }
    cout << "\n";
  }

  return;
}
//****************************************************************************80

int *i4vec_zero_new ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    I4VEC_ZERO_NEW creates and zeroes an I4VEC.
//
//  Discussion:
//
//    An I4VEC is a vector of I4's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 July 2008
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of entries in the vector.
//
//    Output, int I4VEC_ZERO_NEW[N], a vector of zeroes.
//
{
  int *a;
  int i;

  a = new int[n];

  for ( i = 0; i < n; i++ )
  {
    a[i] = 0;
  }
  return a;
}
//****************************************************************************80

void input_print ( int chain_num, int cr_num, string input_filename, 
  double gr_threshold, int jumpstep, double limits[], int gen_num, 
  string ob_filename, int ob_num, int pair_num, int par_num, int printstep, 
  string restart_string )

//****************************************************************************80
//
//  Purpose:
//
//    INPUT_PRINT prints the data from the input file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int CHAIN_NUM, the total number of chains.
//
//    Input, int CR_NUM, the total number of CR values.
//
//    Input, string INPUT_FILENAME, the name of the input file.
//
//    Input, double GR_THRESHOLD, the convergence tolerance for the
//    Gelman-Rubin statistic.
//
//    Input, int JUMPSTEP, forces a "long jump" every
//    JUMPSTEP generations.
//
//    Input, double LIMITS[2*PAR_NUM], lower and upper limits
//    for each parameter.
//
//    Input, int GEN_NUM, the total number of generations.
//
//    Input, string OB_FILENAME, the name of the observation file.
//
//    Input, int OB_NUM, the total number of observations.
//
//    Input, int PAIR_NUM, the number of pairs of 
//    crossover chains.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input, int PRINTSTEP, the interval between generations on 
//    which the Gelman-Rubin statistic will be computed and written to a file.
//
//    Input, string RESTART_STRING, indicates whether a restart
//    data file should be used to initialize the chains.
//
{
  int j;

  cout << "\n";
  cout << "INPUT_PRINT:\n";
  cout << "  Input data read from \"" << input_filename << "\".\n";
  cout << "\n";
  cout << "  Number of parameters\n";
  cout << "  PAR_NUM = " << par_num << "\n";
  cout << "\n";
  cout << "  LIMITS: Lower and upper limits for each parameter:\n";
  cout << "\n";
  cout << "  Index           Lower           Upper\n";
  cout << "\n";
  for ( j = 0; j < par_num; j++ )
  {
    cout << "  " << setw(5) << j
         << "  " << setw(14) << limits[0+j*2]
         << "  " << setw(14) << limits[1+j*2] << "\n";
  }
  cout << "\n";
  cout << "  Number of generations:\n";
  cout << "  GEN_NUM = " << gen_num << "\n";
  cout << "\n";
  cout << "  Number of simultaneous chains:\n";
  cout << "  CHAIN_NUM = " << chain_num << "\n";
  cout << "\n";
  cout << "  Observation filename:\n";
  cout << "  OB_FILENAME = \"" << ob_filename << "\".\n";
  cout << "\n";
  cout << "  Number of observations:\n";
  cout << "  OB_NUM = " << ob_num << "\n";
  cout << "\n";
  cout << "  Number of pairs of chains for crossover:\n";
  cout << "  PAIR_NUM = " << pair_num << "\n";
  cout << "\n";
  cout << "  Number of crossover values:\n";
  cout << "  CR_NUM = " << cr_num << "\n";
  cout << "\n";
  cout << "  Number of steps til a long jump:\n";
  cout << "  JUMPSTEP = " << jumpstep << "\n";
  cout << "\n";
  cout << "  Interval between Gelman-Rubin computations:\n";
  cout << "  PRINTSTEP = " << printstep << "\n";
  cout << "\n";
  cout << "  Gelman-Rubin convergence tolerance:\n";
  cout << "  GR_THRESHOLD = " << gr_threshold << "\n";
  cout << "\n";
  cout << "  Restart string (Y/N):\n";
  cout << "  RESTART_STRING = \"" << restart_string << "\".\n";

  return;
}
//****************************************************************************80

void input_read ( int &chain_num, int &cr_num, string input_filename, 
  double &gr_threshold, int &jumpstep, double limits[], int &gen_num, 
  string *ob_filename, int &ob_num, int &pair_num, int par_num, int &printstep,
  string *restart_string )

//****************************************************************************80
//
//  Purpose:
//
//    INPUT_READ reads the data from the input file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    03 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Output, int &CHAIN_NUM, the total number of chains.
//
//    Output, int &CR_NUM, the total number of CR values.
//
//    Input, string INPUT_FILENAME, the name of the input file.
//
//    Output, double &GR_THRESHOLD, the convergence tolerance for
//    the Gelman-Rubin statistic.
//
//    Output, int &JUMPSTEP, forces a "long jump" every
//    JUMPSTEP generations.
//
//    Output, double LIMITS[2*PAR_NUM], lower and upper limits
//    for each parameter.
//
//    Output, int &GEN_NUM, the total number of generations.
//
//    Output, string *OB_FILENAME, the name of the 
//    observation file.
//
//    Output, int &OB_NUM, the total number of observations.
//
//    Output, int &PAIR_NUM, the number of pairs of 
//    crossover chains.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Output, int &PRINTSTEP, the interval between generations on 
//    which the Gelman-Rubin statistic will be computed and written to a file.
//
//    Output, string *RESTART_STRING, indicates whether a restart
//    data file should be used to initialize the chains.
//
{
  ifstream input;
  int j;
  double l;
  string line;
  double u;

  input.open ( input_filename.c_str ( ) );

  if ( !input )
  {
    cerr << "\n";
    cerr << "INPUT_READ - Fatal error!\n";
    cerr << "  Could not open the file \"" << input_filename << "\".\n";
    exit ( 1 );
  }
//
//  Read and ignore line 1.
//
  getline ( input, line );
//
//  Read and ignore line 2.
//
  getline ( input, line );
//
//  Read and ignore line 3.
//
  getline ( input, line );
//
//  Line 4 is actually PAR_NUM lines,
//  each a pair of values.
//
  for ( j = 0; j < par_num; j++ )
  {
    input >> l >> u;
    limits[0+j*2] = l;
    limits[1+j*2] = u;
  }
  getline ( input, line );
//
//  Read and ignore line 5.
//
  getline ( input, line );
//
//  Line 6 is GEN_NUM.
//
  input >> gen_num;
  getline ( input, line );
//
//  Read and ignore line 7.
//
  getline ( input, line );
//
//  Line 8 is CHAIN_NUM.
//
  input >> chain_num;
  getline ( input, line );
//
//  Read and ignore line 9.
//
  getline ( input, line );
//
//  Line 10 is OB_FILENAME.
//
  getline ( input, *ob_filename );
//
//  Read and ignore line 11.
//
  getline ( input, line );
//
//  Line 12 is OB_NUM.
//
  input >> ob_num;
  getline ( input, line );
//
//  Read and ignore line 13.
//
  getline ( input, line );
//
//  Line 14 is PAIR_NUM.
//
  input >> pair_num;
  getline ( input, line );
//
//  Read and ignore line 15.
//
  getline ( input, line );
//
//  Line 16 is CR_NUM.
//
  input >> cr_num;
  getline ( input, line );
//
//  Read and ignore line 17.
//
  getline ( input, line );
//
//  Line 18 is JUMPSTEP.
//
  input >> jumpstep;
  getline ( input, line );
//
//  Read and ignore line 19.
//
  getline ( input, line );
//
//  Line 20 is PRINTSTEP.
//
  input >> printstep;
  getline ( input, line );
//
//  Read and ignore line 21.
//
  getline ( input, line );
//
//  Line 22 is GR_THRESHOLD.
//
  input >> gr_threshold;
  getline ( input, line );
//
//  Read and ignore line 23.
//
  getline ( input, line );
//
//  Line 24 is RESTART_STRING.
//
  getline ( input, *restart_string );

  input.close ( );

  return;
}
//****************************************************************************80

void input_size ( string input_filename, int &par_num )

//****************************************************************************80
//
//  Purpose:
//
//    INPUT_SIZE reads the parameter size variable from the input file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    03 May 2013
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string INPUT_FILENAME, the name of the input file.
//
//    Output, int &PAR_NUM, the total number of parameters.
//
{
  ifstream input;
  string line;
 
  input.open ( input_filename.c_str ( ) );

  if ( ! input )
  {
    cerr << "\n";
    cerr << "INPUT_SIZE - Fatal error!\n";
    cerr << "  Could not open the file \"" << input_filename << "\".\n";
    exit ( 1 );
  }
//
//  Read and ignore first line.
//
  getline ( input, line );
//
//  PAR_NUM is on the second line.
//
  input >> par_num;
 
  input.close ( );

  return;
}
//****************************************************************************80

void jumprate_choose ( double cr[], int cr_index, int cr_num, int gen_index,
  int jump_dim[], int &jump_num, double &jumprate, double jumprate_table[],
  int jumpstep, int par_num )

//****************************************************************************80
//
//  Purpose:
//
//    JUMPRATE_CHOOSE chooses a jump rate from the jump rate table.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Jasper Vrugt, CJF ter Braak, CGH Diks, Bruce Robinson, James Hyman, 
//    Dave Higdon,
//    Accelerating Markov Chain Monte Carlo Simulation by Differential 
//    Evolution with Self-Adaptive Randomized Subspace Sampling,
//    International Journal of Nonlinear Sciences and Numerical Simulation,
//    Volume 10, Number 3, March 2009, pages 271-288.
//
//  Parameters:
//
//    Input, double CR[CR_NUM], the CR values.
//
//    Input, int CR_INDEX, the index of the CR.
//    1 <= CR_INDEX <= CR_NUM.
//
//    Input, int CR_NUM, the total number of CR values.
//
//    Input, int GEN_INDEX, the current generation.
//    1 <= GEN_INDEX <= GEN_NUM.
//
//    Output, int JUMP_DIM[PAR_NUM], the indexes of the
//    parameters to be updated.
//
//    Output, int &JUMP_NUM, the number of dimensions in which
//    a jump will be made.  0 <= JUMP_NUM <= PAR_NUM.
//
//    Output, double &JUMPRATE, the jump rate.
//
//    Input, double JUMPRATE_TABLE[PAR_NUM], the jump rate table.
//
//    Input, int JUMPSTEP, forces a "long jump" every
//    JUMPSTEP generations.
//
//    Input, int PAR_NUM, the total number of parameters.
//
{
  int i;
  double r;
//
//  Determine the dimensions that will be updated.
//
  jump_num = 0;
  for ( i = 0; i < par_num; i++ )
  {
    jump_dim[i] = 0;
  }

  for ( i = 0; i < par_num; i++ )
  {
    r = r8_uniform_01_sample ( );

    if ( 1.0 - cr[cr_index] < r )
    {
      jump_dim[jump_num] = i;
      jump_num = jump_num + 1;
    }
  }
//
//  Calculate the general jump rate.
//
  if ( jump_num == 0 )
  {
    jumprate = 0.0;
  }
  else
  {
    jumprate = jumprate_table[jump_num-1];
  }
//
//  If parameter dimension is 1, 2, or 3, fix the jump rate to 0.6.
//
  if ( par_num <= 3 )
  {
    jumprate = 0.6;
  }
//
//  Determine if a long jump is forced.
//
  if ( ( gen_index % jumpstep ) == 0 )
  {
    jumprate = 0.98;
  }
  return;
}
//****************************************************************************80

double *jumprate_table_init ( int pair_num, int par_num )

//****************************************************************************80
//
//  Purpose:
//
//    JUMPRATE_TABLE_INIT initializes the jump rate table.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, int PAIR_NUM, the number of pairs of 
//    crossover chains.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Output, double JUMPRATE_TABLE_INIT[PAR_NUM], the jumprate table.
//
{
  double c;
  int i;
  double *jumprate_table;
 
  jumprate_table = new double[par_num];

  c = 2.38 / sqrt ( ( double ) ( 2 * pair_num ) );

  for ( i = 0; i < par_num; i++ )
  {
    jumprate_table[i] = c / sqrt ( ( double ) ( i + 1 ) );
  }

  return jumprate_table;
}
//****************************************************************************80

void jumprate_table_print ( double jumprate_table[], int pair_num, int par_num )

//****************************************************************************80
//
//  Purpose:
//
//    JUMPRATE_TABLE_PRINT prints the jump rate table.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double JUMPRATE_TABLE[PAR_NUM], the jumprate table.
//
//    Input, int PAIR_NUM, the number of pairs of 
//    crossover chains.
//
//    Input, int PAR_NUM, the total number of parameters.
//
{
  int i;
 
  cout << "\n";
  cout << "JUMPRATE_TABLE_PRINT\n";
  cout << "\n";
  cout << "   I        Jumprate\n";
  cout << "\n";
  for ( i = 0; i < par_num; i++ )
  {
    cout << "  " << setw(2) << i
         << "  " << setw(14) << jumprate_table[i] << "\n";
  }
  return;
}
//****************************************************************************80

void mnor_init ( double cov_mat[], double &logdet, int mnor_ind[], 
  double mnor_mean[], int mnor_num, int par_num, double parm[], 
  double prior_den_par_val[] )

//****************************************************************************80
//
//  Purpose:
//
//    MNOR_INIT initializes information for multivariate normal prior densities.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Output, double COV_MAT[MNOR_NUM*MNOR_NUM], the covariance
//    matrix for a multivariate normal density.
//
//    Output, double &LOGDET, the logarithm of the determinant
//    of the covariance matrix.
//
//    Input, int MNOR_IND[MNOR_NUM], the indices of
//    parameters with multivariate density.
//
//    Output, int MNOR_MEAN[MNOR_NUM], the means associated
//    with the multivariate density.
//
//    Input, int MNOR_NUM, the number of parameters with
//    a multivariate normal distribution.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Output, double PARM[MNOR_NUM*(MNOR_NUM+3)/2+1], parameters 
//    for a multivariate normal distribution.
//
//    Input, double PRIOR_DEN_PAR_VAL[2*PAR_NUM], parameter
//    values associated with prior densities.
//
{
  double a;
  int i;
  int j;
//
//  Collect the mean values of the multinormal density.
//
  for ( i = 0; i < mnor_num; i++ )
  {
    j = mnor_ind[i];
    a = prior_den_par_val[0+j*2];
    mnor_mean[i] = a;
  }
//
//  Pack data into PARM.
//
  mnor_process ( cov_mat, logdet, mnor_mean, mnor_num, parm );

  return;
}
//****************************************************************************80

void mnor_process ( double cov_mat[], double &logdet, double mnor_mean[], 
  int mnor_num, double parm[] )

//****************************************************************************80
//
//  Purpose:
//
//    MNOR_PROCESS sets data for the generation of multivariate normal deviates.
//
//  Discussion:
//
//    The final form and purpose of this routine is not decided yet.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN77 version by Barry Brown, James Lovato.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input/output, double COVM[MNOR_NUM*MNOR_NUM].  On input, 
//    the covariance matrix of the multivariate distribution.  On output, 
//    the information has been overwritten.
//
//    Output, double &LOGDET, the logarithm of the determinant of the
//    covariance matrix.
//
//    Input, double MNOR_MEAN[MNOR_NUM], the means of the multivariate 
//    normal distribution.
//
//    Input, int MNOR_NUM, the number of parameters with
//    a multivariate normal distribution.
//
//    Output, double PARM[MNOR_NUM*(MNOR_NUM+3)/2+1], parameters
//    needed to generate multivariate normal deviates.
//
{
  int i;
  int icount;
  int info;
  int j;
//
//  Store MNOR_NUM.
//
  parm[0] = mnor_num;
//
//  Store MEANV.
//
  for ( i = 1; i < mnor_num; i++ )
  {
    parm[i] = mnor_mean[i-1];
  }
//
//  Compute the Cholesky decomposition.
//
  info = r8mat_pofa ( cov_mat, mnor_num, mnor_num );

  if ( info != 0 )
  {
    cerr << "\n";
    cerr << "MNOR_PROCESS - Fatal error!\n";
    cerr << "  Matrix not positive definite.\n";
    exit ( 1 );
  }
//
//  Store the upper half of the Cholesky factor.
//
  icount = mnor_num + 1;

  for ( i = 0; i < mnor_num; i++ )
  {
    for ( j = i; j < mnor_num; j++ )
    {
      parm[icount] = cov_mat[i+j*mnor_num];
      icount = icount + 1;
    }
  }
//
//  Compute the determinant...NOT DONE!
//
  logdet = 1.0;

  return;
}
//****************************************************************************80

void mnor_read ( double cov_mat[], string mnor_filename, int mnor_num )

//****************************************************************************80
//
//  Purpose:
//
//    MNOR_READ reads a file for multivariate normal prior densities.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Output, double COV_MAT[MNOR_NUM*MNOR_NUM], the covariance
//    matrix for a multivariate normal density.
//
//    Input, string MNOR_FILENAME, the name of the file 
//    containing the covariance matrix information.
//
//    Input, int MNOR_NUM, the number of parameters with
//    a multivariate normal distribution.
//
{
  int col;
  int cov_dim;
  ifstream cov;
  int index;
  string line;
  int row;

  cov.open ( mnor_filename.c_str ( ) );

  if ( ! cov )
  {
    cerr << "\n";
    cerr << "MNOR_READ - Fatal error!\n";
    cerr << "  Could not open the file \"" << mnor_filename << "\".\n";
    exit ( 1 );
  }
//
//  Read and ignore the first line.
//
  getline ( cov, line );
//
//  Get the dimension from the second line, then discard the line.
//
  cov >> cov_dim;

  if ( cov_dim != mnor_num )
  {
    cerr << "\n";
    cerr << "MNOR_READ - Fatal error!\n";
    cerr << "  Dimension mismatch.\n";
    cerr << "  Expecting matrix dimension MNOR_NUM = " << mnor_num << "\n";
    cerr << "  Found matrix dimension COV_DIM = " << cov_dim << "\n";
    exit ( 1 );
  }
//
//  Read the matrix.
//
  for ( row = 0; row < mnor_num; row++ )
  {
    for ( col = 0; col < mnor_num; col++ )
    {
      index = row + col * mnor_num;
      cov >> cov_mat[index];
    }
  }

  cov.close ( );

  return;
}
//****************************************************************************80

void ob_read ( double ob_data[], string ob_filename, int ob_num )

//****************************************************************************80
//
//  Purpose:
//
//    OB_READ reads observational data from a file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Output, double OB_DATA[OB_NUM], the data.
//
//    Input, string OB_FILENAME, the name of the file.
//
//    Input, int OB_NUM, the total number of observations.
//
{
  int i;
  string line;
  ifstream ob;

  ob.open ( ob_filename.c_str ( ) );

  if ( ! ob )
  {
    cerr << "\n";
    cerr << "OB_READ - Fatal error!\n";
    cerr << "  Could not open file \"" << ob_filename << "\".\n";
    exit ( 1 );
  }
//
//  Read and ignore the header line.
//
  getline ( ob, line );
//
//  Read OB_NUM lines of data.
//
  for ( i = 0; i < ob_num; i++ )
  {
    ob >> ob_data[i];
  }

  ob.close ( );

  return;
}
//****************************************************************************80

double prior_density ( double cov_mat[], double logdet, int mnor_ind[], 
  double mnor_mean[], int mnor_num, int par_num, int prior_den_par_num[], 
  double prior_den_par_val[], int prior_den_type[], int uni_ind[], int uni_num, 
  double zp[] )

//****************************************************************************80
//
//  Purpose:
//
//    PRIOR_DENSITY evaluates the prior density function.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double COV_MAT[MNOR_NUM*MNOR_NUM], the covariance 
//    matrix for a multivariate normal density.
//
//    Input, double LOGDET, the logarithm of the determinant
//    of the multivariate normal covariance matrix.
//
//    Input, int MNOR_IND[MNOR_NUM], the indices of
//    parameters with a multivariate normal distribution.
//
//    Input, double MNOR_MEAN[PAR_NUM], the mean of
//    parameters with a multivariate normal distribution.
//
//    Input, int MNOR_NUM, the number of parameters with
//    a multivariate normal distribution.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input, int PRIOR_DEN_PAR_NUM[PAR_NUM], the number of
//    variables associated with the prior density for each parameter.
//
//    Input, double PRIOR_DEN_PAR_VAL[2*PAR_NUM], parameter values
//    associated with prior densities.
//
//    Input, int PRIOR_DEN_TYPE[PAR_NUM], the type
//    of the prior density.
//
//    Input, int UNI_IND[UNI_NUM], indices of parameters 
//    associated with univariate distributions.
//
//    Input, int UNI_NUM, the number of parameters with
//    a univariate distribution.
//
//    Input, double ZP[PAR_NUM], the argument of the density
//    function.
//
//    Output, real PRIOR_DENSITY, the value of the prior density function.
//
{
  double a;
  double b;
  int i;
  int j;
  double *tmp;
  int tp;
  double value;

  value = 1.0;
//
//  Univariate density:
//
  for ( i = 0; i < uni_num; i++ )    
  {
    j = uni_ind[i];
    tp = prior_den_type[j];

    if ( 1 <= prior_den_par_num[j] )
    {
      a = prior_den_par_val[0+j*2];
    }
    else
    {
      a = 0.0;
    }

    if ( 2 <= prior_den_par_num[j] )
    {
      b = prior_den_par_val[1+j*2];
    }
    else
    {
      b = 0.0;
    }

    if ( tp == 1 )
    {
      value = value * r8_beta_pdf ( a, b, zp[j] );
    }
    else if ( tp == 2 )
    {
      value = value * r8_chi_pdf ( a, zp[j] );
    }
    else if ( tp == 3 )
    {
      value = value * r8_exponential_pdf ( a, zp[j] );
    }
    else if ( tp == 4 )
    {
      value = value * r8_gamma_pdf ( a, b, zp[j] );
    }
    else if ( tp == 5 )
    {
      value = value * r8_invchi_pdf ( a, zp[j] );
    }
    else if ( tp == 6 )
    {
      value = value * r8_invgam_pdf ( a, b, zp[j] );
    }
    else if ( tp == 7 )
    {
      value = value * r8_normal_pdf ( a, b, zp[j] );
    }
    else if ( tp == 8 )
    {
      value = value * r8_scinvchi_pdf ( a, b, zp[j] );
    }
    else if ( tp == 9 )
    {
      value = value * r8_uniform_pdf ( a, b, zp[j] );
    }
    else
    {
      cerr << "\n";
      cerr << "PRIOR_DENSITY - Fatal error!\n";
      cerr << "  Unrecognized type for prior density.\n";
      exit ( 1 );
    }
  }
//
//  Calculate density values for multinormal density     
//
  if ( 1 < mnor_num )
  {
    tmp = new double[mnor_num];

    for ( i = 0; i < mnor_num; i++ )
    {
      tmp[i] = zp[mnor_ind[i]];
    }
//
//  I would assume R8_MNOR_PDF would be accessing PARM, not these pieces.
//
    value = r8_mnor_pdf ( mnor_num, mnor_mean, cov_mat, logdet, tmp ) ;

    delete [] tmp;
  }

  return value;
}
//****************************************************************************80

void prior_print ( int mnor_ind[], int mnor_num, int par_num, 
  int prior_den_par_num[], double prior_den_par_val[], int prior_den_type[], 
  string prior_filename, int uni_ind[], int uni_num )

//****************************************************************************80
//
//  Purpose:
//
//    PRIOR_PRINT prints information about the prior density.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    03 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, int MNOR_IND[MNOR_NUM], the indices of parameters
//    with a multivariate normal distribution.
//
//    Input, int MNOR_NUM, the number of parameters with
//    a multivariate normal distribution.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input, int PRIOR_DEN_PAR_NUM[PAR_NUM], the number of
//    variables associated with the prior density for each parameter.
//
//    Input, double PRIOR_DEN_PAR_VAL[2*PAR_NUM], parameter values
//    associated with prior densities.
//
//    Input, int PRIOR_DEN_TYPE[PAR_NUM], the type
//    of the prior density.
//
//    Input, string PRIOR_FILENAME, the name of the file
//    containing information about the prior density.
//
//    Input, int UNI_IND[UNI_NUM], the indices of variables
//    with univariate density.
//
//    Input, int UNI_NUM, the number of parameters with
//    a univariate distribution.
//
{
  int i;
  int j;
  int k;
  int pn;
  int pt;
  string type_string;

  cout << "\n";
  cout << "PRIOR_PRINT:\n";
  cout << "  Prior density information read from \"" 
       << prior_filename << "\".\n";
  cout << "\n";
  cout << "  Number of univariate parameters = " << uni_num << "\n";
  cout << "\n";
  cout << "  Full   Uni\n";
  cout << "  Index  Index               Type     #";
  cout << "            Arg1            Arg2\n";
  cout << "\n";

  for ( i = 0; i < uni_num; i++ )
  {
    j = uni_ind[i];
    pt = prior_den_type[j];

    if ( pt == 1 )
    {
      type_string = "beta";
    }
    else if ( pt == 2 )
    {
      type_string = "chi-square";
    }
    else if ( pt == 3 )
    {
      type_string = "exponential";
    }
    else if ( pt == 4 )
    {
      type_string = "gamma";
    }
    else if ( pt == 5 )
    {
      type_string = "inv-chi-square";
    }
    else if ( pt == 6 )
    {
      type_string = "inv-gamma";
    }
    else if ( pt == 7 )
    {
      type_string = "normal";
    }
    else if ( pt == 8 )
    {
      type_string = "scaled-inv-chi-square";
    }
    else if ( pt == 9 )
    {
      type_string = "uniform";
    }
    else if ( pt == 10 )
    {
      type_string = "multinormal";
    }
    else
    {
      type_string = "?";
    }
    pn = prior_den_par_num[j];
    cout << "  " << setw(5) << i
         << "  " << setw(5) << j
         << "  " << setw(20) << type_string
         << "  " << setw(1) << pn;
    for ( k = 0; k < pn; k++ )
    {
      cout << "  " << setw(14) <<  prior_den_par_val[k+j*2];
    }
    cout << "\n";
  }

  cout << "\n";
  cout << "  Number of multivariate parameters = " << mnor_num << "\n";
  if ( 0 < mnor_num )
  {
    cout << "\n";
    cout << "  Index  Parameter index\n";
    cout << "\n";
    for ( i = 0; i < mnor_num; i++ )
    {
      j = mnor_ind[i];
      cout << "  " << setw(5) << i
           << "  " << setw(2) << j << "\n";
    }
  }
  return;
}
//****************************************************************************80

void prior_read ( int mnor_ind[], int &mnor_num, int par_num,
  int prior_den_par_num[], double prior_den_par_val[], int prior_den_type[],
  string prior_filename, int uni_ind[], int &uni_num )

//****************************************************************************80
//
//  Purpose:
//
//    PRIOR_READ reads information about the prior density.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Output, int MNOR_IND[MNOR_NUM], the indices of parameters
//    with a multivariate normal distribution.
//
//    Output, int &MNOR_NUM, the number of parameters with
//    a multivariate normal distribution.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Output, int PRIOR_DEN_PAR_NUM[PAR_NUM], the number of
//    variables associated with the prior density for each parameter.
//
//    Output, double PRIOR_DEN_PAR_VAL[2*PAR_NUM], parameter values
//    associated with prior densities.
//
//    Output, int PRIOR_DEN_TYPE[PAR_NUM], the type
//    of the prior density.
//
//    Input, string PRIOR_FILENAME, the name of the file
//    containing information about the prior density.
//
//    Output, int UNI_IND[UNI_NUM], the indices of variables
//    with univariate density.
//
//    Output, int &UNI_NUM, the number of parameters with
//    a univariate distribution.
//
{
  int k;
  string line;
  int par_index;
  int pn;
  double pv;
  ifstream prior;
  string type_string;

  mnor_num = 0;
  uni_num = 0;

  prior.open ( prior_filename.c_str ( ) );

  if ( ! prior )
  {
    cerr << "\n";
    cerr << "PRIOR_READ - Fatal error!\n";
    cerr << "  Could not open the file \"" << prior_filename << "\".\n";
    exit ( 1 );
  }
//
//  Read and ignore line 1.
//
  getline ( prior, line );
//
//  Read prior information for each parameter dimension.
//
  for ( par_index = 0; par_index < par_num; par_index++ )
  {
    prior >> type_string;

    if ( type_string.compare ( "beta" ) == 0 )
    {
      k = 1;
    }
    else if ( type_string.compare ( "chi-square" ) == 0 )
    {
      k = 2;
    }
    else if ( type_string.compare ( "exponential" ) == 0 )
    {
      k = 3;
    }
    else if ( type_string.compare ( "gamma" ) == 0 )
    {
      k = 4;
    }
    else if ( type_string.compare ( "inv-chi-square" ) == 0 )
    {
      k = 5;
    }
    else if ( type_string.compare ( "inv-gamma" ) == 0 )
    {
      k = 6;
    }
    else if ( type_string.compare ( "normal" ) == 0 )
    {
      k = 7;
    }
    else if ( type_string.compare ( "scaled-inv-chi-square" ) == 0 )
    {
      k = 8;
    }
    else if ( type_string.compare ( "uniform" ) == 0 )
    {
      k = 9;
    }
    else if ( type_string.compare ( "multinormal" ) == 0 )
    {
      k = 10;
    }
    else
    {
      cerr << "\n";
      cerr << "PRIOR_READ - Fatal error!\n";
      cerr << "  Unrecognized input value of type string.\n";
      exit ( 1 );
    }

    prior_den_type[par_index] = k;
//
//  Read parameters for the prior density function.
//
    prior >> pn;
    prior_den_par_num[par_index] = pn;

    if ( 1 <= pn )
    {
      prior >> pv;
      prior_den_par_val[0+par_index*2] = pv;
    }
    if ( 2 <= pn )
    {
      prior >> pv;
      prior_den_par_val[1+par_index*2] = pv;
    }

    if ( k == 10 )
    {
      mnor_ind[mnor_num] = par_index;
      mnor_num = mnor_num + 1;
    }
    else
    {
      uni_ind[uni_num] = par_index;
      uni_num = uni_num + 1;
    }
  }

  prior.close ( );

  return;
}
//****************************************************************************80

double *prior_sample ( int mnor_ind[], int mnor_num, int par_num, 
  double parm[], int prior_den_par_num[], double prior_den_par_val[], 
  int prior_den_type[], int uni_ind[], int uni_num )

//****************************************************************************80
//
//  Purpose:
//
//    PRIOR_SAMPLE samples from the prior distribution.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, int MNOR_IND[MNOR_NUM], indices of parameters
//    associated with multivariate normal distributions.
//
//    Input, int MNOR_NUM, the number of parameters with
//    a multivariate normal distribution.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input, double PARM[MNOR_NUM*(MNOR_NUM+3)/2+1], parameters 
//    for a multivariate normal distribution.
//
//    Input, int PRIOR_DEN_PAR_NUM[PAR_NUM], the number of
//    variables associated with the prior density for each parameter.
//
//    Input, double PRIOR_DEN_PAR_VAL[2*PAR_NUM], parameter values
//    associated with prior densities.
//
//    Input, int PRIOR_DEN_TYPE[PAR_NUM], the type
//    of the prior density.
//
//    Input, int UNI_IND[UNI_NUM], indices of parameters
//    associated with univariate distributions.
//
//    Input, int UNI_NUM, the number of parameters with
//    a univariate distribution.
//
//    Output, double PRIOR_SAMPLE[PAR_NUM], the sample from the distribution.
//
{
  double a;
  double a1;
  double b;
  double b1;
  int i;
  int j;
  int pt;
  double *tmp;
  double *zp;
//
//  Generate samples following the univariate density.
//
  zp = new double[par_num];

  for ( i = 0; i < uni_num; i++ )
  {
    j = uni_ind[i];
    pt = prior_den_type[j];

    if ( 1 <= prior_den_par_num[j] )
    {
      a = prior_den_par_val[0+j*2];
    }
    else
    {
      a = 0.0;
    }

    if ( 2 <= prior_den_par_num[j] )
    {
      b = prior_den_par_val[1+j*2];
    }
    else
    {
      b = 0.0;
    }

    if ( pt == 1 )
    {
      zp[j] = r8_beta_sample ( a, b );
    }
    else if ( pt == 2 )
    {
      zp[j] = r8_chi_sample ( a );
    }
    else if ( pt == 3 )
    {
      zp[j] = r8_exponential_sample ( a );
    }
    else if ( pt == 4 )
    {
      zp[j] = r8_gamma_sample ( a, b );
    }
    else if ( pt == 5 )
    {
      a1 = 0.5;
      b1 = 0.5 * a;
      zp[j] = r8_gamma_sample ( a1, b1 );
      if ( zp[j] != 0.0 )
      {
        zp[j] = 1.0 / zp[j];
      }
    }
    else if ( pt == 6 )
    {
      zp[j] = r8_gamma_sample ( a, b );
      if ( zp[j] != 0.0 )
      {
        zp[j] = 1.0 / zp[j];
      }
    }
    else if ( pt == 7 )
    {
      zp[j] = r8_normal_sample ( a, b );
    }
    else if ( pt == 8 )
    {
      a1 = 0.5 * a * b;
      b1 = 0.5 * a;
      zp[j] = r8_gamma_sample ( a1, b1 );
      if ( zp[j] != 0.0 )
      {
        zp[j] = 1.0 / zp[j];
      }
    }
    else if ( pt == 9 )
    {
      zp[j] = r8_uniform_sample ( a, b );
    }
    else
    {
      cerr << "\n";
      cerr << "PRIOR_SAMPLE - Fatal error!\n";
      cerr << "  Unrecognized type for prior density.\n";
      exit ( 1 );
    }
  }
//
//  Generate samples following a multinormal density.
//
  if ( 1 < mnor_num )
  {      
    tmp = r8vec_mnor_sample ( parm );
  
    for ( i = 0; i < mnor_num; i++ )
    {
      j = mnor_ind[i];
      zp[j] = tmp[i];
    }  
    delete [] tmp;    
  }

  return zp;
}
//****************************************************************************80

double r8_beta_pdf ( double alpha, double beta, double rval )

//****************************************************************************80
//
//  Purpose:
//
//    R8_BETA_PDF evaluates the PDF of a beta distribution.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double ALPHA, BETA, shape parameters.
//    0.0 < ALPHA, BETA.
//
//    Input, double RVAL, the point where the PDF is evaluated.
//
//    Output, double R8_BETA_PDF, the value of the PDF at RVAL.
//
{
  double temp;
  double value;

  if ( alpha <= 0.0 )
  {
    cerr << "\n";
    cerr << "R8_BETA_PDF - Fatal error!\n";
    cerr << "  Parameter ALPHA is not positive.\n";
    exit ( 1 );
  }

  if ( beta <= 0.0 )
  {
    cerr << "\n";
    cerr << "R8_BETA_PDF- Fatal error!\n";
    cerr << "  Parameter BETA is not positive.\n";
    exit ( 1 );
  }

  if ( rval < 0.0 || 1.0 < rval )
  {
    value = 0.0;
  }
  else
  {
    temp = r8_gamma_log ( alpha + beta ) - r8_gamma_log ( alpha ) 
      - r8_gamma_log ( beta );

    value = exp ( temp ) * pow ( rval, alpha - 1.0 )
      * pow ( 1.0 - rval, beta - 1.0 );
  }
  return value;
}
//****************************************************************************80

double r8_beta_sample ( double aa, double bb )

//****************************************************************************80
//
//  Purpose:
//
//    R8_BETA_SAMPLE generates a beta random deviate.
//
//  Discussion:
//
//    This procedure returns a single random deviate from the beta distribution
//    with parameters A and B.  The density is
//
//      x^(a-1) * (1-x)^(b-1) / Beta(a,b) for 0 < x < 1
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN77 version by Barry Brown, James Lovato.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Russell Cheng,
//    Generating Beta Variates with Nonintegral Shape Parameters,
//    Communications of the ACM,
//    Volume 21, Number 4, April 1978, pages 317-322.
//
//  Parameters:
//
//    Input, double AA, the first parameter of the beta distribution.
//    0.0 < AA.
//
//    Input, double BB, the second parameter of the beta distribution.
//    0.0 < BB.
//
//    Output, double R8_BETA_SAMPLE, a beta random variate.
//
{
  double a;
  double alpha;
  double b;
  double beta;
  double delta;
  double gamma;
  double k1;
  double k2;
  static double log4 = 1.3862943611198906188;
  static double log5 = 1.6094379124341003746;
  double r;
  double s;
  double t;
  double u1;
  double u2;
  double v;
  double value;
  double w;
  double y;
  double z;

  if ( aa <= 0.0 )
  {
    cerr << "\n";
    cerr << "R8_BETA_SAMPLE - Fatal error!\n";
    cerr << "  AA <= 0.0\n";
    exit ( 1 );
  }

  if ( bb <= 0.0 ) 
  {
    cerr << "\n";
    cerr << "R8_BETA_SAMPLE - Fatal error!\n";
    cerr << "  BB <= 0.0\n";
    exit ( 1 );
  }
//
//  Algorithm BB
//
  if ( 1.0 < aa && 1.0 < bb )
  {
    a = r8_min ( aa, bb );
    b = r8_max ( aa, bb );
    alpha = a + b;
    beta = sqrt ( ( alpha - 2.0 ) / ( 2.0 * a * b - alpha ) );
    gamma = a + 1.0 / beta;

    for ( ; ; )
    {
      u1 = r8_uniform_01_sample ( );
      u2 = r8_uniform_01_sample ( );
      v = beta * log ( u1 / ( 1.0 - u1 ) );
      w = a * exp ( v );

      z = u1 * u1 * u2;
      r = gamma * v - log4;
      s = a + r - w;

      if ( 5.0 * z <= s + 1.0 + log5 )
      {
        break;
      }

      t = log ( z );
      if ( t <= s )
      {
        break;
      }

      if ( t <= ( r + alpha * log ( alpha / ( b + w ) ) ) )
      {
        break;
      }
    }
  }
//
//  Algorithm BC
//
  else
  {
    a = r8_max ( aa, bb );
    b = r8_min ( aa, bb );
    alpha = a + b;
    beta = 1.0 / b;
    delta = 1.0 + a - b;
    k1 = delta * ( 1.0 / 72.0 + b / 24.0 ) 
      / ( a / b - 7.0 / 9.0 );
    k2 = 0.25 + ( 0.5 + 0.25 / delta ) * b;

    for ( ; ; )
    {
      u1 = r8_uniform_01_sample ( );
      u2 = r8_uniform_01_sample ( );

      if ( u1 < 0.5 )
      {
        y = u1 * u2;
        z = u1 * y;

        if ( k1 <= 0.25 * u2 + z - y )
        {
          continue;
        }
      }
      else
      {
        z = u1 * u1 * u2;

        if ( z <= 0.25 )
        {
          v = beta * log ( u1 / ( 1.0 - u1 ) );
          w = a * exp ( v );

          if ( aa == a )
          {
            value = w / ( b + w );
          }
          else
          {
            value = b / ( b + w );
          }
          return value;
        }

        if ( k2 < z )
        {
          continue;
        }
      }

      v = beta * log ( u1 / ( 1.0 - u1 ) );
      w = a * exp ( v );

      if ( log ( z ) <= alpha * ( log ( alpha / ( b + w ) ) + v ) - log4 )
      {
        break;
      }
    }
  }

  if ( aa == a )
  {
    value = w / ( b + w );
  }
  else
  {
    value = b / ( b + w );
  }
  return value;
}
//****************************************************************************80

double r8_chi_pdf ( double df, double rval )

//****************************************************************************80
//
//  Purpose:
//
//    R8_CHI_PDF evaluates the PDF of a chi-squared distribution.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C  by John Burkardt.
//
//  Parameters:
//
//    Input, double DF, the degrees of freedom.
//    0.0 < DF.
//
//    Input, double RVAL, the point where the PDF is evaluated.
//
//    Output, double R8_CHI_PDF, the value of the PDF at RVAL.
//
{
  double temp1;
  double temp2;
  double value;

  if ( df <= 0.0 )
  {
    cerr << "\n";
    cerr << "R8_CHI_PDF - Fatal error!\n";
    cerr << "  Degrees of freedom must be positive.\n";
    exit ( 1 );
  }
      
  if ( rval <= 0.0 )
  {
    value = 0.0;
  }
  else
  {
    temp2 = df * 0.5;

    temp1 = ( temp2 - 1.0 ) * log ( rval ) - 0.5 * rval 
      - temp2 * log ( 2.0 ) - r8_gamma_log ( temp2 );

    value = exp ( temp1 );
  }
  return value;
}
//****************************************************************************80

double r8_chi_sample ( double df )

//****************************************************************************80
//
//  Purpose:
//
//    R8_CHI_SAMPLE generates a Chi-Square random deviate.
//
//  Discussion:
//
//    This procedure generates a random deviate from the chi square distribution
//    with DF degrees of freedom random variable.
//
//    The algorithm exploits the relation between chisquare and gamma.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN77 version by Barry Brown, James Lovato.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double DF, the degrees of freedom.
//    0.0 < DF.
//
//    Output, double R8_CHI_SAMPLE, a random deviate from the distribution.
//
{
  double arg1;
  double arg2;
  double value;

  if ( df <= 0.0 )
  {
    cerr << "\n";
    cerr << "R8_CHI_SAMPLE - Fatal error!\n";
    cerr << "  DF <= 0.\n";
    cerr << "  Value of DF: " << df << "\n";
    exit ( 1 );
  }

  arg1 = 1.0;
  arg2 = df / 2.0;

  value = 2.0 * r8_gamma_sample ( arg1, arg2 );

  return value;
}
//****************************************************************************80

double r8_epsilon ( )

//****************************************************************************80
//
//  Purpose:
//
//    R8_EPSILON returns the R8 roundoff unit.
//
//  Discussion:
//
//    The roundoff unit is a number R which is a power of 2 with the
//    property that, to the precision of the computer's arithmetic,
//      1 < 1 + R
//    but
//      1 = ( 1 + R / 2 )
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 September 2012
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double R8_EPSILON, the R8 round-off unit.
//
{
  static double value = 2.220446049250313E-016;

  return value;
}
//****************************************************************************80

double r8_exponential_pdf ( double beta, double rval )

//****************************************************************************80
//
//  Purpose:
//
//    R8_EXPONENTIAL_PDF evaluates the PDF of an exponential distribution.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double BETA, the scale value.
//    0.0 < BETA.
//
//    Input, double RVAL, the point where the PDF is evaluated.
//
//    Output, double R8_EXPONENTIAL_PDF, the value of the PDF at RVAL.
//
{
  double value;

  if ( beta <= 0.0 )
  {
    cerr << "\n";
    cerr << "R8_EXPONENTIAL_PDF - Fatal error!\n";
    cerr << "  BETA parameter must be positive.\n";
    exit ( 1 );
  }

  if ( rval < 0.0 )
  {
    value = 0.0;
  }
  else
  {
    value = exp ( - rval / beta ) / beta;
  }

  return value;
}
//****************************************************************************80

double r8_exponential_sample ( double lambda )

//****************************************************************************80
//
//  Purpose:
//
//    R8_EXPONENTIAL_SAMPLE samples the exponential PDF.
//
//  Discussion:
//
//    Note that the parameter LAMBDA is a multiplier.  In some formulations,
//    it is used as a divisor instead.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double LAMBDA, the parameter of the PDF.
//
//    Output, double R8_EXPONENTIAL_SAMPLE, a sample of the PDF.
//
{
  double r;
  double value;

  r = r8_uniform_01_sample ( );

  value = - log ( r ) * lambda;

  return value;
}
//****************************************************************************80

double r8_exponential_01_sample ( )

//****************************************************************************80
//
//  Purpose:
//
//    R8_EXPONENTIAL_01_SAMPLE samples the standard exponential PDF.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double R8_EXPONENTIAL_01_SAMPLE, a sample of the PDF.
//
{
  double r;
  double value;

  r = r8_uniform_01_sample ( );

  value = - log ( r );

  return value;
}
//****************************************************************************80

double r8_gamma_log ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    R8_GAMMA_LOG evaluates the logarithm of the gamma function.
//
//  Discussion:
//
//    This routine calculates the LOG(GAMMA) function for a positive real
//    argument X.  Computation is based on an algorithm outlined in
//    references 1 and 2.  The program uses rational functions that
//    theoretically approximate LOG(GAMMA) to at least 18 significant
//    decimal digits.  The approximation for X > 12 is from reference
//    3, while approximations for X < 12.0 are similar to those in
//    reference 1, but are unpublished.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    19 April 2013
//
//  Author:
//
//    Original FORTRAN77 version by William Cody, Laura Stoltz.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    William Cody, Kenneth Hillstrom,
//    Chebyshev Approximations for the Natural Logarithm of the
//    Gamma Function,
//    Mathematics of Computation,
//    Volume 21, Number 98, April 1967, pages 198-203.
//
//    Kenneth Hillstrom,
//    ANL/AMD Program ANLC366S, DGAMMA/DLGAMA,
//    May 1969.
//
//    John Hart, Ward Cheney, Charles Lawson, Hans Maehly,
//    Charles Mesztenyi, John Rice, Henry Thatcher,
//    Christoph Witzgall,
//    Computer Approximations,
//    Wiley, 1968,
//    LC: QA297.C64.
//
//  Parameters:
//
//    Input, double X, the argument of the function.
//
//    Output, double R8_GAMMA_LOG, the value of the function.
//
{
  double c[7] = {
    -1.910444077728E-03, 
     8.4171387781295E-04, 
    -5.952379913043012E-04,
     7.93650793500350248E-04, 
    -2.777777777777681622553E-03, 
     8.333333333333333331554247E-02, 
     5.7083835261E-03 };
  double corr;
  static double d1 = -5.772156649015328605195174E-01;
  static double d2 = 4.227843350984671393993777E-01;
  static double d4 = 1.791759469228055000094023;
  static double frtbig = 2.25E+76;
  int i;
  double p1[8] = {
    4.945235359296727046734888, 
    2.018112620856775083915565E+02, 
    2.290838373831346393026739E+03, 
    1.131967205903380828685045E+04, 
    2.855724635671635335736389E+04, 
    3.848496228443793359990269E+04, 
    2.637748787624195437963534E+04, 
    7.225813979700288197698961E+03 };
  double p2[8] = { 
    4.974607845568932035012064, 
    5.424138599891070494101986E+02, 
    1.550693864978364947665077E+04, 
    1.847932904445632425417223E+05, 
    1.088204769468828767498470E+06, 
    3.338152967987029735917223E+06, 
    5.106661678927352456275255E+06, 
    3.074109054850539556250927E+06 };
  double p4[8] = {
    1.474502166059939948905062E+04, 
    2.426813369486704502836312E+06, 
    1.214755574045093227939592E+08, 
    2.663432449630976949898078E+09, 
    2.940378956634553899906876E+10, 
    1.702665737765398868392998E+11, 
    4.926125793377430887588120E+11, 
    5.606251856223951465078242E+11 };
  double q1[8] = { 
    6.748212550303777196073036E+01, 
    1.113332393857199323513008E+03, 
    7.738757056935398733233834E+03, 
    2.763987074403340708898585E+04, 
    5.499310206226157329794414E+04, 
    6.161122180066002127833352E+04, 
    3.635127591501940507276287E+04, 
    8.785536302431013170870835E+03 };
  double q2[8] = { 
    1.830328399370592604055942E+02, 
    7.765049321445005871323047E+03, 
    1.331903827966074194402448E+05, 
    1.136705821321969608938755E+06, 
    5.267964117437946917577538E+06, 
    1.346701454311101692290052E+07, 
    1.782736530353274213975932E+07, 
    9.533095591844353613395747E+06 };
  double q4[8] = { 
    2.690530175870899333379843E+03, 
    6.393885654300092398984238E+05, 
    4.135599930241388052042842E+07, 
    1.120872109616147941376570E+09, 
    1.488613728678813811542398E+10, 
    1.016803586272438228077304E+11, 
    3.417476345507377132798597E+11, 
    4.463158187419713286462081E+11 };
  double res;
  static double sqrtpi = 0.9189385332046727417803297;
  static double xbig = 2.55E+305;
  double xden;
  static double xinf = 1.79E+308;
  double xm1;
  double xm2;
  double xm4;
  double xnum;
  double y;
  double ysq;

  y = x;

  if ( 0.0 < y && y <= xbig )
  {
    if ( y <= r8_epsilon ( ) )
    {
      res = - log ( y );
    }
//
//  EPS < X <= 1.5.
//
    else if ( y <= 1.5 )
    {
      if ( y < 0.6796875 )
      {
        corr = -log ( y );
        xm1 = y;
      }
      else
      {
        corr = 0.0;
        xm1 = ( y - 0.5 ) - 0.5;
      }

      if ( y <= 0.5 || 0.6796875 <= y )
      {
        xden = 1.0;
        xnum = 0.0;
        for ( i = 0; i < 8; i++ )
        {
          xnum = xnum * xm1 + p1[i];
          xden = xden * xm1 + q1[i];
        }
        res = corr + ( xm1 * ( d1 + xm1 * ( xnum / xden ) ) );
      }
      else
      {
        xm2 = ( y - 0.5 ) - 0.5;
        xden = 1.0;
        xnum = 0.0;
        for ( i = 0; i < 8; i++ )
        {
          xnum = xnum * xm2 + p2[i];
          xden = xden * xm2 + q2[i];
        }
        res = corr + xm2 * ( d2 + xm2 * ( xnum / xden ) );
      }
    }
//
//  1.5 < X <= 4.0.
//
    else if ( y <= 4.0 )
    {
      xm2 = y - 2.0;
      xden = 1.0;
      xnum = 0.0;
      for ( i = 0; i < 8; i++ )
      {
        xnum = xnum * xm2 + p2[i];
        xden = xden * xm2 + q2[i];
      }
      res = xm2 * ( d2 + xm2 * ( xnum / xden ) );
    }
//
//  4.0 < X <= 12.0.
//
    else if ( y <= 12.0 )
    {
      xm4 = y - 4.0;
      xden = -1.0;
      xnum = 0.0;
      for ( i = 0; i < 8; i++ )
      {
        xnum = xnum * xm4 + p4[i];
        xden = xden * xm4 + q4[i];
      }
      res = d4 + xm4 * ( xnum / xden );
    }
//
//  Evaluate for 12 <= argument.
//
    else
    {
      res = 0.0;

      if ( y <= frtbig )
      {
        res = c[6];
        ysq = y * y;
        for ( i = 0; i < 6; i++ )
        {
          res = res / ysq + c[i];
        }
      }
      res = res / y;
      corr = log ( y );
      res = res + sqrtpi - 0.5 * corr;
      res = res + y * ( corr - 1.0 );
    }
  }
//
//  Return for bad arguments.
//
  else
  {
    res = xinf;
  }

  return res;
}
//****************************************************************************80

double r8_gamma_pdf ( double beta, double alpha, double rval )

//****************************************************************************80
//
//  Purpose:
//
//    R8_GAMMA_PDF evaluates the PDF of a gamma distribution.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double BETA, the rate parameter.
//    0.0 < BETA.
//
//    Input, double ALPHA, the shape parameter.
//    0.0 < ALPHA.
//
//    Input, double RVAL, the point where the PDF is evaluated.
//
//    Output, double R8_GAMMA_PDF, the value of the PDF at RVAL.
//
{
  double temp;
  double value;

  if ( alpha <= 0.0 )
  {
    cerr << "\n";
    cerr << "R8_GAMMA_PDF - Fatal error!\n";
    cerr << "  Parameter ALPHA is not positive.\n";
    exit ( 1 );
  }

  if ( beta <= 0.0 )
  {
    cerr << "\n";
    cerr << "R8_GAMMA_PDF - Fatal error!\n";
    cerr << "  Parameter BETA is not positive.\n";
    exit ( 1 );
  }

  if ( rval <= 0.0 )
  {
    value = 0.0;
  }
  else
  {
    temp = alpha * log ( beta ) + ( alpha - 1.0 ) * log ( rval ) 
      - beta * rval - r8_gamma_log ( alpha );

    value = exp ( temp );
  }

  return value;
}
//****************************************************************************80

double r8_gamma_sample ( double a, double r )

//****************************************************************************80
//
//  Purpose:
//
//    R8_GAMMA_SAMPLE generates a Gamma random deviate.
//
//  Discussion:
//
//    This procedure generates random deviates from the gamma distribution whose
//    density is (A^R)/Gamma(R) * X^(R-1) * Exp(-A*X)
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN77 version by Barry Brown, James Lovato.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Joachim Ahrens, Ulrich Dieter,
//    Generating Gamma Variates by a Modified Rejection Technique,
//    Communications of the ACM,
//    Volume 25, Number 1, January 1982, pages 47-54.
//
//    Joachim Ahrens, Ulrich Dieter,
//    Computer Methods for Sampling from Gamma, Beta, Poisson and
//    Binomial Distributions,
//    Computing,
//    Volume 12, Number 3, September 1974, pages 223-246.
//
//  Parameters:
//
//    Input, double A, the location parameter.
//
//    Input, double R, the shape parameter.
//
//    Output, double R8_GAMMA_SAMPLE, a random deviate from the distribution.
//
{
  double value;

  value = r8_gamma_01_sample ( r ) / a;

  return value;
}
//****************************************************************************80

double r8_gamma_01_sample ( double a )

//****************************************************************************80
//
//  Purpose:
//
//    R8_GAMMA_01_SAMPLE samples the standard Gamma distribution.
//
//  Discussion:
//
//    This procedure corresponds to algorithm GD in the reference.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN77 version by Barry Brown, James Lovato.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Joachim Ahrens, Ulrich Dieter,
//    Generating Gamma Variates by a Modified Rejection Technique,
//    Communications of the ACM,
//    Volume 25, Number 1, January 1982, pages 47-54.
//
//  Parameters:
//
//    Input, double A, the parameter of the standard gamma
//    distribution.  0.0 < A < 1.0.
//
//    Output, double R8_GAMMA_01_SAMPLE, a random deviate from the distribution.
//
{
  double a1 =  0.3333333;
  double a2 = -0.2500030;
  double a3 =  0.2000062;
  double a4 = -0.1662921;
  double a5 =  0.1423657;
  double a6 = -0.1367177;
  double a7 =  0.1233795;
  double b;
  double c;
  double d;
  double e;
  double e1 = 1.0;
  double e2 = 0.4999897;
  double e3 = 0.1668290;
  double e4 = 0.0407753;
  double e5 = 0.0102930;
  double p;
  double q;
  double q0;
  double q1 =  0.04166669;
  double q2 =  0.02083148;
  double q3 =  0.00801191;
  double q4 =  0.00144121;
  double q5 = -0.00007388;
  double q6 =  0.00024511;
  double q7 =  0.00024240;
  double r;
  double s;
  double s2;
  double si;
  double sqrt32 = 5.6568542494923801952;
  double t;
  double u;
  double v;
  double value;
  double w;
  double x;

  if ( 1.0 <= a )
  {
    s2 = a - 0.5;
    s = sqrt ( s2 );
    d = sqrt32 - 12.0 * s;
//
//  Immediate acceptance.
//
    t = r8_normal_01_sample ( );
    x = s + 0.5 * t;
    value = x * x;

    if ( 0.0 <= t )
    {
      return value;
    }
//
//  Squeeze acceptance.
//
    u = r8_uniform_01_sample ( );
    if ( d * u <= t * t * t )
    {
      return value;
    }

    r = 1.0 / a;
    q0 = (((((( 
        q7   * r 
      + q6 ) * r 
      + q5 ) * r 
      + q4 ) * r
      + q3 ) * r
      + q2 ) * r
      + q1 ) * r;
//
//  Approximation depending on size of parameter A.
//
    if ( 13.022 < a )
    {
      b = 1.77;
      si = 0.75;
      c = 0.1515 / s;
    }
    else if ( 3.686 < a )
    {
      b = 1.654 + 0.0076 * s2;
      si = 1.68 / s + 0.275;
      c = 0.062 / s + 0.024;
    }
    else
    {
      b = 0.463 + s + 0.178 * s2;
      si = 1.235;
      c = 0.195 / s - 0.079 + 0.16 * s;
    }
//
//  Quotient test.
//
    if ( 0.0 < x )
    {
      v = 0.5 * t / s;

      if ( 0.25 < fabs ( v ) )
      {
        q = q0 - s * t + 0.25 * t * t + 2.0 * s2 * log ( 1.0 + v );
      }
      else
      {
        q = q0 + 0.5 * t * t * (((((( 
            a7   * v 
          + a6 ) * v 
          + a5 ) * v 
          + a4 ) * v 
          + a3 ) * v 
          + a2 ) * v 
          + a1 ) * v;
      }

      if ( log ( 1.0 - u ) <= q )
      {
        return value;
      }
    }

    for ( ; ; )
    {
      e = r8_exponential_01_sample ( );
      u = 2.0 * r8_uniform_01_sample ( ) - 1.0;
 
      if ( 0.0 <= u )
      {
        t = b + fabs ( si * e );
      }
      else
      {
        t = b - fabs ( si * e );
      }
//
//  Possible rejection.
//
      if ( t < -0.7187449 )
      {
        continue;
      }
//
//  Calculate V and quotient Q.
//
      v = 0.5 * t / s;

      if ( 0.25 < fabs ( v ) )
      {
        q = q0 - s * t + 0.25 * t * t + 2.0 * s2 * log ( 1.0 + v );
      }
      else
      {
        q = q0 + 0.5 * t * t * (((((( 
            a7   * v 
          + a6 ) * v 
          + a5 ) * v 
          + a4 ) * v 
          + a3 ) * v 
          + a2 ) * v 
          + a1 ) * v;
      }
//
//  Hat acceptance.
//
      if ( q <= 0.0 )
      {
        continue;
      }

      if ( 0.5 < q )
      {
        w = exp ( q ) - 1.0;
      }
      else
      {
        w = (((( 
            e5   * q 
          + e4 ) * q 
          + e3 ) * q 
          + e2 ) * q 
          + e1 ) * q;
      }
//
//  May have to sample again.
//
      if ( c * fabs ( u ) <= w * exp ( e - 0.5 * t * t ) )
      {
        break;
      }
    }

    x = s + 0.5 * t;
    value = x * x;
  }
//
//  Method for A < 1.
//
  else if ( a < 1.0 )
  {
    b = 1.0 + 0.3678794 * a;

    for ( ; ; )
    {
      p = b * r8_uniform_01_sample ( );

      if ( p < 1.0 )
      {
        value = exp ( log ( p ) / a );
        if ( value <= r8_exponential_01_sample ( ) )
        {
          break;
        }
      }
      else
      {
        value = - log ( ( b - p ) / a );
        if ( ( 1.0 - a ) * log ( value ) <= r8_exponential_01_sample ( ) )
        {
          break;
        }
      }
    }
  }
  return value;
}
//****************************************************************************80

double r8_invchi_pdf ( double df, double rval )

//****************************************************************************80
//
//  Purpose:
//
//    R8_INVCHI_PDF evaluates the PDF of an inverse chi-squared distribution.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    21 April 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double DF, the degrees of freedom.
//    0.0 < DF.
//
//    Input, double RVAL, the point where the PDF is evaluated.
//
//    Output, double R8_INVCHI_PDF, the value of the PDF at RVAL.
//
{
  double temp1;
  double temp2;
  double value; 

  if ( df <= 0.0 )
  {
    cerr << "\n";
    cerr << "R8_INVCHI_PDF - Fatal error!\n";
    cerr << "  Degrees of freedom must be positive.\n";
    exit ( 1 );
  }
      
  if ( rval <= 0.0 )
  {
    value = 0.0;
  }
  else
  {
    temp2 = df * 0.5;

    temp1 = - temp2 * log ( 2.0 ) - ( temp2 + 1.0 ) * log ( rval ) 
      - 0.5 / rval - r8_gamma_log ( temp2 );

    value = exp ( temp1 );
  }
  return value;
}
//****************************************************************************80

double r8_invgam_pdf ( double beta, double alpha, double rval )

//****************************************************************************80
//
//  Purpose:
//
//    R8_INVGAM_PDF evaluates the PDF of an inverse gamma distribution.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double BETA, the rate parameter.
//    0.0 < BETA.
//
//    Input, double ALPHA, the shape parameter.
//    0.0 < ALPHA.
//
//    Input, double RVAL, the point where the PDF is evaluated.
//
//    Output, double R8_INVGAM_PDF, the value of the PDF at RVAL.
//
{
  double temp;
  double value;

  if ( alpha <= 0.0 )
  {
    cerr << "\n";
    cerr << "R8_INVGAM_PDF - Fatal error!\n";
    cerr << "  Parameter ALPHA is not positive.\n";
    exit ( 1 );
  }

  if ( beta <= 0.0 )
  {
    cerr << "\n";
    cerr << "R8_INVGAM_PDF - Fatal error!\n";
    cerr << "  Parameter BETA is not positive.\n";
    exit ( 1 );
  }

  if ( rval <= 0.0 )
  {
    value = 0.0;
  }
  else
  {
    temp = alpha * log ( beta ) - ( alpha + 1.0 ) * log ( rval ) 
      - beta / rval - r8_gamma_log ( alpha );

    value = exp ( temp );
  }

  return value;
}
//****************************************************************************80

double r8_max ( double x, double y )

//****************************************************************************80
//
//  Purpose:
//
//    R8_MAX returns the maximum of two R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    18 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, Y, the quantities to compare.
//
//    Output, double R8_MAX, the maximum of X and Y.
//
{
  double value;

  if ( y < x )
  {
    value = x;
  }
  else
  {
    value = y;
  }
  return value;
}
//****************************************************************************80

double r8_min ( double x, double y )

//****************************************************************************80
//
//  Purpose:
//
//    R8_MIN returns the minimum of two R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    31 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, Y, the quantities to compare.
//
//    Output, double R8_MIN, the minimum of X and Y.
//
{
  double value;

  if ( y < x )
  {
    value = y;
  }
  else
  {
    value = x;
  }
  return value;
}
//****************************************************************************80

double r8_mnor_pdf ( int ndim, double mu[], double cov[], double logdet, 
  double rvals[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8_MNOR_PDF evaluates the PDF of a multivariate normal distribution.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, int NDIM, the spatial dimension.
//
//    Input, double MU[NDIM], the mean values.
//
//    Input, double COV[NDIM*NDIM], the covariance matrix.
//    (This matrix should be symmetric, positive definite.  It is
//    overwritten during the action of this double.)
//
//    Input, double LOGDET, the logarithm of the determinant
//    of the covariance matrix.
//
//    Input, double RVALS[NDIM], the point where the PDF is evaluated.
//
//    Output, double R8_MNOR_PDF, the value of the PDF at RVALS.
//
{
  double *covinv;
  int i;
  int j;              
  double pi = 3.141592653589793;
  double rsum;
  double value;
//
//  Invert the covariance matrix.
//
  covinv = new double[ndim*ndim];

  for ( i = 0; i < ndim; i++ )
  {
    for ( j = 0; j <= i; j++ )
    {
      covinv[j+i*ndim] = cov[j+i*ndim];
    }
    for ( j = i + 1; j < ndim; j++ )
    {
      covinv[j+i*ndim] = cov[i+j*ndim];
    }
  }
//
//  Calculate the PDF of the multivariate normal distribution.
//
  rsum = 0.0;
  for ( j = 0; j < ndim; j++ )
  {
    for ( i = 0; i < ndim; i++ )
    {
      rsum = rsum + ( rvals[i] - mu[i] ) * ( rvals[j] - mu[j] ) 
        * covinv[i+j*ndim];
    }
  }

  rsum = - 0.5 * ( double ) ( ndim ) * log ( 2.0 * pi ) 
    - 0.5 * logdet - 0.5 * rsum;

  value = exp ( rsum );

  delete [] covinv;

  return value;
}
//****************************************************************************80

double *r8vec_mnor_sample ( double parm[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_MNOR_SAMPLE generates a multivariate normal deviate.
//
//  Discussion:
//
//    The method is:
//    1) Generate P independent standard normal deviates - Ei ~ N(0,1)
//    2) Using Cholesky decomposition find A so that A'*A = COVM
//    3) A' * E + MEANV ~ N(MEANV,COVM)
//
//    Note that PARM contains information needed to generate the
//    deviates, and is set up by SETGMN.
//
//    PARM(1) contains the size of the deviates, P
//    PARM(2:P+1) contains the mean vector.
//    PARM(P+2:P*(P+3)/2+1) contains the upper half of the Cholesky
//    decomposition of the covariance matrix.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN77 version by Barry Brown, James Lovato.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double PARM[P*(P+3)/2+1], parameters set by SETGMN.
//
//    Output, double R8VEC_MNOR_SAMPLE[P], a random deviate from the
//    distribution.
//
{
  double ae;
  int i;
  int icount;
  int j;
  int p;
  double *work;
  double *x;

  p = ( int ) ( parm[0] );
//
//  Generate P independent normal deviates.
//
  work = new double[p];

  for ( i = 0; i < p; i++ )
  {
    work[i] = r8_normal_01_sample ( );
  }
//
//  Compute X = MEANV + A' * WORK
//
  x = new double[p];

  for ( i = 0; i < p; i++ )
  {
    icount = 0;
    ae = 0.0;
    for ( j = 0; j <= i; j++ )
    {
      icount = icount + j;
      ae = ae + parm[i+j*p-icount+p+1] * work[j];
    }
    x[i] = ae + parm[i+1];
  }

  delete [] work;

  return x;
}
//****************************************************************************80

double r8_normal_pdf ( double av, double sd, double rval )

//****************************************************************************80
//
//  Purpose:
//
//    R8_NORMAL_PDF evaluates the PDF of a normal distribution.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double AV, the mean value.
//
//    Input, double SD, the standard deviation.
//    0.0 < SD.
//
//    Input, double RVAL, the point where the PDF is evaluated.
//
//    Output, double R8_NORMAL_PDF, the value of the PDF at RVAL.
//
{
  double pi = 3.141592653589793;
  double rtemp;
  double value;

  if ( sd <= 0.0 )
  {
    cerr << "\n";
    cerr << "R8_NORMAL_PDF - Fatal error!\n";
    cerr << "  Standard deviation must be positive.\n";
    exit ( 1 );
  }

  rtemp = ( rval - av ) * ( rval - av ) * 0.5 / ( sd * sd );

  value = exp ( - rtemp ) / sd / sqrt ( 2.0 * pi );

  return value;
}
//****************************************************************************80

double r8_normal_sample ( double av, double sd )

//****************************************************************************80
//
//  Purpose:
//
//    R8_NORMAL_SAMPLE generates a normal random deviate.
//
//  Discussion:
//
//    This procedure generates a single random deviate from a normal 
//    distribution with mean AV, and standard deviation SD.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN77 version by Barry Brown, James Lovato.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Joachim Ahrens, Ulrich Dieter,
//    Extensions of Forsythe's Method for Random
//    Sampling from the Normal Distribution,
//    Mathematics of Computation,
//    Volume 27, Number 124, October 1973, page 927-937.
//
//  Parameters:
//
//    Input, double AV, the mean.
//
//    Input, double SD, the standard deviation.
//
//    Output, double R8_NORMAL_SAMPLE, a random deviate from the distribution.
//
{
  double value;

  value = sd * r8_normal_01_sample ( ) + av;

  return value;
}
//****************************************************************************80

double r8_normal_01_sample ( )

//****************************************************************************80
//
//  Purpose:
//
//    R8_NORMAL_01_SAMPLE samples the standard normal probability distribution.
//
//  Discussion:
//
//    The standard normal probability distribution function (PDF) has
//    mean 0 and standard deviation 1.
//
//    The Box-Muller method is used, which is efficient, but
//    generates two values at a time.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double R8_NORMAL_01_SAMPLE, a normally distributed random value.
//
{
  const double pi = 3.14159265358979323;
  double r1;
  double r2;
  static int used = 0;
  double x;
  static double y = 0.0;
//
//  If we've used an even number of values so far, generate two more, 
//  return one, and save one.
//
  if ( ( used % 2 )== 0 )
  {
    r1 = r8_uniform_01_sample ( );
    r2 = r8_uniform_01_sample ( );

    x = sqrt ( -2.0 * log ( r1 ) ) * cos ( 2.0 * pi * r2 );
    y = sqrt ( -2.0 * log ( r1 ) ) * sin ( 2.0 * pi * r2 );
  }
  else
  {
    x = y;
  }

  used = used + 1;

  return x;
}
//****************************************************************************80

int r8_round_i4 ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    R8_ROUND_I4 rounds an R8, returning an I4.
//
//  Example:
//
//        X         Value
//
//      1.3         1
//      1.4         1
//      1.5         1 or 2
//      1.6         2
//      0.0         0
//     -0.7        -1
//     -1.1        -1
//     -1.6        -2
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    25 March 2013
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, the value.
//
//    Output, int R8_ROUND_I4, the rounded value.
//
{
  int value;

  if ( x < 0.0 )
  {
    value = - floor ( - x + 0.5 );
  }
  else
  {
    value =   floor (   x + 0.5 );
  }

  return value;
}
//****************************************************************************80

double r8_scinvchi_pdf ( double df, double s, double rval )

//****************************************************************************80
//
//  Purpose:
//
//    R8_SCINVCHI_PDF: PDF for a scaled inverse chi-squared distribution.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double DF, the degrees of freedom.
//    0.0 < DF.
//
//    Input, double S, the scale factor.
//    0.0 < S.
//
//    Input, double RVAL, the point where the PDF is evaluated.
//
//    Output, double R8_SCINVCHI_PDF, the value of the PDF at RVAL.
//    inverse-chi-square distribution.
//
{
  double temp1;
  double temp2;
  double value;

  if ( df <= 0.0 )
  {
    cerr << "\n";
    cerr << "R8_SCINVCHI_PDF - Fatal error!\n";
    cerr << "  Degrees of freedom must be positive.\n";
    exit ( 1 );
  }

  if ( s <= 0.0 )
  {
    cerr << "\n";
    cerr << "R8_SCINVCHI_PDF - Fatal error!\n";
    cerr << "  Scale parameter must be positive.\n";
    exit ( 1 );
  }

  if ( rval <= 0.0 )
  {
    value = 0.0;
  }
  else
  {
    temp2 = df * 0.5;

    temp1 = temp2 * log ( temp2 ) + temp2 * log ( s ) 
      - ( temp2 * s / rval ) 
      - ( temp2 + 1.0 ) * log ( rval ) - r8_gamma_log ( temp2 );

    value = exp ( temp1 );
  }

  return value;
}
//****************************************************************************80

double r8_uniform_pdf ( double lower, double upper, double rval )

//****************************************************************************80
//
//  Purpose:
//
//    R8_UNIFORM_PDF evaluates the PDF of a uniform distribution.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double LOWER, UPPER, the lower and upper range limits.
//    LOWER < UPPER.
//
//    Input, double RVAL, the point where the PDF is evaluated.
//
//    Output, double R8_UNIFORM_PDF, the value of the PDF at RVAL.
//
{
  double value;

  if ( upper <= lower )
  {
    cerr << "\n";
    cerr << "R8_UNIFORM_PDF - Fatal error!\n";
    cerr << "  For uniform PDF, the lower limit must be\n";
    cerr << "  less than the upper limit\n";
    exit ( 1 );
  }

  if ( rval < lower )
  {
    value = 0.0;
  }
  else if ( rval <= upper )
  {
    value = 1.0 / ( upper - lower );
  }
  else
  {
    value = 0.0;
  }

  return value;
}
//****************************************************************************80

double r8_uniform_sample ( double low, double high )

//****************************************************************************80
//
//  Purpose:
//
//    R8_UNIFORM_SAMPLE generates a uniform random deviate.
//
//  Discussion:
//
//    This procedure generates a real deviate uniformly distributed between
//    LOW and HIGH.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN77 version by Barry Brown, James Lovato.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double LOW, HIGH, the lower and upper bounds.
//
//    Output, double R8_UNIFORM_SAMPLE, a random deviate from the distribution.
//
{
  double value;

  value = low + ( high - low ) * r8_uniform_01_sample ( );

  return value;
}
//****************************************************************************80

double r8_uniform_01_sample ( )

//****************************************************************************80
//
//  Purpose:
//
//    R8_UNIFORM_01_SAMPLE generates a uniform random deviate from [0,1].
//
//  Discussion:
//
//    This function should be the only way that the package accesses random
//    numbers.
//
//    Setting OPTION to 0 accesses the R8_UNIFORM_01() function in RNGLIB,
//    for which there are versions in various languages, which should result
//    in the same values being returned.
//
//    Setting OPTION to 1 in the C++ version calls the system
//    RNG "rand()".
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    John Burkardt.
//
//  Parameters:
//
//    Output, double R8_UNIFORM_01_SAMPLE, a random deviate.
//
{
  static int option = 0;
  double value;

  if ( option == 1 )
  {
    value = ( double ) rand ( ) / ( double ) RAND_MAX;
  }
  else
  {
    value = r8_uni_01 ( );
  }

  return value;
}
//****************************************************************************80

double *r8block_zero_new ( int l, int m, int n )

//****************************************************************************80
//
//  Purpose:
//
//    R8BLOCK_ZERO_NEW returns a new zeroed R8BLOCK.
//
//  Discussion:
//
//    An R8BLOCK is a triple dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    13 April 2013
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int L, M, N, the number of rows and columns.
//
//    Output, double R8BLOCK_ZERO_NEW[L*M*N], the new zeroed matrix.
//
{
  double *a;
  int i;
  int j;
  int k;

  a = new double[l*m*n];

  for ( k = 0; k < n; k++ )
  {
    for ( j = 0; j < m; j++ )
    {
      for ( i = 0; i < l; i++ )
      {
        a[i+j*l+k*l*m] = 0.0;
      }
    }
  }
  return a;
}
//****************************************************************************80

int r8mat_pofa ( double a[], int lda, int n )

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_POFA factors a real symmetric positive definite matrix.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN77 version by Jack Dongarra, Cleve Moler, Jim Bunch, 
//    Pete Stewart,
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
//    LINPACK User's Guide,
//    SIAM, (Society for Industrial and Applied Mathematics),
//    3600 University City Science Center,
//    Philadelphia, PA, 19104-2688.
//    ISBN 0-89871-172-X
//
//  Parameters:
//
//    Input/output, double A[LDA*N].  On input, the symmetric matrix
//    to be  factored.  Only the diagonal and upper triangle are used.
//    On output, an upper triangular matrix R so that A = R'*R
//    where R' is the transpose.  The strict lower triangle is unaltered.
//    If INFO /= 0, the factorization is not complete.
//
//    Input, int LDA, the leading dimension of the array A.
//
//    Input, int N, the order of the matrix.
//
//    Output, int R8MAT_POFA, error flag.
//    0, for normal return.
//    K, signals an error condition.  The leading minor of order K is not
//    positive definite.
//
{
  double dot;
  int i;
  int info;
  int j;
  int k;
  double s;
  double t;

  for ( j = 1; j <= n; j++ )
  {
    s = 0.0;

    for ( k = 1; k <= j-1; k++ )
    {
      dot = 0.0;
      for ( i = 0; i < k - 1; i++ )
      {
        dot = dot + a[i+(k-1)*lda] * a[i+(j-1)*lda];
      }
      t = a[k-1+(j-1)*lda] - dot;
      t = t / a[k-1+(k-1)*lda];
      a[k-1+(j-1)*lda] = t;
      s = s + t * t;
    }

    s = a[j-1+(j-1)*lda] - s;

    if ( s <= 0.0 )
    {
      info = j;
      return info;
    }

    a[j-1+(j-1)*lda] = sqrt ( s );
  }

  info = 0;

  return info;
}
//****************************************************************************80

double *r8mat_zero_new ( int m, int n )

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_ZERO_NEW returns a new zeroed R8MAT.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    03 October 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, N, the number of rows and columns.
//
//    Output, double R8MAT_ZERO_NEW[M*N], the new zeroed matrix.
//
{
  double *a;
  int i;
  int j;

  a = new double[m*n];

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      a[i+j*m] = 0.0;
    }
  }
  return a;
}
//****************************************************************************80

double *r8vec_copy_new ( int n, double a1[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_COPY_NEW copies an R8VEC to a new R8VEC.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    03 July 2008
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of entries in the vectors.
//
//    Input, double A1[N], the vector to be copied.
//
//    Output, double R8VEC_COPY_NEW[N], the copy of A1.
//
{
  double *a2;
  int i;

  a2 = new double[n];

  for ( i = 0; i < n; i++ )
  {
    a2[i] = a1[i];
  }
  return a2;
}
//****************************************************************************80

void r8vec_heap_d ( int n, double a[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_HEAP_D reorders an R8VEC into a descending heap.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//    A heap is an array A with the property that, for every index J,
//    A[J] >= A[2*J+1] and A[J] >= A[2*J+2], (as long as the indices
//    2*J+1 and 2*J+2 are legal).
//
//  Diagram:
//
//                  A(0)
//
//            A(1)         A(2)
//
//      A(3)       A(4)  A(5) A(6)
//
//    A(7) A(8)  A(9) A(10)
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    30 April 1999
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Albert Nijenhuis, Herbert Wilf,
//    Combinatorial Algorithms,
//    Academic Press, 1978, second edition,
//    ISBN 0-12-519260-6.
//
//  Parameters:
//
//    Input, int N, the size of the input array.
//
//    Input/output, double A[N].
//    On input, an unsorted array.
//    On output, the array has been reordered into a heap.
//
{
  int i;
  int ifree;
  double key;
  int m;
//
//  Only nodes (N/2)-1 down to 0 can be "parent" nodes.
//
  for ( i = (n/2)-1; 0 <= i; i-- )
  {
//
//  Copy the value out of the parent node.
//  Position IFREE is now "open".
//
    key = a[i];
    ifree = i;

    for ( ; ; )
    {
//
//  Positions 2*IFREE + 1 and 2*IFREE + 2 are the descendants of position
//  IFREE.  (One or both may not exist because they equal or exceed N.)
//
      m = 2 * ifree + 1;
//
//  Does the first position exist?
//
      if ( n <= m )
      {
        break;
      }
      else
      {
//
//  Does the second position exist?
//
        if ( m + 1 < n )
        {
//
//  If both positions exist, take the larger of the two values,
//  and update M if necessary.
//
          if ( a[m] < a[m+1] )
          {
            m = m + 1;
          }
        }
//
//  If the large descendant is larger than KEY, move it up,
//  and update IFREE, the location of the free position, and
//  consider the descendants of THIS position.
//
        if ( key < a[m] )
        {
          a[ifree] = a[m];
          ifree = m;
        }
        else
        {
          break;
        }
      }
    }
//
//  When you have stopped shifting items up, return the item you
//  pulled out back to the heap.
//
    a[ifree] = key;
  }

  return;
}
//****************************************************************************80

void r8vec_sort_heap_a ( int n, double a[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_SORT_HEAP_A ascending sorts an R8VEC using heap sort.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    30 April 1999
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Albert Nijenhuis, Herbert Wilf,
//    Combinatorial Algorithms,
//    Academic Press, 1978, second edition,
//    ISBN 0-12-519260-6.
//
//  Parameters:
//
//    Input, int N, the number of entries in the array.
//
//    Input/output, double A[N].
//    On input, the array to be sorted;
//    On output, the array has been sorted.
//
{
  int n1;
  double temp;

  if ( n <= 1 )
  {
    return;
  }
//
//  1: Put A into descending heap form.
//
  r8vec_heap_d ( n, a );
//
//  2: Sort A.
//
//  The largest object in the heap is in A[0].
//  Move it to position A[N-1].
//
  temp = a[0];
  a[0] = a[n-1];
  a[n-1] = temp;
//
//  Consider the diminished heap of size N1.
//
  for ( n1 = n-1; 2 <= n1; n1-- )
  {
//
//  Restore the heap structure of the initial N1 entries of A.
//
    r8vec_heap_d ( n1, a );
//
//  Take the largest object from A[0] and move it to A[N1-1].
//
    temp = a[0];
    a[0] = a[n1-1];
    a[n1-1] = temp;
  }

  return;
}
//****************************************************************************80

double r8vec_sum ( int n, double a[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_SUM returns the sum of an R8VEC.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    15 October 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of entries in the vector.
//
//    Input, double A[N], the vector.
//
//    Output, double R8VEC_SUM, the sum of the vector.
//
{
  int i;
  double value;

  value = 0.0;
  for ( i = 0; i < n; i++ )
  {
    value = value + a[i];
  }
  return value;
}
//****************************************************************************80

void r8vec_transpose_print ( int n, double a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_TRANSPOSE_PRINT prints an R8VEC "transposed".
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//  Example:
//
//    A = (/ 1.0, 2.1, 3.2, 4.3, 5.4, 6.5, 7.6, 8.7, 9.8, 10.9, 11.0 /)
//    TITLE = 'My vector:  '
//
//    My vector:
//
//        1.0    2.1    3.2    4.3    5.4
//        6.5    7.6    8.7    9.8   10.9
//       11.0
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    12 November 2010
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of components of the vector.
//
//    Input, double A[N], the vector to be printed.
//
//    Input, string TITLE, a title.
//
{
  int i;
  int ihi;
  int ilo;

  cout << "\n";
  cout << title << "\n";
  cout << "\n";

  if ( n <= 0 )
  {
    cout << "  (Empty)\n";
    return;
  }

  for ( ilo = 0; ilo < n; ilo = ilo + 5 )
  {
    ihi = i4_min ( ilo + 5, n );
    for ( i = ilo; i < ihi; i++ )
    {
      cout << "  " << setw(12) << a[i];
    }
    cout << "\n";
  }

  return;
}
//****************************************************************************80

double *r8vec_zero_new ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_ZERO_NEW creates and zeroes an R8VEC.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    10 July 2008
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of entries in the vector.
//
//    Output, double R8VEC_ZERO_NEW[N], a vector of zeroes.
//
{
  double *a;
  int i;

  a = new double[n];

  for ( i = 0; i < n; i++ )
  {
    a[i] = 0.0;
  }
  return a;
}
//****************************************************************************80

void restart_read ( int chain_num, double fit[], int gen_num, int par_num, 
  string restart_read_filename, double z[] )

//****************************************************************************80
//
//  Purpose:
//
//    RESTART_READ reads parameter sample data from a restart file.
//
//  Discussion:
//
//    Only a single generation (presumably the last one) was written to the file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 April 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, int CHAIN_NUM, the total number of chains.
//
//    Output, double FIT[CHAIN_NUM*GEN_NUM], the likelihood of
//    each sample.
//
//    Input, int GEN_NUM, the total number of generations.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input, string RESTART_READ_FILENAME, the name of 
//    the restart file.
//
//    Output, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain 
//    sample data.
//
{
  int chain_index;
  int gen_index = 0;
  int index;
  string line;
  int par_index;
  ifstream restart;

  restart.open ( restart_read_filename.c_str ( ) );

  if ( ! restart )
  {
    cerr << "\n";
    cerr << "RESTART_READ - Fatal error!\n";
    cerr << "  Could not open the file \"" 
         << restart_read_filename << "\".\n";
    exit ( 1 );
  }
//
//  Assume only one generation.
//
  gen_index = 0;
//
//  Read and ignore line 1.
//
  getline ( restart, line );
//
//  Read the final fitness and parameter values for each chain.
//
  for ( chain_index = 0; chain_index < chain_num; chain_index++ )
  {
    index = chain_index 
          + chain_num * gen_index;
    restart >> fit[index];
    for ( par_index = 0; par_index < par_num; par_index++ )
    {
      index = par_index 
            + par_num * chain_index 
            + par_num * chain_num * gen_index;
      restart >> z[index];
    }
  }

  restart.close ( );

  return;
}
//****************************************************************************80

void restart_write ( int chain_num, double fit[], int gen_num, int par_num, 
  string restart_write_filename, double z[] )

//****************************************************************************80
/*
  Purpose:

    RESTART_WRITE writes a restart file.

  Discussion:

    Only data for the final generation is written.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    23 April 2013

  Author:

    Original FORTRAN90 version by Guannan Zhang.
    C++ version by John Burkardt.

  Parameters:

    Input, int CHAIN_NUM, the total number of chains.

    Input, double FIT[CHAIN_NUM*GEN_NUM], the likelihood of
    each sample.

    Input, int GEN_NUM, the total number of generations.

    Input, int PAR_NUM, the total number of parameters.

    Input, string RESTART_WRITE_FILENAME, the name of the 
    restart file.

    Input, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain 
    sample data.
*/
{
  int c;
  int p;
  ofstream restart;

  restart.open ( restart_write_filename.c_str ( ) );   

  if ( !restart )
  {
    cerr << "\n";
    cerr << "RESTART_WRITE - Fatal error!\n";
    cerr << "  Could not open \"" << restart_write_filename << "\".\n";
    exit ( 1 );
  }

  restart << "DREAM.C:Parameter_values_for_restart.\n";

  for ( c = 0; c < chain_num; c++ )
  {
    restart << "  " << c
            << "  "<< fit[c+(gen_num-1)*chain_num];
    for ( p = 0; p < par_num; p++ )
    {
      restart << "  " << z[p+c*par_num+(gen_num-1)*par_num*chain_num];
    }
    restart << "\n"; 
  }

  restart.close ( );

  cout << "\n";
  cout << "RESTART_WRITE:\n";
  cout << "  Created restart file \"" << restart_write_filename << "\".\n";

  return;
}
//****************************************************************************80

double *sample_candidate ( int chain_index, int chain_num, double cr[], 
  int cr_index, int cr_num, int gen_index, int gen_num, 
  double jumprate_table[], int jumpstep, double limits[], int pair_num, 
  int par_num, double z[] )

//****************************************************************************80
//
//  Purpose:
//
//    SAMPLE_CANDIDATE generates candidate parameter samples.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Jasper Vrugt, CJF ter Braak, CGH Diks, Bruce Robinson, James Hyman, 
//    Dave Higdon,
//    Accelerating Markov Chain Monte Carlo Simulation by Differential 
//    Evolution with Self-Adaptive Randomized Subspace Sampling,
//    International Journal of Nonlinear Sciences and Numerical Simulation,
//    Volume 10, Number 3, March 2009, pages 271-288.
//
//  Parameters:
//
//    Input, int CHAIN_INDEX, the chain index.
//    0 <= CHAIN_INDEX < CHAIN_NUM.
//
//    Input, int CHAIN_NUM, the total number of chains.
//
//    Input, double CR[CR_NUM], the CR values.
//
//    Input, int CR_INDEX, the index of the chosen CR value.
//    0 <= CR_INDEX < CR_NUM.
//
//    Input, int CR_NUM, the total number of CR values.
//
//    Input, int GEN_INDEX, the current generation.
//    0 <= GEN_INDEX < GEN_NUM.
//
//    Input, int GEN_NUM, the total number of generations.
//
//    Input, double JUMPRATE_TABLE[PAR_NUM], the jumprate table.
//
//    Input, int JUMPSTEP, forces a "long jump" every
//    JUMPSTEP generations.
//
//    Input, double LIMITS[2*PAR_NUM], limits for the parameters.
//
//    Input, int PAIR_NUM, the number of pairs of 
//    crossover chains.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain 
//    sample data.
//
//    Output, double SAMPLE_CANDIDATE[PAR_NUM], a candidate parameter sample.
//
//  Local parameters:
//
//    Input, int JUMP_DIM[JUMP_NUM], the dimensions in which
//    a jump is to be made.
//
//    Local, int JUMP_NUM, the number of dimensions in which
//    a jump will be made.  0 <= JUMP_NUM <= PAR_NUM.
//
//    Local, double JUMPRATE, the jump rate.
//
{
  double av;
  double b;
  double *diff;
  double *eps;
  int i;
  int *jump_dim;
  int jump_num;
  double jumprate;
  double *noise_e;
  int pair[2];
  int *r;
  double r2;
  double sd;
  double *zp;
//
//  Used to calculate E following a uniform distribution on (-B,+B).
//  Because B is currently zero, the noise term is suppressed.
//
  b = 0.0;
//
//  Pick N pairs of other chains for crossover.
//
  r = new int[2*pair_num];

  for ( i = 0; i < pair_num; i++ )
  {
    while ( 1 )
    {
      r2 = r8_uniform_01_sample ( );
      pair[0] = ( int ) ( r2 * ( double ) chain_num );
      r2 = r8_uniform_01_sample ( );
      pair[1] = ( int ) ( r2 * ( double ) chain_num );

      if ( pair[0] != pair[1] &&
           pair[0] != chain_index && 
           pair[1] != chain_index )
      {
        break;
      }
    }
    r[0+i*2] = pair[0];
    r[1+i*2] = pair[1];
  }
//
//  Determine the jump rate.
//
  jump_dim = new int[par_num];

  jumprate_choose ( cr, cr_index, cr_num, gen_index, jump_dim, jump_num, 
    jumprate, jumprate_table, jumpstep, par_num );
//
//  Calculate E in equation 4 of Vrugt.
//
  noise_e = new double[par_num];

  for ( i = 0; i < par_num; i++ )
  {
    noise_e[i] = b * ( 2.0 * r8_uniform_01_sample ( ) - 1.0 );
  }
//
//  Get epsilon value from multinormal distribution                      
//
  eps = new double[par_num];

  av = 0.0;
  sd = 1.0E-10;
  for ( i = 0; i < par_num; i++ )
  {
    eps[i] = r8_normal_sample ( av, sd );
  }
//
//  Generate the candidate sample ZP based on equation 4 of Vrugt.
//
  diff = diff_compute ( chain_num, gen_index, gen_num, jump_dim, jump_num, 
    pair_num, par_num, r, z );

  zp = new double[par_num];

  for ( i = 0; i < par_num; i++ )
  {
    zp[i] = z[i+chain_index*par_num+(gen_index-1)*par_num*chain_num];
  }
  for ( i = 0; i < par_num; i++ )
  {
    zp[i] = zp[i] + ( 1.0 + noise_e[i] ) * jumprate * diff[i] + eps[i];
  }
//
//  Enforce limits on the sample ZP.
//
  sample_limits ( limits, par_num, zp );

  delete [] diff;
  delete [] eps;
  delete [] jump_dim;
  delete [] noise_e;
  delete [] r;

  return zp;
}
//****************************************************************************80

void sample_limits ( double limits[], int par_num, double zp[] )

//****************************************************************************80
//
//  Purpose:
//
//    SAMPLE_LIMITS enforces limits on a sample variable.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double LIMITS[2*PAR_NUM], the parameter limits.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input/output, double ZP[PAR_NUM], a variable, whose entries,
//    if necessary, will be "folded" so that they lie within the limits.
//
{
  int i;
  double w;

  for ( i = 0; i < par_num; i++ )
  {
    w = limits[1+i*2] - limits[0+i*2];

    if ( w == 0.0 )
    {
      zp[i] = limits[0+i*2];
    }
    else if ( w < 0.0 )
    {
      cerr << "\n";
      cerr << "SAMPLE_LIMITS - Fatal error!\n";
      cerr << "  Upper limit less than lower limit.\n";
      exit ( 1 );
    }
    else
    {
      while ( zp[i] < limits[0+i*2] )
      {
        zp[i] = zp[i] + w;
      }
      while ( limits[1+i*2] < zp[i] )
      {
        zp[i] = zp[i] - w;
      }
    }
  }
  return;
}

//****************************************************************************80

double *std_compute ( int chain_num, int gen_index, int gen_num, int par_num, 
  double z[] )

//****************************************************************************80
//
//  Purpose:
//
//    STD_COMPUTE computes the current standard deviations, for each parameter.
//
//  Discussion:
//
//    The computation encompasses all chains and generations up to the
//    current ones.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 May 2013
//
//  Author:
//
//    Original FORTRAN90 version by Guannan Zhang.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, int CHAIN_NUM, the total number of chains.
//
//    Input, int GEN_INDEX, the current generation.
//    0 <= GEN_INDEX < GEN_NUM.
//
//    Input, int GEN_NUM, the total number of generations.
//
//    Input, int PAR_NUM, the total number of parameters.
//
//    Input, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain 
//    sample data.
//
//    Output, double STD_COMPUTE[PAR_NUM], the standard deviations.
//
{
  int i;
  int j;
  int k;
  double mean;
  double *std;

  std = new double[par_num];

  for ( i = 0; i < par_num; i++ )
  {
    mean = 0.0;
    for ( k = 0; k <= gen_index; k++ )
    {
      for ( j = 0; j < chain_num; j++ )
      {
        mean = mean + z[i+j*par_num+k*par_num*chain_num];
      }
    }
    mean = mean / ( double ) ( chain_num ) / ( double ) ( gen_index );

    std[i] = 0.0;
    for ( k = 0; k <= gen_index; k++ )
    {
      for ( j = 0; j < chain_num; j++ )
      {
        std[i] = std[i] + pow ( z[i+j*par_num+k*par_num*chain_num] - mean, 2 );
      }
    }
    std[i] = sqrt ( std[i] / ( double ) ( chain_num * gen_index - 1 ) );
  }

  return std;
}