diff --git a/docs/README.jmd b/docs/README.jmd
index 62107927..44c0cd0f 100644
--- a/docs/README.jmd
+++ b/docs/README.jmd
@@ -135,6 +135,10 @@ vline!([1/2 - π/tmax, 1/2, 1/2 + π/tmax], c=:black, ls=:dash, subplot=2)
 @doc Proposals.RSlice
 ```
 ---
+```julia; echo=false
+@doc Proposals.HSlice
+```
+---
 ### Convergence
 ```julia; echo=false
 @doc dlogz_convergence
diff --git a/src/NestedSamplers.jl b/src/NestedSamplers.jl
index 7f5a8bdc..a855230b 100644
--- a/src/NestedSamplers.jl
+++ b/src/NestedSamplers.jl
@@ -18,7 +18,7 @@ import AbstractMCMC: AbstractSampler,
                         sample_end!,
                         bundle_samples,
                         mcmcsample
-using Distributions: quantile, UnivariateDistribution
+using Distributions: quantile, UnivariateDistribution, Uniform
 using MCMCChains: Chains
 import StatsBase
 using StatsFuns: logaddexp,
diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl
index 7be93611..92998e60 100644
--- a/src/proposals/Proposals.jl
+++ b/src/proposals/Proposals.jl
@@ -9,6 +9,7 @@ The available implementations are
 * [`Proposals.RStagger`](@ref) - random staggering away to a new point given an existing one
 * [`Proposals.Slice`](@ref) - slicing away to a new point given an existing one
 * [`Proposals.RSlice`](@ref) - random slicing away to a new point given an existing one
+* [`Proposals.HSlice`](@ref) - Hamiltonian slicing away to a new point using a series of random trajectories given an existing point
 """
 module Proposals
 
@@ -17,6 +18,8 @@ using ..Bounds
 using Random
 using LinearAlgebra
 using Parameters
+using Distributions
+using StatsBase
 
 export AbstractProposal
 
@@ -27,7 +30,7 @@ The abstract type for live point proposal algorithms.
 
 # Interface
 
-Each `AbstractProposal` must have this function, 
+Each `AbstractProposal` must have this function,
 ```julia
 (::AbstractProposal)(::AbstractRNG, point, loglstar, bounds, loglikelihood, prior_transform)
 ```
@@ -48,12 +51,12 @@ Propose a new live point by uniformly sampling within the bounding volume.
 struct Uniform <: AbstractProposal end
 
 function (::Uniform)(rng::AbstractRNG,
-    point::AbstractVector,
-    logl_star,
-    bounds::AbstractBoundingSpace,
-    loglike,
-    prior_transform)
-    
+                     point::AbstractVector,
+                     logl_star,
+                     bounds::AbstractBoundingSpace,
+                     loglike,
+                     prior_transform)
+
     ncall = 0
     while true
         u = rand(rng, bounds)
@@ -110,7 +113,7 @@ function (prop::RWalk)(rng::AbstractRNG,
             u_prop = @. point + prop.scale * du
             # inside unit-cube
             unitcheck(u_prop) && break
-            
+
             fail += 1
             nfail += 1
             # check if stuck generating bad numbers
@@ -133,7 +136,7 @@ function (prop::RWalk)(rng::AbstractRNG,
         end
         nc += 1
         ncall += 1
-        
+
         # check if stuck generating bad points
         if nc > 50 * prop.walks
             @warn "Random walk proposals appear to be extremely inefficient. Adjusting the scale-factor accordingly"
@@ -141,7 +144,7 @@ function (prop::RWalk)(rng::AbstractRNG,
             nc = accept = reject = 0
         end
     end
-    
+
     # update proposal scale using acceptance ratio
     update_scale!(prop, accept, reject, n)
 
@@ -160,7 +163,7 @@ end
 """
     Proposals.RStagger(;ratio=0.5, walks=25, scale=1)
 
-Propose a new live point by random staggering away from an existing live point. 
+Propose a new live point by random staggering away from an existing live point.
 This differs from the random walk proposal in that the step size here is exponentially adjusted
 to reach a target acceptance rate _during_ each proposal, in addition to _between_
 proposals.
