function [ jump_dim, jump_num, jumprate ] = jumprate_choose ( cr, cr_index, ...
  cr_num, gen_index, jumprate_table, jumpstep, par_num )

%*****************************************************************************80
%
%% JUMPRATE_CHOOSE chooses a jump rate from the jump rate table.
%
%  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, real CR(CR_NUM), the CR values.
%
%    Input, integer CR_INDEX, the index of the CR.
%    1 <= CR_INDEX <= CR_NUM.
%
%    Input, integer CR_NUM, the total number of CR values.
%    1 <= CR_NUM.
%
%    Input, integer GEN_INDEX, the current generation.
%    1 <= GEN_INDEX <= GEN_NUM.
%
%    Input, real JUMPRATE_TABLE(PAR_NUM), the jump rate table.
%
%    Input, integer JUMPSTEP, forces a "long jump" every
%    JUMPSTEP generations.
%
%    Input, integer ( kind = 4 ) PAR_NUM, the total number of parameters.
%    1 <= PAR_NUM.
%
%  Output parameters:
%
%    Output, integer JUMP_DIM(PAR_NUM), the indexes of the
%    parameters to be updated.
%
%    Output, integer JUMP_NUM, the number of dimensions in which
%    a jump will be made.  0 <= JUMP_NUM <= PAR_NUM.
%
%    Output, real JUMPRATE, the jump rate.
%

%
%  Determine the dimensions that will be updated.
%
  jump_num = 0;
  jump_dim = zeros ( par_num, 1 );

  for i = 1 : par_num

    r = r8_uniform_01_sample ( );

    if ( 1.0 - cr(cr_index) < r )
      jump_num = jump_num + 1;
      jump_dim(jump_num) = i;
    end

  end
%
%  Calculate the general jump rate.
%
  if ( jump_num == 0 )
    jumprate = 0.0;
  else
    jumprate = jumprate_table(jump_num);
  end
%
%  If parameter dimension is 1, 2, or 3, fix the jump rate to 0.6.
%
  if ( par_num <= 3 )
    jumprate = 0.6;
  end
%
%  Determine if a long jump is forced.
%
  if ( mod ( gen_index - 1, jumpstep ) == 0 )
    jumprate = 0.98;
  end

  return
end
