function [ 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 )
    
%*****************************************************************************80
%
%% DREAM_ALGM gets a candidate parameter sample.
% 
%  Licensing:
%
%    This code is distributed under the GNU LGPL 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 parameters:
%
%    Input, integer CHAIN_NUM, the total number of chains.
%    3 <= CHAIN_NUM.
%
%    Input, integer CR_NUM, the total number of CR values.
%    1 <= CR_NUM.
%
%    Input, real FIT(CHAIN_NUM,GEN_NUM), the likelihood of
%    each sample.
%
%    Input, integer GEN_NUM, the total number of generations.
%    2 <= GEN_NUM.
%
%    Input, real GR(PAR_NUM,GR_NUM), 
%    the Gelman-Rubin R statistic.
%
%    Input, logical GR_CONV, the Gelman-Rubin convergence flag.
%
%    Input, integer GR_COUNT, counts the number of 
%    generations at which the Gelman-Rubin statistic has been computed.
%
%    Input, integer GR_NUM, the number of times the Gelman-Rubin
%    statistic may be computed.
%
%    Input, real GR_THRESHOLD, the convergence tolerance for
%    the Gelman-Rubin statistic.
%
%    Input, real JUMPRATE_TABLE(PAR_NUM), the jumprate table.
%
%    Input, integer JUMPSTEP, forces a "long jump" every
%    JUMPSTEP generations.
%
%    Input, real LIMITS(2,PAR_NUM), lower and upper bounds
%    for each parameter.
%
%    Input, integer PAIR_NUM, the number of pairs of 
%    crossover chains.
%    0 <= PAIR_NUM.
%
%    Input, integer PAR_NUM, the total number of parameters.
%    1 <= PAR_NUM.
%
%    Input, integer PRINTSTEP, the interval between generations on 
%    which the Gelman-Rubin statistic will be computed and written to a file.
%
%    Input, function value = problem_prior_density ( par_num, zp )
%
%    Input, function value = problem_sample_likelihood ( par_num, zp )
%
%    Input, real Z(PAR_NUM,CHAIN_NUM,GEN_NUM), the Markov chain 
%    sample data.
%
%  Output parameters:
%
%    Output, logical GR_CONV, the Gelman-Rubin convergence flag.
%
%    Output, integer GR_COUNT, counts the number of 
%    generations at which the Gelman-Rubin statistic has been computed.
%
%    Output, real Z(PAR_NUM,CHAIN_NUM,GEN_NUM), the Markov chain 
%    sample data.
%
%  Local parameters:
%
%    Local, integer CHAIN_INDEX, the index of the current chain.
%    1 <= CHAIN_INDEX <= CHAIN_NUM.
%
%    Local, real CR(CR_NUM), the CR values.
%
%    Local, real CR_DIS(CR_NUM), the CR distances.
%
%    Local, integer CR_INDEX, the index of the selected CR value.
%    1 <= CR_INDEX <= CR_NUM.
%
%    Local, real CR_PROB(CR_NUM), the probability of each CR.
%
%    Local, real CR_UPS(CR_NUM), the number of updates for each CR.
%
%    Local, integer GEN_INDEX, the index of the current generation.
%    1 <= GEN_INDEX <= GEN_NUM.
%
%    Local, real ZP(PAR_NUM), a candidate sample.
%
%    Local, integer ZP_ACCEPT, the number of candidates accepted.
%
%    Local, real ZP_ACCEPT_RATE, the rate at which generated
%    candidates were accepted.
%
%    Local, integer ZP_COUNT, the number of candidates generated.
%
%    Local, real ZP_RATIO, the Metropolis ratio for a candidate.
%
  zp_count = 0;
  zp_accept = 0;
%
%  Initialize the CR values.
%
  [ cr, cr_dis, cr_prob, cr_ups ] = cr_init ( cr_num );

  for gen_index = 2 : gen_num

    for chain_index = 1 : chain_num
%
%  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 = problem_sample_likelihood ( par_num, zp );

      zp_old(1:par_num) = z(1:par_num,chain_index,gen_index-1);
      zp_old_fit = fit(chain_index,gen_index-1);
%
%  Compute the Metropolis ratio for ZP versus ZP_OLD.
%
      pd1 = problem_prior_density ( par_num, zp );
      pd2 = problem_prior_density ( par_num, zp_old );

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

      zp_ratio = 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 )
        z(1:par_num,chain_index,gen_index) = zp(1:par_num);
        fit(chain_index,gen_index) = zp_fit;
        zp_accept = zp_accept + 1;
      else
        z(1:par_num,chain_index,gen_index) = zp_old(1:par_num);
        fit(chain_index,gen_index) = zp_old_fit;
      end
%
%  Update the CR distance.
%
      if ( ~ gr_conv )
        if ( 1 < cr_num )  
          [ cr_dis, cr_ups ] = cr_dis_update ( chain_index, chain_num, ...
            cr_dis, cr_index, cr_num, cr_ups, gen_index, gen_num, par_num, z );
        end
      end

    end
%
%  Update the multinomial distribution of CR.
%
    if ( ~ gr_conv )
      if ( 1 < cr_num )
        if ( mod ( gen_index, 10 ) == 0 )
          cr_prob = cr_prob_update ( cr_dis, cr_num, cr_ups );
        end
      end
    end
%
%  Every PRINTSTEP interval,
%  * compute the Gelman Rubin R statistic for this generation,
%    and determine if convergence has occurred.
%
    if ( mod ( gen_index, printstep ) == 0 )
      [ gr, gr_conv, gr_count ] = gr_compute ( chain_num, gen_index, ...
        gen_num, gr_count, gr_num, gr_threshold, par_num, z );
    end
%
%  Check for outlier chains.
%
    if ( ~ gr_conv )
      if ( mod ( gen_index, 10 ) == 0 )
        [ fit, z ] = chain_outliers ( chain_num, gen_index, gen_num, ...
          par_num, fit, z );
      end
    end

  end     
%
%  Compute the acceptance rate.
%
  zp_accept_rate = zp_accept / zp_count;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Candidates generated: %d\n', zp_count );
  fprintf ( 1, '  Candidates accepted:  %d\n', zp_accept );
  fprintf ( 1, '  The acceptance rate is %g\n', zp_accept_rate );

  return
end

