function dream ( problem_size, problem_value, problem_prior_density, ...
  problem_prior_sample, problem_sample_likelihood )

%*****************************************************************************80
%
%% dream() is the main program for the DREAM algorithm.
%
%  Discussion:
%
%    The DREAM program was originally developed by Guannan Zhang, of
%    Oak Ridge National Laboratory (ORNL); it has been incorporated into 
%    the DAKOTA package of Sandia National Laboratory, and is
%    intended to form part of the ORNL package known as TASMANIA.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 January 2018
%
%  Author:
%
%    Original FORTRAN90 version by Guannan Zhang.
%    MATLAB 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.
%
%  Input:
%
%    function [ chain_num, cr_num, gen_num, pair_num, par_num ] ...
%    = problem_size (  )
%
%    function [ chain_filename, gr_filename, gr_threshold, jumpstep, ...
%    limits, printstep, restart_read_filename, restart_write_filename ] ...
%    = problem_value ( par_num )
%
%    function value = problem_prior_density ( par_num, zp )
%
%    function zp = problem_prior_sample ( par_num )
%
%    function value = problem_sample_likelihood ( par_num, zp )
%
%  Local:
%
%    string CHAIN_FILENAME, the "base" filename
%    to be used for the chain files.  If this is the empty string '',
%    then the chain files will not be written.  This name should 
%    include a string of 0s which will be replaced by the chain 
%    indices.  For example, "chain000.txt" would work as long as the
%    number of chains was 1000 or less.
%
%    integer CHAIN_NUM, the total number of chains.
%    3 <= CHAIN_NUM.
%
%    integer CR_NUM, the total number of CR values.
%    1 <= CR_NUM.
%
%    real FIT(CHAIN_NUM,GEN_NUM), the likelihood of
%    each sample.
%
%    integer GEN_NUM, the total number of generations.
%    2 <= GEN_NUM.
%
%    real GR(PAR_NUM,GR_NUM), 
%    the Gelman-Rubin R statistic.
%
%    logical GR_CONV, the Gelman-Rubin convergence flag.
%
%    integer GR_COUNT, counts the number of generations
%    at which the Gelman-Rubin statistic has been computed.
%
%    string GR_FILENAME, the name of the file
%    in which values of the Gelman-Rubin statistic will be recorded,
%    or '' if no such file is to be created.
%
%    integer GR_NUM, the number of times the Gelman-Rubin
%    statistic may be computed.
%
%    real GR_THRESHOLD, the convergence tolerance for
%    the Gelman-Rubin statistic.
%
%    real JUMPRATE_TABLE(PAR_NUM), the jumprate table.
%
%    integer JUMPSTEP, forces a "long jump" every
%    JUMPSTEP generations.
%
%    real LIMITS(2,PAR_NUM), lower and upper bounds
%    for each parameter.
%
%    integer PAIR_NUM, the number of pairs of 
%    crossover chains.
%    0 <= PAIR_NUM.
%
%    integer PAR_NUM, the total number of parameters.
%    1 <= PAR_NUM.
%
%    integer PRINTSTEP, the interval between generations on 
%    which the Gelman-Rubin statistic will be computed and written to a file.
%
%    string RESTART_READ_FILENAME, the name of the file
%    containing restart information, or '' if this is not a restart run.
%
%    string RESTART_WRITE_FILENAME, the name of the file
%    to be written, containing restart information, or '' if a restart file 
%    is not to be written.
%
%    real Z(PAR_NUM,CHAIN_NUM,GEN_NUM), the Markov chain 
%    sample data.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'dream():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  MCMC acceleration by Differential Evolution.\n' );
%
%  Get the problem sizes.
%
  [ chain_num, cr_num, gen_num, pair_num, par_num ] = problem_size ( );
