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