From bc005c8545718806cc846f1b6a4c7c0c301a62fb Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Thu, 9 Jul 2020 22:06:46 +0530 Subject: [PATCH 01/26] [WIP] Implement Proposals.HSlice --- src/proposals/Proposals.jl | 41 +++++++++++++++++++++++++++++++++++++- 1 file changed, 40 insertions(+), 1 deletion(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index 3b842895..d083e7e1 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -247,6 +247,45 @@ function (prop::RStagger)(rng::AbstractRNG, return u, v, logl, ncall end - + +""" + Proposals.HSlice(;slices=5, scale=1) +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 +- `max_move` is the limit for `ncall` +- `compute_jac` a true/false statement for whether the Jacobian is needed. +""" +@with_kw mutable struct HSlice <: AbstractProposal + slices = 5 + scale = 1.0 + grad = nothing + max_move = 100 + 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" +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 + +end + + end # module Proposals From 5114b2e73be66268a65dd0d63a5eefd5c48de301 Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Fri, 10 Jul 2020 21:49:11 +0530 Subject: [PATCH 02/26] update Proposals.HSlice --- src/proposals/Proposals.jl | 146 ++++++++++++++++++++++++++++++++++++- 1 file changed, 145 insertions(+), 1 deletion(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index d083e7e1..1bdba60c 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -15,6 +15,7 @@ using ..Bounds using Random: AbstractRNG using LinearAlgebra using Parameters +using Distributions: Uniform export AbstractProposal @@ -270,6 +271,7 @@ After a series of reflections is established, a new live point is proposed by sl @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, @@ -284,7 +286,149 @@ function (prop::HSlice)(rng::AbstractRNG, n = length(point) jitter = 0.25 # 25% jitter nc = nmove = nreflect = ncontract = 0 - + local # incomplete + + # 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.tran_axes(bounds) + axis = dot(axes, drhat) * prop.scale * 0.01 + + # creating starting window + vel = axis # current velocity + u_l = @. u - Uniform(1.0 - jitter, 1.0 + jitter) * vel + u_r = @. u + Uniform(1.0 - jitter, 1.0 + jitter) * vel + append!(nodes_l, u_l) + append!(nodes_m, u) + append!(nodes_r, u_r) + + # progress right (i.e. forwards in time) + reverse = false + reflect = false + u_r = u + ncall = 0 + + while ncall <= 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 .+= Uniform(1.0 - jitter, 1.0 + jitter) * 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 ## incomplete line 938 + break + end + end + + # define the rest of chord + if ## incomplete + + end + + # check if turned around + if reverse + break + end + + # reflect off the boundary + u_r = u_out + logl_r = logl_out + if ## incomplete + # 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 = ## incomplete + + if compute_jac + end + end + end end From c7d179587668a591dc6a25d8bf9b4352b092c766 Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Sat, 11 Jul 2020 04:23:58 +0530 Subject: [PATCH 03/26] update Proposals.HSlice --- src/proposals/Proposals.jl | 303 ++++++++++++++++++++++++++++++++++++- 1 file changed, 300 insertions(+), 3 deletions(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index 1bdba60c..80015860 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -328,7 +328,7 @@ function (prop::HSlice)(rng::AbstractRNG, while true # step forward - u_r .+= Uniform(1.0 - jitter, 1.0 + jitter) * vel + u_r += Uniform(1.0 - jitter, 1.0 + jitter) * vel # evaluate point if unitcheck(u_r) @@ -426,10 +426,307 @@ function (prop::HSlice)(rng::AbstractRNG, h = ## incomplete if 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 + 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 = u + ncall = 0 + + while ncall <= 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 += Uniform(1.0 - jitter, 1.0 + jitter) * 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 ## incomplete + break + end + end + + # define the rest of chord + if ## incomplete + + end + + # check if turned around + if reverse + break + end + + # reflect off the boundary + u_l = u_out + logl_l = logl_out + + if grad ## incomplete + + # 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 end - end + else + # if gradient is provided, evaluate it + h = grad(v_l) + if 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 + 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 ## incomplete + + # remove initial fallback chord + ## incomplete end -end + + ## incomplete + + # slice sample from all chords simultaneously, this is equivalent to slice sampling in *time* along trajectory + axlen_init = axlen + + while true + + # safety check + if ## incomplete + + end + + # select chord + axprob = ## incomplete + idx = ## incomplete + + # 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 + end + + # 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 + ## incomplete + end + ## check all the loops, & end statements + ## also check where the placement of the update statement + ## also check placememt of return statement + end + + # update the Hamiltonian slice proposal scale based on the relative amount of time spent moving vs reflecting + ## ncontract ... check this formula step + ## check all formulas here, also check where to write prop.xyz and where not to write + fmove = (1.0 * nmove) / (nmove + nreflect + ncontract + 2) + norm = ## incomplete + prop.scale *= ## incomplete + + return u_prop, v_prop, logl_prop, nc + end # end of function HSlice end # module Proposals From 2af40c19f5cbcdff4e99187985acf89a1259f4de Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Sun, 12 Jul 2020 04:18:34 +0530 Subject: [PATCH 04/26] update Proposals.HSlice --- src/proposals/Proposals.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index 80015860..8318d765 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -657,10 +657,12 @@ function (prop::HSlice)(rng::AbstractRNG, end # initialize lengths of cords - if ## incomplete + if length(nodes_l) > 1 # remove initial fallback chord - ## incomplete + popfirst!(nodes_l) + popfirst!(nodes_m) + popfirst!(nodes_r) end ## incomplete From 142c3f45b8623fd45458d31266755a16b75f50ef Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Tue, 21 Jul 2020 22:43:46 +0530 Subject: [PATCH 05/26] Update Proposals.jl --- src/proposals/Proposals.jl | 484 +------------------------------------ 1 file changed, 1 insertion(+), 483 deletions(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index 8318d765..17817e40 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -12,10 +12,9 @@ module Proposals using ..Bounds -using Random: AbstractRNG +using Random using LinearAlgebra using Parameters -using Distributions: Uniform export AbstractProposal @@ -249,486 +248,5 @@ function (prop::RStagger)(rng::AbstractRNG, return u, v, logl, ncall end -""" - Proposals.HSlice(;slices=5, scale=1) -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 -- `max_move` is the limit for `ncall` -- `compute_jac` a true/false statement for whether the Jacobian is needed. -""" -@with_kw mutable struct HSlice <: AbstractProposal - slices = 5 - scale = 1.0 - grad = nothing - max_move = 100 - 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 # incomplete - - # 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.tran_axes(bounds) - axis = dot(axes, drhat) * prop.scale * 0.01 - - # creating starting window - vel = axis # current velocity - u_l = @. u - Uniform(1.0 - jitter, 1.0 + jitter) * vel - u_r = @. u + Uniform(1.0 - jitter, 1.0 + jitter) * vel - append!(nodes_l, u_l) - append!(nodes_m, u) - append!(nodes_r, u_r) - - # progress right (i.e. forwards in time) - reverse = false - reflect = false - u_r = u - ncall = 0 - - while ncall <= 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 += Uniform(1.0 - jitter, 1.0 + jitter) * 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 ## incomplete line 938 - break - end - end - - # define the rest of chord - if ## incomplete - - end - - # check if turned around - if reverse - break - end - - # reflect off the boundary - u_r = u_out - logl_r = logl_out - if ## incomplete - # 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 = ## incomplete - - if 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 - 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 = u - ncall = 0 - - while ncall <= 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 += Uniform(1.0 - jitter, 1.0 + jitter) * 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 ## incomplete - break - end - end - - # define the rest of chord - if ## incomplete - - end - - # check if turned around - if reverse - break - end - - # reflect off the boundary - u_l = u_out - logl_l = logl_out - - if grad ## incomplete - - # 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 - end - else - # if gradient is provided, evaluate it - h = grad(v_l) - if 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 - 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 - - ## incomplete - - # slice sample from all chords simultaneously, this is equivalent to slice sampling in *time* along trajectory - axlen_init = axlen - - while true - - # safety check - if ## incomplete - - end - - # select chord - axprob = ## incomplete - idx = ## incomplete - - # 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 - end - - # 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 - ## incomplete - end - ## check all the loops, & end statements - ## also check where the placement of the update statement - ## also check placememt of return statement - end - - # update the Hamiltonian slice proposal scale based on the relative amount of time spent moving vs reflecting - ## ncontract ... check this formula step - ## check all formulas here, also check where to write prop.xyz and where not to write - fmove = (1.0 * nmove) / (nmove + nreflect + ncontract + 2) - norm = ## incomplete - prop.scale *= ## incomplete - - return u_prop, v_prop, logl_prop, nc - end # end of function HSlice - - end # module Proposals From 0bd8438a79c9ed13b5a82b4248ce8a349dcee965 Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Tue, 21 Jul 2020 23:24:56 +0530 Subject: [PATCH 06/26] Update Proposals.jl --- src/proposals/Proposals.jl | 510 ++++++++++++++++++++++++++++++++++++- 1 file changed, 496 insertions(+), 14 deletions(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index 8fc82378..13da9a62 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -15,6 +15,7 @@ using ..Bounds using Random using LinearAlgebra using Parameters +using Distributions: Uniform export AbstractProposal @@ -42,11 +43,11 @@ 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 @@ -241,7 +242,7 @@ function (prop::RStagger)(rng::AbstractRNG, 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. ## Parameters @@ -295,7 +296,7 @@ function (prop::Slice)(rng::AbstractRNG, 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 @@ -311,12 +312,12 @@ This is a standard _random_ implementation where each slice is along a random di 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 @@ -428,4 +429,485 @@ function sample_slice(rng, axis, u, logl_star, loglike, prior_transform, nc, nex end # end of sample within limits while end -end # module Proposals \ No newline at end of file +""" + Proposals.HSlice(;slices=5, scale=1.0, grad = nothing, max_move = nothing, 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 +- `max_move` is the limit for `ncall` +- `compute_jac` a true/false statement for whether the Jacobian is needed. +""" +@with_kw mutable struct HSlice <: AbstractProposal + slices = 5 + scale = 1.0 + grad = nothing + max_move = 100 + 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 # incomplete + + # 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.tran_axes(bounds) + axis = dot(axes, drhat) * prop.scale * 0.01 + + # creating starting window + vel = axis # current velocity + u_l = @. u - Uniform(1.0 - jitter, 1.0 + jitter) * vel + u_r = @. u + Uniform(1.0 - jitter, 1.0 + jitter) * vel + append!(nodes_l, u_l) + append!(nodes_m, u) + append!(nodes_r, u_r) + + # progress right (i.e. forwards in time) + reverse = false + reflect = false + u_r = u + ncall = 0 + + while ncall <= 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 += Uniform(1.0 - jitter, 1.0 + jitter) * 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 ## incomplete line 938 + break + end + end + + # define the rest of chord + if ## incomplete + + end + + # check if turned around + if reverse + break + end + + # reflect off the boundary + u_r = u_out + logl_r = logl_out + if ## incomplete + # 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 = ## incomplete + + if 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 + 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 = u + ncall = 0 + + while ncall <= 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 += Uniform(1.0 - jitter, 1.0 + jitter) * 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 ## incomplete + break + end + end + + # define the rest of chord + if ## incomplete + + end + + # check if turned around + if reverse + break + end + + # reflect off the boundary + u_l = u_out + logl_l = logl_out + + if grad ## incomplete + + # 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 + end + else + # if gradient is provided, evaluate it + h = grad(v_l) + if 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 + 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 + + ## incomplete + + # slice sample from all chords simultaneously, this is equivalent to slice sampling in *time* along trajectory + axlen_init = axlen + + while true + + # safety check + if ## incomplete + + end + + # select chord + axprob = ## incomplete + idx = ## incomplete + + # 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 + end + + # 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 + ## incomplete + end + ## check all the loops, & end statements + ## also check where the placement of the update statement + ## also check placememt of return statement + end + + # update the Hamiltonian slice proposal scale based on the relative amount of time spent moving vs reflecting + ## ncontract ... check this formula step + ## check all formulas here, also check where to write prop.xyz and where not to write + fmove = (1.0 * nmove) / (nmove + nreflect + ncontract + 2) + norm = ## incomplete + prop.scale *= ## incomplete + + return u_prop, v_prop, logl_prop, nc +end # end of function HSlice + +end # module Proposals From 95036eb0e56acc309e57281d1d648c00e67a302f Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Thu, 23 Jul 2020 18:50:12 +0530 Subject: [PATCH 07/26] Update Proposals.jl --- src/proposals/Proposals.jl | 375 ++++++++++++++++++++----------------- 1 file changed, 198 insertions(+), 177 deletions(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index 13da9a62..66d0993e 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -7,6 +7,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 @@ -23,7 +24,7 @@ export AbstractProposal NestedSamplers.AbstractProposal 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,7 +49,7 @@ function (::Uniform)(rng::AbstractRNG, bounds::AbstractBoundingSpace, loglike, prior_transform) - + ncall = 0 while true u = rand(rng, bounds) @@ -103,7 +104,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 @@ -126,7 +127,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" @@ -134,7 +135,7 @@ function (prop::RWalk)(rng::AbstractRNG, nc = accept = reject = 0 end end - + # update proposal scale using acceptance ratio update_scale!(prop, accept, reject, n) @@ -152,7 +153,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. @@ -184,7 +185,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 @@ -195,7 +196,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 @@ -204,7 +205,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) @@ -218,7 +219,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 @@ -226,7 +227,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" @@ -234,13 +235,13 @@ 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.0) Propose a new live point by a series of random slices away from an existing live point. @@ -252,7 +253,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 @@ -267,33 +268,33 @@ 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.0) @@ -306,7 +307,7 @@ 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 @@ -324,25 +325,25 @@ function (prop::RSlice)(rng::AbstractRNG, 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 @@ -353,7 +354,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) @@ -361,38 +362,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) @@ -403,7 +404,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 @@ -426,13 +427,13 @@ 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, 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 +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 @@ -445,10 +446,10 @@ After a series of reflections is established, a new live point is proposed by sl @with_kw mutable struct HSlice <: AbstractProposal slices = 5 scale = 1.0 - grad = nothing + grad = nothing max_move = 100 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" @@ -461,13 +462,13 @@ function (prop::HSlice)(rng::AbstractRNG, loglike, prior_transform; kwargs...) - + # setup n = length(point) jitter = 0.25 # 25% jitter nc = nmove = nreflect = ncontract = 0 local # incomplete - + # Hamiltonian slice sampling loop for it in 1:prop.slices # define the left, inner, and right nodes for a given chord @@ -475,41 +476,41 @@ function (prop::HSlice)(rng::AbstractRNG, 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.tran_axes(bounds) + axes = Bounds.axes(bounds) axis = dot(axes, drhat) * prop.scale * 0.01 - + # creating starting window vel = axis # current velocity - u_l = @. u - Uniform(1.0 - jitter, 1.0 + jitter) * vel - u_r = @. u + Uniform(1.0 - jitter, 1.0 + jitter) * vel + u_l = @. point - Uniform(1.0 - jitter, 1.0 + jitter) * vel + u_r = @. point + Uniform(1.0 - jitter, 1.0 + jitter) * vel append!(nodes_l, u_l) - append!(nodes_m, u) + append!(nodes_m, point) append!(nodes_r, u_r) - + # progress right (i.e. forwards in time) reverse = false reflect = false - u_r = u + u_r = point ncall = 0 - + while ncall <= 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 += Uniform(1.0 - jitter, 1.0 + jitter) * vel - + # evaluate point if unitcheck(u_r) v_r = prior_transform(u_r) @@ -519,14 +520,14 @@ function (prop::HSlice)(rng::AbstractRNG, nmove += 1 else logl_r = -Inf - end - + 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 + # 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 @@ -544,34 +545,41 @@ function (prop::HSlice)(rng::AbstractRNG, else reflect = false append!(u_in, u_r) - end - + end + # check if the edge is bracketed - if ## incomplete line 938 + if u_out != nothing break - end + end end - + # define the rest of chord - if ## incomplete - + if length(nodes_l) == length(nodes_r) + 1 + try + u_in = u_in[rand(1:length(u_in))] # pick point randomly + ## incomplete + u_in = point + ## incomplete + end + append!(nodes_m, u_in) + append!(nodes_r, u_out) end - + # check if turned around if reverse break - end - + end + # reflect off the boundary u_r = u_out logl_r = logl_out - if ## incomplete + 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) @@ -580,9 +588,9 @@ function (prop::HSlice)(rng::AbstractRNG, else logl_r_r = -Inf reverse = true # cannot compute gradient - end + end nc += 1 - + # left side u_r_l[i] += 1e-10 if unitcheck(u_r_l) @@ -591,28 +599,28 @@ function (prop::HSlice)(rng::AbstractRNG, else logl_r_l = -Inf reverse = true # cannot compute gradient - end - + end + if reverse break # give up because have to turn around - end + end nc += 1 - + # compute dlnl/du h[i] = (logl_r_r - logl_r_l) / 2e-10 end - else + else # if the gradient is provided, evaluate it - h = ## incomplete - + h = grad(v_r) ## check this step, include formula for grad + if 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) @@ -620,69 +628,69 @@ function (prop::HSlice)(rng::AbstractRNG, else reverse = true # cannot compute Jacobian v_r_r = v_r # assume no movement - end - + end + # left side u_r_l[i] -= 1e-10 if unitcheck(u_r_l) v_r_l = prior_transform(u_r_l) - else + else reverse = true # cannot compute Jacobian v_r_r = v_r # assume no movement - end - + 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 + end + + jac = jac ## ?? 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 + # 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 end - + # progress left (i.e. backwards in time) reverse = false reflect = false - vel = axis # current velocity - u_l = u + vel = -axis # current velocity + u_l = point ncall = 0 - + while ncall <= 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 += Uniform(1.0 - jitter, 1.0 + jitter) * vel - + # evaluate point if unitcheck(u_l) v_l = prior_transform(u_l) @@ -692,8 +700,8 @@ function (prop::HSlice)(rng::AbstractRNG, nmove += 1 else logl_l = -Inf - end - + end + # check if the log-likelihood constraint are satisfied (i.e. in or out of bounds) if logl_l < logl_star if reflect @@ -705,162 +713,175 @@ function (prop::HSlice)(rng::AbstractRNG, # if already in bounds, then safe u_out = u_l logl_out = logl_l - end - + 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 + end + else reflect = false append!(u_in, u_l) - end - + end + # check if the edge is bracketed - if u_out ## incomplete + if u_out != nothing break - end - end - + end + end + # define the rest of chord - if ## incomplete - - end - + if length(nodes_r) == length(nodes_l) + 1 + try + u_in = u_in[rand(1:length(u_in))] # pick point randomly + ## incomplete + u_in = point + ## incomplete + end + append!(nodes_m, u_in) + append!(nodes_l, u_out) + end + # check if turned around if reverse break - end - + end + # reflect off the boundary u_l = u_out logl_l = logl_out - - if grad ## incomplete - + + if 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 + 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 + else logl_l_l = -Inf reverse = true # cannot compute gradient - end - + 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 + h[i] = (logl_l_r - logl_l_l) / 2e-10 + end end else # if gradient is provided, evaluate it h = grad(v_l) if 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 + else reverse = true # cannot compute Jacobian v_l_r = v_l # assume no movement - end - + end + # left side u_l_l[i] -= 1e-10 if unitcheck(u_l_l) v_l_l = prior_transform(u_l_l) - else + else reverse = true # cannot compute Jacobian v_l_r = v_l # assume no movement - end - + 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 + jac = jac ## ?? h = dot(jac, h) # apply Jacobian - end + end nc += 1 - end - + 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 + # 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 end - + # initialize lengths of cords if length(nodes_l) > 1 - + # remove initial fallback chord popfirst!(nodes_l) popfirst!(nodes_m) popfirst!(nodes_r) end - - ## incomplete - + + ## incomplete LINE 1175 of dynesty + Nchords = length(nodes_l) + axlen = zeros(Float64, Nchords) + + for i ## incomplete line 1179 of dynesty + 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 ## incomplete - - end - + + end + # select chord axprob = ## incomplete idx = ## incomplete - + # define chord u_l = nodes_l[idx] u_m = nodes_m[idx] @@ -873,16 +894,16 @@ function (prop::HSlice)(rng::AbstractRNG, logl_prop = loglike(v_prop) else logl_prop = -Inf - end + end nc += 1 ncontract += 1 - + # if succeed, move to the new position if logl_prop >= logl_star u = u_prop break end - + # 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 (+/-) @@ -896,17 +917,17 @@ function (prop::HSlice)(rng::AbstractRNG, ## incomplete end ## check all the loops, & end statements - ## also check where the placement of the update statement + ## also check where the placement of the update statement ## also check placememt of return statement - end - + end + # update the Hamiltonian slice proposal scale based on the relative amount of time spent moving vs reflecting ## ncontract ... check this formula step ## check all formulas here, also check where to write prop.xyz and where not to write fmove = (1.0 * nmove) / (nmove + nreflect + ncontract + 2) norm = ## incomplete prop.scale *= ## incomplete - + return u_prop, v_prop, logl_prop, nc end # end of function HSlice From a82b7aca06997d36ccbcfb9ff9f3a94c9339b42d Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Thu, 23 Jul 2020 18:58:46 +0530 Subject: [PATCH 08/26] include Proposals.HSlice --- docs/README.jmd | 4 ++++ 1 file changed, 4 insertions(+) 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 From 99828f07b8deb0ce7e268b27516e97c3a75c7b4d Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Thu, 23 Jul 2020 19:05:42 +0530 Subject: [PATCH 09/26] include Uniform in Distributions --- src/NestedSamplers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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, From 3c39ae91079727169a9c43b889fe711aff71c68a Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Thu, 23 Jul 2020 19:23:27 +0530 Subject: [PATCH 10/26] Update Proposals.jl --- src/proposals/Proposals.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index 66d0993e..9115ac14 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -487,8 +487,8 @@ function (prop::HSlice)(rng::AbstractRNG, # creating starting window vel = axis # current velocity - u_l = @. point - Uniform(1.0 - jitter, 1.0 + jitter) * vel - u_r = @. point + Uniform(1.0 - jitter, 1.0 + jitter) * vel + u_l = @. point - rand(rng, Uniform(1.0 - jitter, 1.0 + jitter)) * vel + u_r = @. point + rand(rng, Uniform(1.0 - jitter, 1.0 + jitter)) * vel append!(nodes_l, u_l) append!(nodes_m, point) append!(nodes_r, u_r) @@ -509,7 +509,7 @@ function (prop::HSlice)(rng::AbstractRNG, while true # step forward - u_r += Uniform(1.0 - jitter, 1.0 + jitter) * vel + u_r += rand(rng, Uniform(1.0 - jitter, 1.0 + jitter)) * vel # evaluate point if unitcheck(u_r) @@ -689,7 +689,7 @@ function (prop::HSlice)(rng::AbstractRNG, while true # step forward - u_l += Uniform(1.0 - jitter, 1.0 + jitter) * vel + u_l += rand(rng, Uniform(1.0 - jitter, 1.0 + jitter)) * vel # evaluate point if unitcheck(u_l) From 3a0063337547af50632f83c3867891a2fca3c45d Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Thu, 23 Jul 2020 20:39:45 +0530 Subject: [PATCH 11/26] include HSlice --- src/staticsampler.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) 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) From 2b9ec69739e2bea328170ec1af84b0e9f6c809a3 Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Thu, 23 Jul 2020 20:48:28 +0530 Subject: [PATCH 12/26] include HSlice --- test/proposals/proposals.jl | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/test/proposals/proposals.jl b/test/proposals/proposals.jl index 0661fe8c..7311df4e 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,16 @@ 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.compute_jac == false + + @test_throws AssertionError Proposals.HSlice(slices=-2) + @test_throws AssertionError Proposals.HSlice(scale=-3) +end From eb48700f64087a1d9b100501b629dc8915eefd8f Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Thu, 23 Jul 2020 20:54:55 +0530 Subject: [PATCH 13/26] include HSlice --- test/sampler.jl | 6 ++++++ 1 file changed, 6 insertions(+) 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 From d40e37e5589d764f267db38afd32cc71a5040052 Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Thu, 23 Jul 2020 20:56:51 +0530 Subject: [PATCH 14/26] include Proposals.HSlice --- test/sampling.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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) From df89c30e39b938dffc1dfd3fcf3122a5aee058d6 Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Fri, 24 Jul 2020 03:45:19 +0530 Subject: [PATCH 15/26] update Proposals.HSlice --- src/proposals/Proposals.jl | 257 ++++++++++++++++++------------------- 1 file changed, 128 insertions(+), 129 deletions(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index 9115ac14..b23afa25 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -17,6 +17,7 @@ using Random using LinearAlgebra using Parameters using Distributions: Uniform +using StatsBase export AbstractProposal @@ -431,7 +432,7 @@ function sample_slice(rng, axis, u, logl_star, loglike, prior_transform, nc, nex end """ - Proposals.HSlice(;slices=5, scale=1.0, grad = nothing, max_move = nothing, compute_jac = false) + 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. @@ -439,15 +440,17 @@ After a series of reflections is established, a new live point is proposed by sl ## 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 -- `max_move` is the limit for `ncall` -- `compute_jac` a true/false statement for whether the Jacobian is needed. +- `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" @@ -467,7 +470,7 @@ function (prop::HSlice)(rng::AbstractRNG, n = length(point) jitter = 0.25 # 25% jitter nc = nmove = nreflect = ncontract = 0 - local # incomplete + local drhat, 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 @@ -499,7 +502,7 @@ function (prop::HSlice)(rng::AbstractRNG, u_r = point ncall = 0 - while ncall <= max_move + while ncall <= prop.max_move # iterate until the edge of the distribution is bracketed append!(nodes_l, u_r) @@ -556,10 +559,10 @@ function (prop::HSlice)(rng::AbstractRNG, # 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 - ## incomplete + u_in = u_in[rand(1:end)] # pick point randomly ## or use u_in = u_in[rand(1:length(u_in))] + catch ## is it an exception? u_in = point - ## incomplete + ## no pass statement here? end append!(nodes_m, u_in) append!(nodes_r, u_out) @@ -611,9 +614,9 @@ function (prop::HSlice)(rng::AbstractRNG, end else # if the gradient is provided, evaluate it - h = grad(v_r) ## check this step, include formula for grad + h = prop.grad(v_r) ## check this step, include formula for grad - if compute_jac + if prop.compute_jac jac = [] # evaluate and apply Jacobian dv/du if gradient is defined as d(lnL)/dv instead of d(lnL)/du @@ -646,7 +649,7 @@ function (prop::HSlice)(rng::AbstractRNG, append!(jac, ((v_r_r - v_r_l) / 2e-10)) end - jac = jac ## ?? + jac = jac ## ?? Line 1010 in dy h = dot(jac, h) # apply Jacobian end nc += 1 @@ -678,7 +681,7 @@ function (prop::HSlice)(rng::AbstractRNG, u_l = point ncall = 0 - while ncall <= max_move + 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 @@ -735,10 +738,10 @@ function (prop::HSlice)(rng::AbstractRNG, # 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 - ## incomplete + u_in = u_in[rand(1:end)] # pick point randomly ## or use u_in = u_in[rand(1:length(u_in))] + catch ## is it an exception? u_in = point - ## incomplete + ## no pass statement here? end append!(nodes_m, u_in) append!(nodes_l, u_out) @@ -753,7 +756,7 @@ function (prop::HSlice)(rng::AbstractRNG, u_l = u_out logl_l = logl_out - if grad == nothing + if prop.grad == nothing # if the gradient is not provided, attempt to approximate it numerically using 2nd-order methods h = zeros(n) @@ -790,143 +793,139 @@ function (prop::HSlice)(rng::AbstractRNG, # compute dlnl/du h[i] = (logl_l_r - logl_l_l) / 2e-10 end - end - else - # if gradient is provided, evaluate it - h = grad(v_l) - if compute_jac - jac = [] + 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 + # 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 + # 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 + # 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 + if reverse + break # give up because have to turn around + end - append!(jac, ((v_l_r - v_l_l) / 2e-10)) + append!(jac, ((v_l_r - v_l_l) / 2e-10)) + end + jac = jac ## ?? Line 1148 in dy + h = dot(jac, h) # apply Jacobian end - jac = jac ## ?? - h = dot(jac, h) # apply Jacobian + nc += 1 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) + # 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 + # 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 - end - - # initialize lengths of cords - if length(nodes_l) > 1 - - # remove initial fallback chord - popfirst!(nodes_l) - popfirst!(nodes_m) - popfirst!(nodes_r) - end - ## incomplete LINE 1175 of dynesty - Nchords = length(nodes_l) - axlen = zeros(Float64, Nchords) + # initialize lengths of cords + if length(nodes_l) > 1 - for i ## incomplete line 1179 of dynesty - axlen[i] = norm(nr - nl) - end + # remove initial fallback chord + popfirst!(nodes_l) + popfirst!(nodes_m) + popfirst!(nodes_r) + end - # slice sample from all chords simultaneously, this is equivalent to slice sampling in *time* along trajectory - axlen_init = axlen + nodes_l, nodes_m, nodes_r = (nodes_l, nodes_m, nodes_r) ## check this assignment + Nchords = length(nodes_l) + axlen = zeros(Float64, Nchords) - while true + for i, (nl, nm, nr) in enumerate(zip(nodes_l, nodes_m, nodes_r)) ## check if this is iterating correctly! + axlen[i] = norm(nr - nl) + end - # safety check - if ## incomplete + # slice sample from all chords simultaneously, this is equivalent to slice sampling in *time* along trajectory + axlen_init = axlen - end + while true - # select chord - axprob = ## incomplete - idx = ## incomplete + # safety check + if any(axlen < 1e-5 * axlen_init) + error("Hamiltonian slice sampling appears to be stuck!") + end - # 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 + # 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 - end + # if succeed, move to the new position + if logl_prop >= logl_star + u = u_prop + break + end - # 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 + # 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 - ## incomplete + 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 - ## check all the loops, & end statements - ## also check where the placement of the update statement - ## also check placememt of return statement + end end # update the Hamiltonian slice proposal scale based on the relative amount of time spent moving vs reflecting - ## ncontract ... check this formula step - ## check all formulas here, also check where to write prop.xyz and where not to write + ## ?? ncontract: why is it set to 0, in line 229 of nessam fmove = (1.0 * nmove) / (nmove + nreflect + ncontract + 2) - norm = ## incomplete - prop.scale *= ## incomplete + norm = max(prop.fmove, 1.0 - prop.fmove) + prop.scale *= exp((fmove - prop.fmove) / norm) return u_prop, v_prop, logl_prop, nc end # end of function HSlice From 4f145b85c05128c57d405a0da93b7a336c58efe1 Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Fri, 24 Jul 2020 03:52:58 +0530 Subject: [PATCH 16/26] test fmove in .HSlice --- test/proposals/proposals.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/proposals/proposals.jl b/test/proposals/proposals.jl index 7311df4e..0c9dd899 100644 --- a/test/proposals/proposals.jl +++ b/test/proposals/proposals.jl @@ -91,6 +91,7 @@ end @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) From f8dcf408eaf583f950aad8e8b6595eebac2a4382 Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Fri, 24 Jul 2020 20:52:36 +0530 Subject: [PATCH 17/26] edit pick point randomly in u_in --- src/proposals/Proposals.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index b23afa25..9025f4d5 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -559,7 +559,7 @@ function (prop::HSlice)(rng::AbstractRNG, # define the rest of chord if length(nodes_l) == length(nodes_r) + 1 try - u_in = u_in[rand(1:end)] # pick point randomly ## or use u_in = u_in[rand(1:length(u_in))] + 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? @@ -738,7 +738,7 @@ function (prop::HSlice)(rng::AbstractRNG, # define the rest of chord if length(nodes_r) == length(nodes_l) + 1 try - u_in = u_in[rand(1:end)] # pick point randomly ## or use u_in = u_in[rand(1:length(u_in))] + 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? From 022480bceb775adc923c945e54feaf5e90243b18 Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Wed, 29 Jul 2020 19:47:47 +0530 Subject: [PATCH 18/26] add parentheses (i, (nl, nm, nr)), end while --- src/proposals/Proposals.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index 9025f4d5..4e0c189b 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -861,11 +861,11 @@ function (prop::HSlice)(rng::AbstractRNG, popfirst!(nodes_r) end - nodes_l, nodes_m, nodes_r = (nodes_l, nodes_m, nodes_r) ## check this assignment + 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)) ## check if this is iterating correctly! + for (i, (nl, nm, nr)) in enumerate(zip(nodes_l, nodes_m, nodes_r)) axlen[i] = norm(nr - nl) end @@ -903,8 +903,6 @@ function (prop::HSlice)(rng::AbstractRNG, if logl_prop >= logl_star u = u_prop break - end - # 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 (+/-) From 07f984d24a7cd923ef9f86d4490fb79d2bd04509 Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Wed, 29 Jul 2020 20:32:33 +0530 Subject: [PATCH 19/26] remove Uniform from using Distributions --- src/proposals/Proposals.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index 4e0c189b..74a6661e 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -16,7 +16,7 @@ using ..Bounds using Random using LinearAlgebra using Parameters -using Distributions: Uniform +using Distributions using StatsBase export AbstractProposal From fceacc221a466956c05c9df93712ad1243bc7996 Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Wed, 29 Jul 2020 22:05:31 +0530 Subject: [PATCH 20/26] include using LinearAlgebra --- test/sampling.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/sampling.jl b/test/sampling.jl index 2f5152d1..20002d18 100644 --- a/test/sampling.jl +++ b/test/sampling.jl @@ -3,6 +3,7 @@ using AbstractMCMC using MCMCChains: Chains using StatsFuns using StatsBase +using LinearAlgebra @testset "Bundles" begin logl(x::AbstractVector) = exp(-x[1]^2 / 2) / √(2π) From 4923930a399b5156c290b0a25f9202b47975580f Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Wed, 29 Jul 2020 22:09:12 +0530 Subject: [PATCH 21/26] Update sampling.jl --- test/sampling.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/sampling.jl b/test/sampling.jl index 20002d18..2f5152d1 100644 --- a/test/sampling.jl +++ b/test/sampling.jl @@ -3,7 +3,6 @@ using AbstractMCMC using MCMCChains: Chains using StatsFuns using StatsBase -using LinearAlgebra @testset "Bundles" begin logl(x::AbstractVector) = exp(-x[1]^2 / 2) / √(2π) From 89cb7dc4bdc127e6aff068da61041a8b47a33cb7 Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Wed, 29 Jul 2020 22:17:39 +0530 Subject: [PATCH 22/26] Update Proposals.jl --- src/proposals/Proposals.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index 74a6661e..ad420d5b 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -922,8 +922,7 @@ function (prop::HSlice)(rng::AbstractRNG, # 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) - norm = max(prop.fmove, 1.0 - prop.fmove) - prop.scale *= exp((fmove - prop.fmove) / norm) + 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 From e0a1c1f5ce954f8f11cb1b1d59e142bd4cfd33a9 Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Wed, 29 Jul 2020 23:17:45 +0530 Subject: [PATCH 23/26] Update Proposals.jl --- src/proposals/Proposals.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index ad420d5b..9245d1cc 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -470,7 +470,7 @@ function (prop::HSlice)(rng::AbstractRNG, n = length(point) jitter = 0.25 # 25% jitter nc = nmove = nreflect = ncontract = 0 - local drhat, nodes_l, nodes_m, nodes_r, u_l, u_r, v_l, v_r, u_out, u_in, vel, jac, reverse, reflect, ncall + 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 From 97c0cdcd1604da4e95e4f97ee8829ab022d86dcc Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Wed, 29 Jul 2020 23:35:35 +0530 Subject: [PATCH 24/26] Update Proposals.jl --- src/proposals/Proposals.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index 9245d1cc..7eae0461 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -486,7 +486,7 @@ function (prop::HSlice)(rng::AbstractRNG, # transform and scale based on past tuning axes = Bounds.axes(bounds) - axis = dot(axes, drhat) * prop.scale * 0.01 + axis = (prop.scale * 0.01) .* (axes * drhat) # creating starting window vel = axis # current velocity From 520b563a764c54ec88a6046a63d695edd4592df0 Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Thu, 30 Jul 2020 00:15:09 +0530 Subject: [PATCH 25/26] Update Proposals.jl --- src/proposals/Proposals.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index 7eae0461..cd6f8aa9 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -490,8 +490,8 @@ function (prop::HSlice)(rng::AbstractRNG, # creating starting window vel = axis # current velocity - u_l = @. point - rand(rng, Uniform(1.0 - jitter, 1.0 + jitter)) * vel - u_r = @. point + rand(rng, Uniform(1.0 - jitter, 1.0 + jitter)) * vel + 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) @@ -512,7 +512,7 @@ function (prop::HSlice)(rng::AbstractRNG, while true # step forward - u_r += rand(rng, Uniform(1.0 - jitter, 1.0 + jitter)) * vel + u_r .+= rand(rng, Uniform(1.0 - jitter, 1.0 + jitter), n) * vel # evaluate point if unitcheck(u_r) @@ -656,7 +656,7 @@ function (prop::HSlice)(rng::AbstractRNG, end # compute specular reflection off boundary - vel_ref = vel - 2 * h * dot(vel, h) / norm(h)^2 + vel_ref = @. vel - 2 * h * dot(vel, h) / norm(h)^2 dotprod = dot(vel_ref, vel) dotprod /= norm(vel_ref) * norm(vel) @@ -692,7 +692,7 @@ function (prop::HSlice)(rng::AbstractRNG, while true # step forward - u_l += rand(rng, Uniform(1.0 - jitter, 1.0 + jitter)) * vel + u_l .+= rand(rng, Uniform(1.0 - jitter, 1.0 + jitter), n) * vel # evaluate point if unitcheck(u_l) @@ -834,7 +834,7 @@ function (prop::HSlice)(rng::AbstractRNG, nc += 1 end # compute specular reflection off boundary - vel_ref = vel - 2 * h * dot(vel, h) / norm(h)^2 + vel_ref = @. vel - 2 * h * dot(vel, h) / norm(h)^2 dotprod = dot(vel_ref, vel) dotprod /= norm(vel_ref) * norm(vel) From 54782977fbaab9f944f56cdc45af477bde2ec191 Mon Sep 17 00:00:00 2001 From: Saranjeet Kaur Date: Thu, 30 Jul 2020 01:03:43 +0530 Subject: [PATCH 26/26] whitespaces --- src/proposals/Proposals.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index cd6f8aa9..92998e60 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -1,6 +1,8 @@ """ NestedSamplers.Proposals + This module contains the different algorithms for proposing new points within a bounding volume in unit space. + The available implementations are * [`Proposals.Uniform`](@ref) - samples uniformly within the bounding volume * [`Proposals.RWalk`](@ref) - random walks to a new point given an existing one @@ -23,8 +25,11 @@ export AbstractProposal """ NestedSamplers.AbstractProposal + The abstract type for live point proposal algorithms. + # Interface + Each `AbstractProposal` must have this function, ```julia (::AbstractProposal)(::AbstractRNG, point, loglstar, bounds, loglikelihood, prior_transform) @@ -40,6 +45,7 @@ unitcheck(us) = all(u -> 0 < u < 1, us) """ Proposals.Uniform() + Propose a new live point by uniformly sampling within the bounding volume. """ struct Uniform <: AbstractProposal end @@ -66,7 +72,9 @@ Base.show(io::IO, p::Uniform) = print(io, "NestedSamplers.Proposals.Uniform") """ Proposals.RWalk(;ratio=0.5, walks=25, scale=1) + Propose a new live point by random walking away from an existing live point. + ## Parameters - `ratio` is the target acceptance ratio - `walks` is the minimum number of steps to take @@ -154,10 +162,12 @@ end """ Proposals.RStagger(;ratio=0.5, walks=25, scale=1) + 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. + ## Parameters - `ratio` is the target acceptance ratio - `walks` is the minimum number of steps to take @@ -245,8 +255,10 @@ end """ 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. + ## Parameters - `slices` is the minimum number of slices - `scale` is the proposal distribution scale, which will update _between_ proposals. @@ -299,8 +311,10 @@ end # end of function Slice """ 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. @@ -433,10 +447,12 @@ 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