# include <stdlib.h>
# include <stdio.h>
# include <math.h>
# include <time.h>
# include <string.h>

# include "dream1.h"
# include "rnglib.h"

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

/******************************************************************************/

int main ( )

/******************************************************************************/
/*
  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, char *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, char *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, char *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, char *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, char *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, char *PRIOR_FILENAME, the name of the file
    containing information about the prior density functions.

    Local, char *RESTART_READ_FILENAME, the name of the file
    containing restart information.

    Local, char *RESTART_WRITE_FILENAME, the name of the file
    to be written, containing restart information.

    Local, int RESTART_FLAG, is TRUE if the computation is to be initialized
    from a restart file.

    Local, char *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.
*/
{
  char *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;
  char *gr_filename = "gr.txt";
  int gr_num;
  double gr_threshold;
  char *input_filename = "input.txt";
  double *jumprate_table;
  int jumpstep;
  double *limits;
  double logdet;
  char *mnor_filename = "covmatrix.txt";
  int *mnor_ind;
  double *mnor_mean;
  int mnor_num;
  char ob_filename[255];
  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;
  char *prior_filename = "prior.txt";
  char *restart_read_filename = "restart_read.txt";
  char *restart_write_filename = "restart_write.txt";
  int restart_flag;
  char restart_string[10];
  int *uni_ind;
  int uni_num;
  double *z;

  timestamp ( );
  printf ( "\n" );
  printf ( "DREAM1_TEST\n" );
  printf ( "  C version\n" );
  printf ( "  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 = 1;
  }
  else
  {
    restart_flag = 0;
  }
/*
  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 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "DREAM1_TEST - Fatal error!\n" );
    fprintf ( stderr, "  Only one parameter follows multinormal distribution.\n" );
    fprintf ( stderr, "  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 );

  printf ( "\n" );
  printf ( "GR_PRINT:\n" );
  printf ( "  GR_CONV  = %d\n", gr_conv );
  printf ( "  GR_COUNT = %d\n", gr_count );
  printf ( "  GR_NUM   = %d\n", gr_num );
/*
  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 )
  {
    free ( cov_mat );
  }
  free ( fit );
  free ( gr );
  free ( jumprate_table );
  free ( limits );
  if ( 0 < mnor_num )
  {
    free ( mnor_ind );
  }
  if ( 0 < mnor_num )
  {
    free ( mnor_mean );
  }
  if ( 0 < mnor_num )
  {
    free ( parm );
  }
  free ( prior_den_par_num );
  free ( prior_den_par_val );
  free ( prior_den_type );
  free ( uni_ind );
  free ( z );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "DREAM1_TEST\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

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

double sample_likelihood1 ( int par_num, double zp[] )

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

    SAMPLE_LIKELIHOOD1 computes the log likelihood function.

  Discussion:

    This function corresponds to 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 int k = 2;
  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;
}
/******************************************************************************/

double sample_likelihood2 ( int par_num, double zp[] )

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

    SAMPLE_LIKELIHOOD2 computes the log likelihood function.

  Discussion:

    This function corresponds to 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;
}