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

using namespace std;

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

int main ( );
double sample_likelihood1 ( int par_num, double zp[] );
double sample_likelihood2 ( int par_num, double zp[] );

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

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    MAIN is the main program for DREAM1_TEST.
//
//  Discussion:
//
//    DREAM1_TEST tests the DREAM1 library.
//
//  Working notes:
//
//    1) the multivariate normal option has not been tested, either on its
//       own, or as part of a more complicated prior.  In particular, the
//       PARM array is set up, but does not seem to be used.
//
//  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.
//
//  Local parameters:
//
//    Local, string CHAIN_FILENAME, the "base" filename
//    to be used for the chain files.  This name should include a string
//    of 0's which will be replaced by the chain indices.  For example,
//    "chain000.txt" would work.
//
//    Local, int CHAIN_NUM, the total number of chains.
//
//    Local, double COV_MAT[MNOR_NUM*MNOR_NUM], the covariance
//    matrix for a multivariate normal density.
//
//    Local, int CR_NUM, the total number of CR values.
//
//    Local, double FIT[CHAIN_NUM*GEN_NUM], the likelihood of
//    each sample.
//
//    Local, int GEN_NUM, the total number of generations.
//
//    Local, double GR[PAR_NUM*GR_NUM], 
//    the Gelman-Rubin R statistic.
//
//    Local, logical GR_CONV, the Gelman-Rubin convergence flag.
//
//    Local, int GR_COUNT, counts the number of generations
//    at which the Gelman-Rubin statistic has been computed.
//
//    Local, string GR_FILENAME, the name of the file
//    in which values of the Gelman-Rubin statistic will be recorded.
//
//    Local, int GR_NUM, the number of times the Gelman-Rubin
//    statistic may be computed.
//
//    Local, double GR_THRESHOLD, the convergence tolerance for
//    the Gelman-Rubin statistic.
//
//    Local, string INPUT_FILENAME, the name of the file
//    containing input data.
//
//    Local, double JUMPRATE_TABLE[PAR_NUM], the jumprate table.
//
//    Local, int JUMPSTEP, forces a "long jump" every
//    JUMPSTEP generations.
//
//    Local, double LIMITS[2*PAR_NUM], lower and upper bounds
//    for each parameter.
//
//    Local, double LOGDET, the logarithm of the determinant
//    of the covariance matrix.
//
//    Local, string MNOR_FILENAME, the name of the file
//    containing multivariate normal data.
//
//    Local, int MNOR_IND[MNOR_NUM], indices of parameters
//    associated with multivariate normal distributions.
//
//    Local, int MNOR_MEAN[MNOR_NUM], the means associated
//    with the multivariate density.
//
//    Local, int MNOR_NUM, the number of parameters with
//    a multivariate normal distribution.
//
//    Local, string OB_FILENAME, the name of the file
//    containing observational data.
//
//    Local, int OB_NUM, the number of observations.
//
//    Local, int PAIR_NUM, the number of pairs of 
//    crossover chains.
//
//    Local, int PAR_NUM, the total number of parameters.
//
//    Local, double PARM[MNOR_NUM*(MNOR_NUM+3)/2+1], parameters 
//    for a multivariate normal distribution.
//
//    Local, 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.
//
//    Local, double PRIOR_DEN_PAR_VAL[2*PAR_NUM], parameter values
//    associated with prior densities.
//
//    Local, int PRIOR_DEN_TYPE[PAR_NUM], the type
//    of the prior density.
//    1, beta
//    2, chi square
//    3, exponential
//    4, gamma
//    5, inverse chi square
//    6, inverse gamma
//    7, normal
//    8, scaled inverse chi square
//    9, uniform
//    10, multinormal
//
//    Local, string PRIOR_FILENAME, the name of the file
//    containing information about the prior density functions.
//
//    Local, string RESTART_READ_FILENAME, the name of the file
//    containing restart information.
//
//    Local, string RESTART_WRITE_FILENAME, the name of the file
//    to be written, containing restart information.
//
//    Local, bool RESTART_FLAG, is TRUE if the computation is to be initialized
//    from a restart file.
//
//    Local, string RESTART_STRING, is 'Y' or 'y' if the
//    computation is to be initialized from a restart file.
//
//    Local, int UNI_IND[UNI_NUM], indices of parameters
//    associated with univariate distributions.
//
//    Local, int UNI_NUM, the number of parameters with
//    a univariate distribution.
//
//    Local, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain 
//    sample data.
//
{
  string chain_filename = "chain000.txt";
  int chain_num;
  double *cov_mat;
  int cr_num;
  double *fit;
  int gen_num;
  double *gr;
  int gr_conv;
  int gr_count;
  string gr_filename = "gr.txt";
  int gr_num;
  double gr_threshold;
  string input_filename = "input.txt";
  double *jumprate_table;
  int jumpstep;
  double *limits;
  double logdet;
  string mnor_filename = "covmatrix.txt";
  int *mnor_ind;
  double *mnor_mean;
  int mnor_num;
  string ob_filename;
  int ob_num;
  int pair_num;
  int par_num;
  double *parm;
  int printstep;
  int *prior_den_par_num;
  double *prior_den_par_val;
  int *prior_den_type;
  string prior_filename = "prior.txt";
  string restart_read_filename = "restart_read.txt";
  string restart_write_filename = "restart_write.txt";
  bool restart_flag;
  string restart_string;
  int *uni_ind;
  int uni_num;
  double *z;

  timestamp ( );
  cout << "\n";
  cout << "DREAM1_TEST\n";
  cout << "  C++ version\n";
  cout << "  Test the DREAM1 library.\n";
//
//  Get the parameter size from the input file.
//
  input_size ( input_filename, par_num );
//
//  Read the data from the input file.
//
  limits = r8mat_zero_new ( 2, par_num );

  input_read ( chain_num, cr_num, input_filename, 
    gr_threshold, jumpstep, limits, gen_num, 
    &ob_filename, ob_num, pair_num, par_num, printstep, 
    &restart_string );
//
//  Convert string to logical.
//
  if ( restart_string[0] == 'y' || restart_string[0] == 'Y' )
  {
    restart_flag = true;
  }
  else
  {
    restart_flag = false;
  }
//
//  Print the data as a job record.
//
  input_print ( chain_num, cr_num, input_filename, gr_threshold, 
    jumpstep, limits, gen_num, ob_filename, ob_num, pair_num, par_num, 
    printstep, restart_string );
//
//  Allocate and zero out memory.
//
  gr_num = gen_num / printstep;

  fit = r8mat_zero_new ( chain_num, gen_num );
  gr = r8mat_zero_new ( par_num, gr_num );
  mnor_ind = i4vec_zero_new ( par_num );
  prior_den_par_num = i4vec_zero_new ( par_num );
  prior_den_par_val = r8mat_zero_new ( 2, par_num );
  prior_den_type = i4vec_zero_new ( par_num );
  uni_ind = i4vec_zero_new ( par_num );
  z = r8block_zero_new ( par_num, chain_num, gen_num );
//
//  Define the prior distributions.
//
  prior_read ( mnor_ind, mnor_num, par_num, prior_den_par_num, 
    prior_den_par_val, prior_den_type, prior_filename, uni_ind, uni_num );

  prior_print ( mnor_ind, mnor_num, par_num, prior_den_par_num, 
    prior_den_par_val, prior_den_type, prior_filename, uni_ind, uni_num );
//
//  Define the multivariate prior distribution.
//
  if ( 1 == mnor_num )
  {
    cerr << "\n";
    cerr << "DREAM1_TEST - Fatal error!\n";
    cerr << "  Only one parameter follows multinormal distribution.\n";
    cerr << "  Please use the normal distribution instead.\n";
    exit ( 1 );
  }
  else if ( 1 < mnor_num )
  {
    cov_mat = r8mat_zero_new ( mnor_num, mnor_num );
    mnor_ind = i4vec_zero_new ( mnor_num );
    mnor_mean = r8vec_zero_new ( mnor_num );
    parm = r8vec_zero_new ( mnor_num*(mnor_num+1)/2+1 );

    mnor_read ( cov_mat, mnor_filename, mnor_num );

    mnor_init ( cov_mat, logdet, mnor_ind, mnor_mean, 
      mnor_num, par_num, parm, prior_den_par_val );
  }
//
//  Set the jump rate table.
//
  jumprate_table = jumprate_table_init ( pair_num, par_num );

  jumprate_table_print ( jumprate_table, pair_num, par_num );
//
//  Initialize the Gelman-Rubin data.
//
  gr_init ( gr, gr_conv, gr_count, gr_num, par_num );

  cout << "\n";
  cout << "GR_PRINT:\n";
  cout << "  GR_CONV  = " << gr_conv << "\n";
  cout << "  GR_COUNT = " << gr_count << "\n";
  cout << "  GR_NUM   = " << gr_num << "\n";
//
//  Set the first generation of the chains from restart data, or by sampling.
//
  if ( restart_flag )
  {
    restart_read ( chain_num, fit, gen_num, par_num, restart_read_filename, z );
  }
  else
  {
    chain_init ( chain_num, fit, gen_num, mnor_ind, mnor_num, par_num, 
      parm, prior_den_par_num, prior_den_par_val, prior_den_type, 
      sample_likelihood2, uni_ind, uni_num, z );
  }
  chain_init_print ( chain_num, fit, gen_num, par_num, restart_flag, 
    restart_read_filename, z );
//
//  Carry out the DREAM algorithm.
//
  dream_algm ( chain_num, cov_mat, cr_num, fit, gen_num, gr, 
    gr_conv, gr_count, gr_num, gr_threshold, jumprate_table, 
    jumpstep, limits, logdet, mnor_ind, mnor_mean, mnor_num, pair_num, 
    par_num, printstep, prior_den_par_num, prior_den_par_val, prior_den_type, 
    sample_likelihood2, uni_ind, uni_num, z );
//
//  Save Gelman-Rubin statistic values to a file.
//
  gr_write ( gr, gr_filename, gr_num, par_num, printstep );
//
//  Save parameter values for all chains at last generation.
//
  restart_write ( chain_num, fit, gen_num, par_num, restart_write_filename, z );
//
//  Write each chain to a separate file.
//
  chain_write ( chain_filename, chain_num, fit, gen_num, par_num, z );
//
//  Free memory.
//
  if ( 0 < mnor_num )
  {
    delete [] cov_mat;
  }
  delete [] fit;
  delete [] gr;
  delete [] jumprate_table;
  delete [] limits;
  if ( 0 < mnor_num )
  {
    delete [] mnor_ind;
  }
  if ( 0 < mnor_num )
  {
    delete [] mnor_mean;
  }
  if ( 0 < mnor_num )
  {
    delete [] parm;
  }
  delete [] prior_den_par_num;
  delete [] prior_den_par_val;
  delete [] prior_den_type;
  delete [] uni_ind;
  delete [] z;
//
//  Terminate.
//
  cout << "\n";
  cout << "DREAM1_TEST\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

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

double sample_likelihood1 ( int par_num, double zp[] )

//****************************************************************************80
//
//  Purpose:
//
//    SAMPLE_LIKELIHOOD1 computes the log likelihood function.
//
//  Discussion:
//
//    This is a 10D bimodal 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 PAR_NUM, the total number of parameters.
//
//    Input, double ZP[PAR_NUM], a sample.
//
//    Output, double SAMPLE_LIKELIHOOD1, the log likelihood function 
//    for the sample.
//
{
  double arg;
  double arg1;
  double arg2;
  int i;
  const double pi = 3.141592653589793;
  double value;

  arg1 = 0.0;
  arg2 = 0.0;
  for ( i = 0; i < par_num; i++ )
  {
    arg1 = arg1 + pow ( zp[i] + 5.0, 2 );
    arg2 = arg2 + pow ( zp[i] - 5.0, 2 );
  }
  arg = 
      1.0 / 3.0 / ( 2.0 * pi ) * exp ( -0.5 * arg1 ) 
    + 2.0 / 3.0 / ( 2.0 * pi ) * exp ( -0.5 * arg2 );

  value = log ( arg );

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

double sample_likelihood2 ( int par_num, double zp[] )

//****************************************************************************80
//
//  Purpose:
//
//    SAMPLE_LIKELIHOOD2 computes the log likelihood function.
//
//  Discussion:
//
//    This is a one mode Gaussian.
//
//  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 PAR_NUM, the total number of parameters.
//
//    Input, double ZP[PAR_NUM], a sample.
//
//    Output, double SAMPLE_LIKELIHOOD2, the log likelihood function 
//    for the sample.
//
{
  double arg1;
  double arg2;
  int i;
  const double pi = 3.141592653589793;
  double value;

  arg1 = 0.0;
  arg2 = 0.0;
  for ( i = 0; i < par_num; i++ )
  {
    arg1 = arg1 + pow ( zp[i] + 5.0, 2 );
    arg2 = arg2 + pow ( zp[i] - 5.0, 2 );
  }

  arg1 = log ( 1.0 / 3.0 / ( 2.0 * pi ) ) - 0.5 * arg1;
  arg2 = log ( 2.0 / 3.0 / ( 2.0 * pi ) ) - 0.5 * arg2;

  value = r8_max ( arg1, arg2 ); 

  return value;
}