function [samples, energies, diagn] = hmc(f, x, options, gradf, varargin)%HMC	Hybrid Monte Carlo sampling.%%	Description%	SAMPLES = HMC(F, X, OPTIONS, GRADF) uses a  hybrid Monte Carlo%	algorithm to sample from the distribution P ~ EXP(-F), where F is the%	first argument to HMC. The Markov chain starts at the point X, and%	the function GRADF is the gradient of the `energy' function F.%%	HMC(F, X, OPTIONS, GRADF, P1, P2, ...) allows additional arguments to%	be passed to F() and GRADF().%%	[SAMPLES, ENERGIES, DIAGN] = HMC(F, X, OPTIONS, GRADF) also returns a%	log of the energy values (i.e. negative log probabilities) for the%	samples in ENERGIES and DIAGN, a structure containing diagnostic%	information (position, momentum and acceptance threshold) for each%	step of the chain in DIAGN.POS, DIAGN.MOM and DIAGN.ACC respectively.%	All candidate states (including rejected ones) are stored in%	DIAGN.POS.%%	[SAMPLES, ENERGIES, DIAGN] = HMC(F, X, OPTIONS, GRADF) also returns%	the ENERGIES (i.e. negative log probabilities) corresponding to the%	samples.  The DIAGN structure contains three fields:%%	POS the position vectors of the dynamic process.%%	MOM the momentum vectors of the dynamic process.%%	ACC the acceptance thresholds.%%	S = HMC('STATE') returns a state structure that contains the state of%	the two random number generators RAND and RANDN and the momentum of%	the dynamic process.  These are contained in fields  randstate,%	randnstate and mom respectively.  The momentum state is only used for%	a persistent momentum update.%%	HMC('STATE', S) resets the state to S.  If S is an integer, then it%	is passed to RAND and RANDN and the momentum variable is randomised.%	If S is a structure returned by HMC('STATE') then it resets the%	generator to exactly the same state.%%	The optional parameters in the OPTIONS vector have the following%	interpretations.%%	OPTIONS(1) is set to 1 to display the energy values and rejection%	threshold at each step of the Markov chain. If the value is 2, then%	the position vectors at each step are also displayed.%%	OPTIONS(5) is set to 1 if momentum persistence is used; default 0,%	for complete replacement of momentum variables.%%	OPTIONS(7) defines the trajectory length (i.e. the number of leap-%	frog steps at each iteration).  Minimum value 1.%%	OPTIONS(9) is set to 1 to check the user defined gradient function.%%	OPTIONS(14) is the number of samples retained from the Markov chain;%	default 100.%%	OPTIONS(15) is the number of samples omitted from the start of the%	chain; default 0.%%	OPTIONS(17) defines the momentum used when a persistent update of%	(leap-frog) momentum is used.  This is bounded to the interval [0,%	1).%%	OPTIONS(18) is the step size used in leap-frogs; default 1/trajectory%	length.%%	See also%	METROP%%	Copyright (c) Ian T Nabney (1996-2001)% Global variable to store state of momentum variables: set by set_state% Used to initialise variable if setglobal HMC_MOMif nargin <= 2  if ~strcmp(f, 'state')    error('Unknown argument to hmc');  end  switch nargin    case 1      samples = get_state(f);      return;    case 2      set_state(f, x);      return;  endenddisplay = options(1);if (round(options(5) == 1))  persistence = 1;  % Set alpha to lie in [0, 1)  alpha = max(0, options(17));  alpha = min(1, alpha);  salpha = sqrt(1-alpha*alpha);else  persistence = 0;endL = max(1, options(7)); % At least one step in leap-froggingif options(14) > 0  nsamples = options(14);else  nsamples = 100;	% Defaultendif options(15) >= 0  nomit = options(15);else  nomit = 0;endif options(18) > 0  step_size = options(18);	% Step size.else  step_size = 1/L;		% Default  endx = x(:)';		% Force x to be a row vectornparams = length(x);% Set up strings for evaluating potential function and its gradient.f = fcnchk(f, length(varargin));gradf = fcnchk(gradf, length(varargin));% Check the gradient evaluation.if (options(9))  % Check gradients  feval('gradchek', x, f, gradf, varargin{:});endsamples = zeros(nsamples, nparams);	% Matrix of returned samples.if nargout >= 2  en_save = 1;  energies = zeros(nsamples, 1);else  en_save = 0;endif nargout >= 3  diagnostics = 1;  diagn_pos = zeros(nsamples, nparams);  diagn_mom = zeros(nsamples, nparams);  diagn_acc = zeros(nsamples, 1);else  diagnostics = 0;endn = - nomit + 1;Eold = feval(f, x, varargin{:});	% Evaluate starting energy.nreject = 0;if (~persistence | isempty(HMC_MOM))  p = randn(1, nparams);		% Initialise momenta at randomelse  p = HMC_MOM;				% Initialise momenta from stored stateendlambda = 1;% Main loop.while n <= nsamples  xold = x;		    % Store starting position.  pold = p;		    % Store starting momenta  Hold = Eold + 0.5*(p*p'); % Recalculate Hamiltonian as momenta have changed  if ~persistence    % Choose a direction at random    if (rand < 0.5)      lambda = -1;    else      lambda = 1;    end  end  % Perturb step length.  epsilon = lambda*step_size*(1.0 + 0.1*randn(1));  % First half-step of leapfrog.  p = p - 0.5*epsilon*feval(gradf, x, varargin{:});  x = x + epsilon*p;    % Full leapfrog steps.  for m = 1 : L - 1    p = p - epsilon*feval(gradf, x, varargin{:});    x = x + epsilon*p;  end  % Final half-step of leapfrog.  p = p - 0.5*epsilon*feval(gradf, x, varargin{:});  % Now apply Metropolis algorithm.  Enew = feval(f, x, varargin{:});	% Evaluate new energy.  p = -p;				% Negate momentum  Hnew = Enew + 0.5*p*p';		% Evaluate new Hamiltonian.  a = exp(Hold - Hnew);			% Acceptance threshold.  if (diagnostics & n > 0)    diagn_pos(n,:) = x;    diagn_mom(n,:) = p;    diagn_acc(n,:) = a;  end  if (display > 1)    fprintf(1, 'New position is\n');    disp(x);  end  if a > rand(1)			% Accept the new state.    Eold = Enew;			% Update energy    if (display > 0)      fprintf(1, 'Finished step %4d  Threshold: %g\n', n, a);    end  else					% Reject the new state.    if n > 0       nreject = nreject + 1;    end    x = xold;				% Reset position     p = pold;   			% Reset momenta    if (display > 0)      fprintf(1, '  Sample rejected %4d.  Threshold: %g\n', n, a);    end  end  if n > 0    samples(n,:) = x;			% Store sample.    if en_save       energies(n) = Eold;		% Store energy.    end  end  % Set momenta for next iteration  if persistence    p = -p;    % Adjust momenta by a small random amount.    p = alpha.*p + salpha.*randn(1, nparams);  else    p = randn(1, nparams);	% Replace all momenta.  end  n = n + 1;endif (display > 0)  fprintf(1, '\nFraction of samples rejected:  %g\n', ...    nreject/(nsamples));endif diagnostics  diagn.pos = diagn_pos;  diagn.mom = diagn_mom;  diagn.acc = diagn_acc;end% Store final momentum value in global so that it can be retrieved laterHMC_MOM = p;return% Return complete state of sampler (including momentum)function state = get_state(f)global HMC_MOMstate.randstate = rand('state');state.randnstate = randn('state');state.mom = HMC_MOM;return% Set complete state of sampler (including momentum) or just set randn% and rand with integer argument.function set_state(f, x)global HMC_MOMif isnumeric(x)  rand('state', x);  randn('state', x);  HMC_MOM = [];else  if ~isstruct(x)    error('Second argument to hmc must be number or state structure');  end  if (~isfield(x, 'randstate') | ~isfield(x, 'randnstate') ...      | ~isfield(x, 'mom'))    error('Second argument to hmc must contain correct fields')  end  rand('state', x.randstate);  randn('state', x.randnstate);  HMC_MOM = x.mom;endreturn