@@ -193,7 +196,7 @@ function (prop::RStagger)(rng::AbstractRNG,
     accept = reject = fail = nfail = nc = ncall = 0
     stagger = 1
     local du, u_prop, logl_prop, u, v, logl
-    
+
     while nc < prop.walks || iszero(accept)
         # get proposed point
         while true
@@ -204,7 +207,7 @@ function (prop::RStagger)(rng::AbstractRNG,
             u_prop = @. point + prop.scale * stagger * du
             # inside unit-cube
             unitcheck(u_prop) && break
-            
+
             fail += 1
             nfail += 1
             # check if stuck generating bad numbers
@@ -213,7 +216,7 @@ function (prop::RStagger)(rng::AbstractRNG,
                 fail = 0
                 prop.scale *= exp(-1/n)
             end
-        end 
+        end
         # check proposed point
         v_prop = prior_transform(u_prop)
         logl_prop = loglike(v_prop)
@@ -227,7 +230,7 @@ function (prop::RStagger)(rng::AbstractRNG,
         end
         nc += 1
         ncall += 1
-        
+
         # adjust _stagger_ to target an acceptance ratio of `prop.ratio`
         ratio = accept / (accept + reject)
         if ratio > prop.ratio
@@ -235,7 +238,7 @@ function (prop::RStagger)(rng::AbstractRNG,
         elseif ratio < prop.ratio
             stagger /= exp(1 / reject)
         end
-        
+
         # check if stuck generating bad points
         if nc > 50 * prop.walks
             @warn "Random walk proposals appear to be extremely inefficient. Adjusting the scale-factor accordingly"
@@ -243,15 +246,15 @@ function (prop::RStagger)(rng::AbstractRNG,
             nc = accept = reject = 0
         end
     end
-            
+
     # update proposal scale using acceptance ratio
     update_scale!(prop, accept, reject, n)
-            
+
     return u, v, logl, ncall
 end
-      
+
 """
-    Proposals.Slice(;slices=5, scale=1)
+    Proposals.Slice(;slices=5, scale=1.0)
 
 Propose a new live point by a series of random slices away from an existing live point.
 This is a standard _Gibbs-like_ implementation where a single multivariate slice is a combination of `slices` univariate slices through each axis.
@@ -263,7 +266,7 @@ This is a standard _Gibbs-like_ implementation where a single multivariate slice
 @with_kw mutable struct Slice <: AbstractProposal
     slices = 5
     scale = 1.0
-    
+
     @assert slices ≥ 1 "Number of slices must be greater than or equal to 1"
     @assert scale ≥ 0 "Proposal scale must be non-negative"
 end
@@ -278,38 +281,40 @@ function (prop::Slice)(rng::AbstractRNG,
     # setup
     n = length(point)
     nc = nexpand = ncontract = 0
-    local u, v, logl 
-    
+    local u, v, logl
+
     # modifying axes and computing lengths
     axes = Bounds.axes(bounds)
     axes = prop.scale .* axes'
     # slice sampling loop
     for it in 1:prop.slices
-        
+
         # shuffle axis update order
         idxs = shuffle!(rng, collect(Base.axes(axes, 1)))
-        
+
         # slice sample along a random direction
         for idx in idxs
-            
+
             # select axis
             axis = axes[idx, :]
-            
-            u, v, logl, nc, nexpand, ncontract = sample_slice(rng, axis, point, logl_star, loglike, 
+
+            u, v, logl, nc, nexpand, ncontract = sample_slice(rng, axis, point, logl_star, loglike,
                                                               prior_transform, nc, nexpand, ncontract)
-        end # end of slice sample along a random direction             
-    end # end of slice sampling loop    
-    
+        end # end of slice sample along a random direction
+    end # end of slice sampling loop
+
     # update slice proposal scale based on the relative size of the slices compared to the initial guess
     prop.scale = prop.scale * nexpand / (2.0 * ncontract)
-         
+
     return u, v, logl, nc
-end   # end of function Slice  
+end   # end of function Slice
 
 """
-    Proposals.RSlice(;slices=5, scale=1)
+    Proposals.RSlice(;slices=5, scale=1.0)
+
 Propose a new live point by a series of random slices away from an existing live point.
 This is a standard _random_ implementation where each slice is along a random direction based on the provided axes.
+
 ## Parameters
 - `slices` is the minimum number of slices
 - `scale` is the proposal distribution scale, which will update _between_ proposals.
@@ -317,43 +322,43 @@ This is a standard _random_ implementation where each slice is along a random di
 @with_kw mutable struct RSlice <: AbstractProposal
     slices = 5
     scale = 1.0
-    
+
     @assert slices ≥ 1 "Number of slices must be greater than or equal to 1"
     @assert scale ≥ 0 "Proposal scale must be non-negative"
 end
 
 function (prop::RSlice)(rng::AbstractRNG,
-                       point::AbstractVector,
-                       logl_star,
-                       bounds::AbstractBoundingSpace,
-                       loglike,
-                       prior_transform;
-                       kwargs...)
+                        point::AbstractVector,
+                        logl_star,
+                        bounds::AbstractBoundingSpace,
+                        loglike,
+                        prior_transform;
+                        kwargs...)
     # setup
     n = length(point)
     nc = nexpand = ncontract = 0
     local u, v, logl
     # random slice sampling loop
     for it in 1:prop.slices
-    
-        # propose a direction on the unit n-sphere 
+
+        # propose a direction on the unit n-sphere
         drhat = randn(rng, n)
         drhat /= norm(drhat)
-        
+
         # transform and scale into parameter space
         axis = prop.scale .* (Bounds.axes(bounds) * drhat)
         u, v, logl, nc, nexpand, ncontract = sample_slice(rng, axis, point, logl_star, loglike,
                                                           prior_transform, nc, nexpand, ncontract)
     end # end of random slice sampling loop
-    
+
     # update random slice proposal scale based on the relative size of the slices compared to the initial guess
     prop.scale = prop.scale * nexpand / (2.0 * ncontract)
-    
+
     return u, v, logl, nc
 end # end of function RSlice
 
 # Method for slice sampling
-function sample_slice(rng, axis, u, logl_star, loglike, prior_transform, nc, nexpand, ncontract) 
+function sample_slice(rng, axis, u, logl_star, loglike, prior_transform, nc, nexpand, ncontract)
     # define starting window
     r = rand(rng)  # initial scale/offset
     u_l = @. u - r * axis  # left bound
@@ -364,7 +369,7 @@ function sample_slice(rng, axis, u, logl_star, loglike, prior_transform, nc, nex
         logl_l = -Inf
     end
     nc += 1
-    nexpand += 1 
+    nexpand += 1
 
     u_r = u_l .+ axis # right bound
     if unitcheck(u_r)
@@ -372,38 +377,38 @@ function sample_slice(rng, axis, u, logl_star, loglike, prior_transform, nc, nex
         logl_r = loglike(v_r)
     else
         logl_r = -Inf
-    end    
+    end
     nc += 1
-    nexpand += 1 
+    nexpand += 1
 
     # stepping out left and right bounds
     while logl_l ≥ logl_star
         u_l .-= axis
-        if unitcheck(u_l)   
+        if unitcheck(u_l)
             v_l = prior_transform(u_l)
             logl_l = loglike(v_l)
         else
             logl_l = -Inf
         end
         nc += 1
-        nexpand += 1 
+        nexpand += 1
     end
 
     while logl_r ≥ logl_star
         u_r .+= axis
-        if unitcheck(u_r)   
+        if unitcheck(u_r)
             v_r = prior_transform(u_r)
             logl_r = loglike(v_r)
         else
             logl_r = -Inf
         end
         nc += 1
-        nexpand += 1 
+        nexpand += 1
     end
 
     # sample within limits. If the sample is not valid, shrink the limits until the `logl_star` bound is hit
     window_init = norm(u_r - u_l)  # initial window size
-    while true 
+    while true
         # define slice and window
         u_hat = u_r - u_l
         window = norm(u_hat)
@@ -414,7 +419,7 @@ function sample_slice(rng, axis, u, logl_star, loglike, prior_transform, nc, nex
         # propose a new position
         r = rand(rng)
         u_prop = @. u_l + r * u_hat
-        if unitcheck(u_prop) 
+        if unitcheck(u_prop)
             v_prop = prior_transform(u_prop)
             logl_prop = loglike(v_prop)
         else
@@ -437,7 +442,505 @@ function sample_slice(rng, axis, u, logl_star, loglike, prior_transform, nc, nex
                 error("Slice sampler has failed to find a valid point.")
             end
         end
-    end # end of sample within limits while    
-end    
-    
+    end # end of sample within limits while
+end
+
+"""
+    Proposals.HSlice(;slices=5, scale=1.0, grad = nothing, max_move = nothing, fmove = 0.9 compute_jac = false)
+
+Propose a new live point by "Hamiltonian" Slice Sampling using a series of random trajectories away from an existing live point.
+Each trajectory is based on the provided axes and samples are determined by moving forwards/ backwards in time until the trajectory hits an edge
+and approximately reflecting off the boundaries.
+After a series of reflections is established, a new live point is proposed by slice sampling across the entire path.
+
+## Parameters
+- `slices` is the minimum number of slices
+- `scale` is the proposal distribution scale, which will update _between_ proposals
+- `grad` is the gradient of the log-likelihood with respect to the unit cube
+- `max_move` is the limit for `ncall`, which is the maximum number of timesteps allowed per proposal forwards and backwards in time
+- `fmove` is the target fraction of samples that are proposed along a trajectory (i.e. not reflecting)
+- `compute_jac` a true/false statement for whether to compute and apply the Jacobian `dv/du` from the target space `v` to the unit cube `u` when evaluating the `grad`.
+"""
+@with_kw mutable struct HSlice <: AbstractProposal
+    slices = 5
+    scale = 1.0
+    grad = nothing
+    max_move = 100
+    fmove = 0.9
+    compute_jac = false
+
+    @assert slices ≥ 1 "Number of slices must be greater than or equal to 1"
+    @assert scale ≥ 0 "Proposal scale must be non-negative"
+    @assert max_move ≥ 1 "The limit for ncall must be greater than or equal to 1"
+end
+
+function (prop::HSlice)(rng::AbstractRNG,
+                        point::AbstractVector,
+                        logl_star,
+                        bounds::AbstractBoundingSpace,
+                        loglike,
+                        prior_transform;
+                        kwargs...)
+
+    # setup
+    n = length(point)
+    jitter = 0.25 # 25% jitter
+    nc = nmove = nreflect = ncontract = 0
+    local nodes_l, nodes_m, nodes_r, u_l, u_r, v_l, v_r, u_out, u_in, vel, jac, reverse, reflect, ncall
+
+    # Hamiltonian slice sampling loop
+    for it in 1:prop.slices
+        # define the left, inner, and right nodes for a given chord
+        # slice sampling will be done using these chords
+        nodes_l = []
+        nodes_m = []
+        nodes_r = []
+
+        # propose a direction on the unit n-sphere
+        drhat = randn(rng, n)
+        drhat /= norm(drhat)
+
+        # transform and scale based on past tuning
+        axes = Bounds.axes(bounds)
+        axis = (prop.scale * 0.01) .* (axes * drhat)
+
+        # creating starting window
+        vel = axis    # current velocity
+        u_l = @. point - rand(rng, Uniform(1.0 - jitter, 1.0 + jitter), n) * vel
+        u_r = @. point + rand(rng, Uniform(1.0 - jitter, 1.0 + jitter), n) * vel
+        append!(nodes_l, u_l)
+        append!(nodes_m, point)
+        append!(nodes_r, u_r)
+
+        # progress right (i.e. forwards in time)
+        reverse = false
+        reflect = false
+        u_r = point
+        ncall = 0
+
+        while ncall <= prop.max_move
+
+            # iterate until the edge of the distribution is bracketed
+            append!(nodes_l, u_r)
+            u_out = nothing
+            u_in = []
+
+            while true
+
+                # step forward
+                u_r .+= rand(rng, Uniform(1.0 - jitter, 1.0 + jitter), n) * vel
+
+                # evaluate point
+                if unitcheck(u_r)
+                    v_r = prior_transform(u_r)
+                    logl_r = loglike(v_r)
+                    nc += 1
+                    ncall += 1
+                    nmove += 1
+                else
+                    logl_r = -Inf
+                end
+
+                # check if the log-likelihood constraint is satisfied
+                # (i.e. if in or out of bounds)
+
+                if logl_r < logl_star
+                    if reflect
+                        # if out of bounds and just reflected, then reverse the direction and terminate immediately
+                        reverse = true
+                        pop!(nodes_l)   # remove since chord does not exist
+                        break
+                    else
+                        # if already in bounds, then safe
+                        u_out = u_r
+                        logl_out = logl_r
+                    end
+                    # check if gradients can be computed assuming termination is with the current `u_out`
+                    if isfinite(logl_out)
+                        reverse = false
+                    else
+                        reverse = true
+                    end
+                else
+                    reflect = false
+                    append!(u_in, u_r)
+                end
+
+                # check if the edge is bracketed
+                if u_out != nothing
+                    break
+                end
+            end
+
+            # define the rest of chord
+            if length(nodes_l) == length(nodes_r) + 1
+                try
+                    u_in = u_in[rand(1:length(u_in))]    # pick point randomly ## or use u_in = u_in[rand(1:end)]?
+                catch   ## is it an exception?
+                    u_in = point
+                    ## no pass statement here?
+                end
+                append!(nodes_m, u_in)
+                append!(nodes_r, u_out)
+            end
+
+            # check if turned around
+            if reverse
+                break
+            end
+
+            # reflect off the boundary
+            u_r = u_out
+            logl_r = logl_out
+            if prop.grad == nothing
+                # if the gradient is not provided, approximate it numerically using 2nd-order methods
+                h = zeros(n)
+                for i in 1:n
+                    u_r_l = u_r
+                    u_r_r = u_r
+
+                    # right side
+                    u_r_r[i] += 1e-10
+                    if unitcheck(u_r_r)
+                        v_r_r = prior_transform(u_r_r)
+                        logl_r_r = loglike(v_r_r)
+                    else
+                        logl_r_r = -Inf
+                        reverse = true    # cannot compute gradient
+                    end
+                    nc += 1
+
+                    # left side
+                    u_r_l[i] += 1e-10
+                    if unitcheck(u_r_l)
+                        v_r_l = prior_transform(u_r_l)
+                        logl_r_l = loglike(v_r_l)
+                    else
+                        logl_r_l = -Inf
+                        reverse = true    # cannot compute gradient
+                    end
+
+                    if reverse
+                        break    # give up because have to turn around
+                    end
+                    nc += 1
+
+                    # compute dlnl/du
+                    h[i] = (logl_r_r - logl_r_l) / 2e-10
+                end
+            else
+                # if the gradient is provided, evaluate it
+                h = prop.grad(v_r)   ## check this step, include formula for grad
+
+                if prop.compute_jac
+                    jac = []
+
+                    # evaluate and apply Jacobian dv/du if gradient is defined as d(lnL)/dv instead of d(lnL)/du
+                    for i in 1:n
+                        u_r_l = u_r
+                        u_r_r = u_r
+
+                        # right side
+                        u_r_r[i] += 1e-10
+                        if unitcheck(u_r_r)
+                            v_r_r = prior_transform(u_r_r)
+                        else
+                            reverse = true    # cannot compute Jacobian
+                            v_r_r = v_r    # assume no movement
+                        end
+
+                        # left side
+                        u_r_l[i] -= 1e-10
+                        if unitcheck(u_r_l)
+                            v_r_l = prior_transform(u_r_l)
+                        else
+                            reverse = true    # cannot compute Jacobian
+                            v_r_r = v_r    # assume no movement
+                        end
+
+                        if reverse
+                            break    # give up because have to turn around
+                        end
+
+                        append!(jac, ((v_r_r - v_r_l) / 2e-10))
+                    end
+
+                    jac = jac  ## ?? Line 1010 in dy
+                    h = dot(jac, h)    # apply Jacobian
+                end
+                nc += 1
+            end
+
+            # compute specular reflection off boundary
+            vel_ref = @. vel - 2 * h * dot(vel, h) / norm(h)^2    
+            dotprod = dot(vel_ref, vel)
+            dotprod /= norm(vel_ref) * norm(vel)
+
+            # check angle of reflection
+            if dotprod < -0.99
+                # the reflection angle is sufficiently small that it might as well be a reflection
+                reverse = true
+                break
+            else
+                # if reflection angle is sufficiently large, proceed as normal to the new position
+                vel = vel_ref
+                u_out = nothing
+                reflect = true
+                nreflect += 1
+            end
+        end
+
+        # progress left (i.e. backwards in time)
+        reverse = false
+        reflect = false
+        vel = -axis    # current velocity
+        u_l = point
+        ncall = 0
+
+        while ncall <= prop.max_move
+
+            # iterate until the edge of the distribution is bracketed
+            # a doubling approach is used to try and locate the bounds faster
+            append!(nodes_r, u_l)
+            u_out = nothing
+            u_in = []
+
+            while true
+
+                # step forward
+                u_l .+= rand(rng, Uniform(1.0 - jitter, 1.0 + jitter), n) * vel
+
+                # evaluate point
+                if unitcheck(u_l)
+                    v_l = prior_transform(u_l)
+                    logl_l = loglike(v_l)
+                    nc += 1
+                    ncall += 1
+                    nmove += 1
+                else
+                    logl_l = -Inf
+                end
+
+                # check if the log-likelihood constraint are satisfied (i.e. in or out of bounds)
+                if logl_l < logl_star
+                    if reflect
+                        # if out of bounds and just reflected, then reverse direction and terminate immediately
+                        reverse = true
+                        pop!(nodes_r)    # remove since chord does not exist
+                        break
+                    else
+                        # if already in bounds, then safe
+                        u_out = u_l
+                        logl_out = logl_l
+                    end
+
+                    # check if gradients can be computed assuming there was termination with the current `u_out`
+                    if isfinite(logl_out)
+                        reverse = false
+                    else
+                        reverse = true
+                    end
+                else
+                    reflect = false
+                    append!(u_in, u_l)
+                end
+
+                # check if the edge is bracketed
+                if u_out != nothing
+                    break
+                end
+            end
+
+            # define the rest of chord
+            if length(nodes_r) == length(nodes_l) + 1
+                try
+                    u_in = u_in[rand(1:length(u_in))]  # pick point randomly  ## or use u_in = u_in[rand(1:end)] ?
+                catch   ## is it an exception?
+                    u_in = point
+                    ## no pass statement here?
+                end
+                append!(nodes_m, u_in)
+                append!(nodes_l, u_out)
+            end
+
+            # check if turned around
+            if reverse
+                break
+            end
+
+            # reflect off the boundary
+            u_l = u_out
+            logl_l = logl_out
+
+            if prop.grad == nothing
+
+                # if the gradient is not provided, attempt to approximate it numerically using 2nd-order methods
+                h = zeros(n)
+                for i in 1:n
+                    u_l_l = u_l
+                    u_l_r = u_l
+
+                    # right side
+                    u_l_r[i] += 1e-10
+                    if unitcheck(u_l_r)
+                        v_l_r = prior_transform(u_l_r)
+                        logl_l_r = loglike(v_l_r)
+                    else
+                        logl_l_r = -Inf
+                        reverse = true    # cannot compute gradient
+                    end
+                    nc += 1
+
+                    # left side
+                    u_l_l[i] -= 1e-10
+                    if unitcheck(u_l_l)
+                        v_l_l = prior_transform(u_l_l)
+                        logl_l_l = loglike(v_l_l)
+                    else
+                        logl_l_l = -Inf
+                        reverse = true    # cannot compute gradient
+                    end
+
+                    if reverse
+                        break    # give up because have to turn around
+                    end
+                    nc += 1
+
+                    # compute dlnl/du
+                    h[i] = (logl_l_r - logl_l_l) / 2e-10
+                end
+            else
+                # if gradient is provided, evaluate it
+                h = prop.grad(v_l)   ## ? write the formula for grad
+                if prop.compute_jac
+                    jac = []
+
+                    # evaluate and apply Jacobian dv/du if gradient is defined as d(lnL)/dv instead of d(lnL)/du
+                    for i in 1:n
+                        u_l_l = u_l
+                        u_l_r = u_l
+
+                        # right side
+                        u_l_r[i] += 1e-10
+                        if unitcheck(u_l_r)
+                            v_l_r = prior_transform(u_l_r)
+                        else
+                            reverse = true    # cannot compute Jacobian
+                            v_l_r = v_l    # assume no movement
+                        end
+
+                        # left side
+                        u_l_l[i] -= 1e-10
+                        if unitcheck(u_l_l)
+                            v_l_l = prior_transform(u_l_l)
+                        else
+                            reverse = true    # cannot compute Jacobian
+                            v_l_r = v_l    # assume no movement
+                        end
+
+                        if reverse
+                            break    # give up because have to turn around
+                        end
+
+                        append!(jac, ((v_l_r - v_l_l) / 2e-10))
+                    end
+                    jac = jac  ## ?? Line 1148 in dy
+                    h = dot(jac, h)    # apply Jacobian
+                end
+                nc += 1
+            end
+            # compute specular reflection off boundary
+            vel_ref = @. vel - 2 * h * dot(vel, h) / norm(h)^2
+            dotprod = dot(vel_ref, vel)
+            dotprod /= norm(vel_ref) * norm(vel)
+
+            # check angle of reflection
+            if dotprod < -0.99
+                # the reflection angle is sufficiently small that it might as well be a reflection
+                reverse = true
+                break
+            else
+                # if the reflection angle is sufficiently large, proceed as normal to the new position
+                vel = vel_ref
+                u_out = nothing
+                reflect = true
+                nreflect += 1
+            end
+        end
+
+        # initialize lengths of cords
+        if length(nodes_l) > 1
+
+            # remove initial fallback chord
+            popfirst!(nodes_l)
+            popfirst!(nodes_m)
+            popfirst!(nodes_r)
+        end
+
+        nodes_l, nodes_m, nodes_r = (nodes_l, nodes_m, nodes_r)
+        Nchords = length(nodes_l)
+        axlen = zeros(Float64, Nchords)
+
+        for (i, (nl, nm, nr)) in enumerate(zip(nodes_l, nodes_m, nodes_r))
+            axlen[i] = norm(nr - nl)
+        end
+
+        # slice sample from all chords simultaneously, this is equivalent to slice sampling in *time* along trajectory
+        axlen_init = axlen
+
+        while true
+
+            # safety check
+            if any(axlen < 1e-5 * axlen_init)
+                error("Hamiltonian slice sampling appears to be stuck!")
+            end
+
+            # select chord
+            axprob = axlen / sum(axlen)
+            idx = sample(1:Nchords, ProbabilityWeights(axprob))
+
+            # define chord
+            u_l = nodes_l[idx]
+            u_m = nodes_m[idx]
+            u_r = nodes_r[idx]
+            u_hat = u_r - u_l
+            rprop = rand(rng)
+            u_prop = @. u_l + rprop * u_hat    # scale from left
+            if unitcheck(u_prop)
+                v_prop = prior_transform(u_prop)
+                logl_prop = loglike(v_prop)
+            else
+                logl_prop = -Inf
+            end
+            nc += 1
+            ncontract += 1
+
+            # if succeed, move to the new position
+            if logl_prop >= logl_star
+                u = u_prop
+                break
+            # if fail, check if the new point is to the left/right of the point interior to the bounds (`u_m`) and update the bounds accordingly
+            else
+                s = dot(u_prop - u_m, u_hat)    # check sign (+/-)
+                if s < 0    # left
+                    nodes_l[idx] = u_prop
+                    axlen[idx] *= 1 - rprop
+                elseif s > 0    # right
+                    nodes_r[idx] = u_prop
+                    axlen[idx] *= rprop
+                else
+                    error("Slice sampler has failed to find a valid point.")
+                end
+            end
+        end
+    end
+
+    # update the Hamiltonian slice proposal scale based on the relative amount of time spent moving vs reflecting
+    ## ?? ncontract: why is it set to 0, in line 229 of nessam
+    fmove = (1.0 * nmove) / (nmove + nreflect + ncontract + 2)
+    prop.scale *= exp((fmove - prop.fmove) / (max(prop.fmove, 1.0 - prop.fmove)))
+
+    return u_prop, v_prop, logl_prop, nc
+end    # end of function HSlice
+
 end # module Proposals
diff --git a/src/staticsampler.jl b/src/staticsampler.jl
index ee78c734..74be3854 100644
--- a/src/staticsampler.jl
+++ b/src/staticsampler.jl
@@ -39,7 +39,7 @@ Static nested sampler with `nactive` active points and `ndims` parameters.
 `bounds` declares the Type of [`Bounds.AbstractBoundingSpace`](@ref) to use in the prior volume. The available bounds are described by [`Bounds`](@ref). `proposal` declares the algorithm used for proposing new points. The available proposals are described in [`Proposals`](@ref). If `proposal` is `:auto`, will choose the proposal based on `ndims`
 * `ndims < 10` - [`Proposals.Uniform`](@ref)
 * `10 ≤ ndims ≤ 20` - [`Proposals.RWalk`](@ref)
-* `ndims > 20` - [`Proposals.Slice`](@ref)
+* `ndims > 20` - [`Proposals.HSlice`](@ref) if a `grad` (gradient) is provided and [`Proposals.Slice`](@ref) otherwise.
 
 The original nested sampling algorithm is roughly equivalent to using `Bounds.Ellipsoid` with `Proposals.Uniform`. The MultiNest algorithm is roughly equivalent to `Bounds.MultiEllipsoid` with `Proposals.Uniform`. The PolyChord algorithm is roughly equivalent to using `Proposals.RSlice`.
 
@@ -50,6 +50,7 @@ The original nested sampling algorithm is roughly equivalent to using `Bounds.El
     * `Proposals.RWalk` and `Proposals.RStagger` - `0.15 * walks`
     * `Proposals.Slice` - `0.9 * ndims * slices`
     * `Proposals.RSlice` - `2 * slices`
+    * `Proposals.HSlice` - `25.0 * slices`
 * `min_ncall` - The minimum number of iterations before trying to fit the first bound
 * `min_eff` - The maximum efficiency before trying to fit the first bound
 """
@@ -71,7 +72,11 @@ function Nested(ndims,
         elseif 10 ≤ ndims ≤ 20
             Proposals.RWalk() 
         else
-            Proposals.Slice()
+            if grad == nothing
+                Proposals.Slice()
+            else
+                Proposals.HSlice()
+            end
         end
     end
 
@@ -107,6 +112,7 @@ default_update_interval(p::Proposals.RWalk, ndims) = 0.15 * p.walks
 default_update_interval(p::Proposals.RStagger, ndims) = 0.15 * p.walks
 default_update_interval(p::Proposals.Slice, ndims) = 0.9 * ndims * p.slices
 default_update_interval(p::Proposals.RSlice, ndims) = 2.0 * p.slices
+default_update_interval(p::Proposals.HSlice, ndims) = 25.0 * p.slices
 
 
 function Base.show(io::IO, n::Nested)
diff --git a/test/proposals/proposals.jl b/test/proposals/proposals.jl
index 0661fe8c..0c9dd899 100644
--- a/test/proposals/proposals.jl
+++ b/test/proposals/proposals.jl
@@ -3,7 +3,8 @@ const PROPOSALS = [
     Proposals.RWalk(),
     Proposals.RStagger(),
     Proposals.Slice(),
-    Proposals.RSlice()
+    Proposals.RSlice(),
+    Proposals.HSlice()
 ]
 
 const BOUNDS = [
@@ -82,4 +83,17 @@ end
 
     @test_throws AssertionError Proposals.RSlice(slices=-2)
     @test_throws AssertionError Proposals.RSlice(scale=-3)
-end
\ No newline at end of file
+end
+
+@testset "HSlice" begin
+    prop = Proposals.HSlice()
+    @test prop.slices == 5
+    @test prop.scale == 1
+    @test prop.grad == nothing
+    @test prop.max_move == 100
+    @test prop.fmove == 0.9
+    @test prop.compute_jac == false
+
+    @test_throws AssertionError Proposals.HSlice(slices=-2)
+    @test_throws AssertionError Proposals.HSlice(scale=-3)
+end
diff --git a/test/sampler.jl b/test/sampler.jl
index 381c2e95..547bebad 100644
--- a/test/sampler.jl
+++ b/test/sampler.jl
@@ -10,6 +10,8 @@ using NestedSamplers: default_update_interval
     @test default_update_interval(Proposals.Slice(slices=10), 25) == 225
     @test default_update_interval(Proposals.RSlice(), 30) == 10  
     @test default_update_interval(Proposals.RSlice(slices=10), 25) == 20
+    @test default_update_interval(Proposals.HSlice(), 30) == 125
+    @test default_update_interval(Proposals.HSlice(slices=10), 25) == 250
 end
 
 spl = Nested(3, 100)
@@ -37,3 +39,7 @@ spl = Nested(10, 1000)
 spl = Nested(30, 1500)
 @test spl.proposal isa Proposals.Slice
 @test spl.update_interval == 202500
+
+spl = Nested(30, 1500)
+@test spl.proposal isa Proposals.HSlice
+@test spl.update_interval == 202500
diff --git a/test/sampling.jl b/test/sampling.jl
index 93776aa1..2f5152d1 100644
--- a/test/sampling.jl
+++ b/test/sampling.jl
@@ -28,7 +28,7 @@ end
 # end
 
 @testset "Gaussian - $bound, $P" for bound in [Bounds.NoBounds, Bounds.Ellipsoid, Bounds.MultiEllipsoid],
-                        P in [Proposals.Uniform, Proposals.RWalk, Proposals.RStagger, Proposals.Slice, Proposals.RSlice]
+                        P in [Proposals.Uniform, Proposals.RWalk, Proposals.RStagger, Proposals.Slice, Proposals.RSlice, Proposals.HSlice]
     σ = 0.1
     μ1 = ones(2)
     μ2 = -ones(2)