%
%  Decide if the problem sizes are acceptable.
%
  if ( chain_num < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DREAM - Fatal error!\n' );
    fprintf ( 1, '  CHAIN_NUM < 3.\n' );
    return
  end

  if ( cr_num < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DREAM - Fatal error!\n' );
    fprintf ( 1, '  CR_NUM < 1.\n' );
    return
  end

  if ( gen_num < 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DREAM - Fatal error!\n' );
    fprintf ( 1, '  GEN_NUM < 2.\n' );
    return
  end

  if ( pair_num < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DREAM - Fatal error!\n' );
    fprintf ( 1, '  PAIR_NUM < 0.\n' );
    return
  end

  if ( par_num < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DREAM - Fatal error!\n' );
    fprintf ( 1, '  PAR_NUM < 1.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'DREAM:\n' );
  fprintf ( 1, '  Accepted data from problem_size():\n' );
  fprintf ( 1, '  CHAIN_NUM = %d\n', chain_num );
  fprintf ( 1, '  CR_NUM    = %d\n', cr_num );
  fprintf ( 1, '  GEN_NUM   = %d\n', gen_num );
  fprintf ( 1, '  PAIR_NUM  = %d\n', pair_num );
  fprintf ( 1, '  PAR_NUM   = %d\n', par_num );
%
%  Initialize RNGLIB, if used.
%  Set the current generator index to #1.
%
  initialize ( );
  g = 1;
  cgn_set ( g );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Current RNGLIB generator index = %d\n', g );
%
%  Get the problem parameter values.
%
  chain_filename = '';
  gr_filename = '';
  limits = zeros ( 2, par_num );
  restart_read_filename = '';
  restart_write_filename = '';

  [ chain_filename, gr_filename, gr_threshold, jumpstep, limits, ...
    printstep, restart_read_filename, restart_write_filename ] ...
    = problem_value ( par_num );
%
%  Print the problem sizes and parameters.
%
  input_print ( chain_filename, chain_num, cr_num, gr_filename, ...
    gr_threshold, jumpstep, limits, gen_num, pair_num, par_num, printstep, ...
    restart_read_filename, restart_write_filename );
%
%  Allocate and zero out memory.
%
  fit = zeros(chain_num,gen_num);
  z = zeros(par_num,chain_num,gen_num);
%
%  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_num = floor ( gen_num / printstep );

  [ gr, gr_conv, gr_count ] = gr_init ( gr_num, par_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'GR_PRINT\n' );
  fprintf ( 1, '  GR_CONV  = %d\n', gr_conv );
  fprintf ( 1, '  GR_COUNT = %d\n', gr_count );
  fprintf ( 1, '  GR_NUM   = %d\n', gr_num );
%
%  Set the first generation of the chains from restart data, or by sampling.
%
  if ( 0 < length ( restart_read_filename ) )
    [ fit, z ] = restart_read ( chain_num, gen_num, par_num, ...
      restart_read_filename );
  else
    [ fit, z ] = chain_init ( chain_num, gen_num, par_num, ...
      problem_prior_sample, problem_sample_likelihood );
  end

  chain_init_print ( chain_num, fit, gen_num, par_num, ...
    restart_read_filename, z );
%
%  Carry out the DREAM algorithm.
%
  [ gr_conv, gr_count, z ] = dream_algm ( chain_num, cr_num, fit, gen_num, ...
    gr, gr_conv, gr_count, gr_num, gr_threshold, jumprate_table, jumpstep, ...
    limits, pair_num, par_num, printstep, problem_prior_density, ...
    problem_sample_likelihood, z );
%
%  Save Gelman-Rubin statistics to a file.
%
  if ( 0 < length ( gr_filename ) )
    gr_write ( gr, gr_filename, gr_num, par_num, printstep );
  end
%
%  Save parameter values for all chains at last generation.
%
  if ( 0 < length ( restart_write_filename ) )
    restart_write ( chain_num, fit, gen_num, par_num, ...
      restart_write_filename, z );
  end
%
%  Write each chain to a separate file.
%
  if ( 0 < length ( chain_filename ) )
    chain_write ( chain_filename, chain_num, fit, gen_num, par_num, z );
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'dream():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
 
  return
end
