function samples = slice_sample(N, burn, logdist, xx, widths, step_out, varargin)
%SLICE_SAMPLE simple axis-aligned implementation of slice sampling for vectors
%
%     samples = slice_sample(N, burn, logdist, xx, widths, step_out, varargin)
%
% Inputs:
%             N  1x1  Number of samples to gather
%          burn  1x1  after burning period of this length
%       logdist  @fn  function logprobstar = logdist(xx, varargin{:})
%            xx  Dx1  initial state (or array with D elements)
%        widths  Dx1  or 1x1, step sizes for slice sampling
%      step_out bool  set to true if widths may sometimes be far too small
%      varargin   -   any extra arguments are passed on to logdist
%
% Outputs:
%      samples  DxN   samples stored in columns (regardless of original shape)
%
% Iain Murray May 2004, tweaks June 2009, a diagnostic added Feb 2010
% See Pseudo-code in David MacKay's text book p375

% startup stuff
D = numel(xx);
samples = zeros(D, N);
if numel(widths) == 1
    widths = repmat(widths, D, 1);
end
log_Px = feval(logdist, xx, varargin{:});

% Main loop
for ii = 1:(N+burn)
    fprintf('Iteration %d                 \r', ii - burn);
    log_uprime = log(rand) + log_Px;

    % Sweep through axes (simplest thing)
    for dd = 1:D
        x_l = xx;
        x_r = xx;
        xprime = xx;

        % Create a horizontal interval (x_l, x_r) enclosing xx
        rr = rand;
        x_l(dd) = xx(dd) - rr*widths(dd);
        x_r(dd) = xx(dd) + (1-rr)*widths(dd);
        if step_out
            % Typo in early editions of book. Book said compare to u, but it should say u'
            while (feval(logdist, x_l, varargin{:}) > log_uprime)
                x_l(dd) = x_l(dd) - widths(dd);
            end
            while (feval(logdist, x_r, varargin{:}) > log_uprime)
                x_r(dd) = x_r(dd) + widths(dd);
            end
        end

        % Inner loop:
        % Propose xprimes and shrink interval until good one found
        zz = 0;
        while 1
            zz = zz + 1;
            fprintf('Iteration %d   Step %d       \r', ii - burn, zz);
            xprime(dd) = rand()*(x_r(dd) - x_l(dd)) + x_l(dd);
            log_Px = feval(logdist, xprime, varargin{:});
            if log_Px > log_uprime
                break % this is the only way to leave the while loop
            else
                % Shrink in
                if xprime(dd) > xx(dd)
                    x_r(dd) = xprime(dd);
                elseif xprime(dd) < xx(dd)
                    x_l(dd) = xprime(dd);
                else
                    error('BUG DETECTED: Shrunk to current position and still not acceptable.');
                end
            end
        end
        xx(dd) = xprime(dd);
    end

    % Record samples
    if ii > burn
        samples(:, ii - burn) = xx(:);
    end
end
fprintf('\n');