Skip to content

Commit

Permalink
Move npoints, agepoints, and Tpoints inside MCMC function
Browse files Browse the repository at this point in the history
  • Loading branch information
@brenhinkeller
brenhinkeller committed Nov 5, 2024
1 parent 1f0dde9 commit 9a442f3
Show file tree
Hide file tree
Showing 5 changed files with 51 additions and 58 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "Thermochron"
uuid = "b9a8379e-ee1a-4388-b7ca-7852af1cef08"
authors = ["C. Brenhin Keller and collaborators"]
version = "0.9.4"
version = "0.10.0"

[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Expand Down
17 changes: 4 additions & 13 deletions examples/ZrnHeInversionVartCryst.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@
age = copy(data.HeAge)
age_sigma = copy(ds.HeAge_Ma_sigma_raw)
age_sigma_empirical = copy(data.HeAge_sigma)
min_rel_uncert = 10/100 # 10% minmum relative age uncertainty
min_rel_uncert = 5/100 # 5% minmum relative age uncertainty

if any(ta)
# Standard deviation of a Gaussian kernel in eU space, representing the
Expand All @@ -66,6 +66,7 @@
for i ∈ findall(ta)
W = normpdf.(eU[i], σeU, eU[ta])
σ_external = nanstd(age[ta], W) # Weighted standard deviation
σ_external /= sqrt(2) # Assume half of variance is unknown external uncertainty
σ_internal = age_sigma[i]
age_sigma_empirical[i] = sqrt(σ_external^2 + σ_internal^2)
age_sigma_empirical[i] = max(age_sigma_empirical[i], min_rel_uncert*age[i])
Expand Down Expand Up @@ -169,19 +170,9 @@

## --- Invert for maximum likelihood t-T path

# This is where the "transdimensional" part comes in
agepoints = zeros(Float64, model.maxpoints) # Array of fixed size to hold all optional age points
Tpoints = zeros(Float64, model.maxpoints) # Array of fixed size to hold all optional age points

# Fill some intermediate points to give the MCMC something to work with
Tr = 150 # Residence temperature
npoints = model.minpoints
agepoints[1:npoints] .= range(0, model.tinit, npoints)
Tpoints[1:npoints] .= Tr # Degrees C

# Run Markov Chain
@time (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, npoints, agepoints, Tpoints, boundary, constraint, detail)
# @time (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC_varkinetics(data, model, npoints, agepoints, Tpoints, boundary, constraint, detail)
@time (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, boundary, constraint, detail)
# @time (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC_varkinetics(data, model, boundary, constraint, detail)
@info """tpointdist & Tpointdist collected, size: $(size(Tpointdist))
Mean log-likelihood: $(nanmean(view(lldist, model.burnin:model.nsteps)))
Mean acceptance rate: $(nanmean(view(acceptancedist, model.burnin:model.nsteps)))
Expand Down
44 changes: 30 additions & 14 deletions src/inversion.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,10 +12,8 @@
tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist = MCMC(data, model, npoints, agepoints, Tpoints, constraint, boundary)
```
"""
function MCMC(data::NamedTuple, model::NamedTuple, npoints::Int, agepoints::DenseVector{T}, Tpoints::DenseVector{T}, boundary::Boundary{T}, constraint::Constraint{T}=Constraint(T), detail::DetailInterval{T}=DetailInterval(T)) where T <: AbstractFloat
# Sanitize inputs
@assert firstindex(agepoints) === 1
@assert firstindex(Tpoints) === 1
function MCMC(data::NamedTuple, model::NamedTuple, boundary::Boundary{T}, constraint::Constraint{T}=Constraint(T), detail::DetailInterval{T}=DetailInterval(T)) where T <: AbstractFloat
# Process inputs
halfwidth = T.(data.halfwidth)::DenseVector{T}
crystAge = T.(data.crystAge)::DenseVector{T}
HeAge = T.(data.HeAge)::DenseVector{T}
Expand All @@ -24,26 +22,36 @@
Th = T.(data.Th)::DenseVector{T}
Sm = (haskey(data, :Sm) ? T.(data.Th) : zeros(T, size(U)))::DenseVector{T}
nsteps = (haskey(model, :nsteps) ? model.nsteps : 10^6)::Int
maxpoints = (haskey(model, :maxpoints) ? model.maxpoints : 50)::Int
minpoints = (haskey(model, :minpoints) ? model.minpoints : 1)::Int
maxpoints = (haskey(model, :maxpoints) ? model.maxpoints : 50)::Int
npoints = (haskey(model, :npoints) ? model.npoints : minpoints)::Int
totalpoints = maxpoints + boundary.npoints + constraint.npoints::Int
simplified = (haskey(model, :simplified) ? model.simplified : false)::Bool
boundarytype = (haskey(model, :boundarytype) ? model.boundarytype : :hard)::Symbol
dynamicjumping = (haskey(model, :dynamicjumping) ? model.dynamicjumping : false)::Bool
dTmax = T(haskey(model, :dTmax) ? model.dTmax : 10.)::T
dTmax_sigma = T(haskey(model, :dTmax_sigma) ? model.dTmax_sigma : 5.)::T
dTmax = T(haskey(model, :dTmax) ? model.dTmax : 10)::T
dTmax_sigma = T(haskey(model, :dTmax_sigma) ? model.dTmax_sigma : 5)::T
agesteps = T.(model.agesteps)::DenseVector{T}
tsteps = T.(model.tsteps)::DenseVector{T}
tinit = T(model.tinit)::T
tnow = zero(T) # It's always time zero today
Tinit = T(model.Tinit)::T
Tnow = T(model.Tnow)::T
Tr = T(haskey(model, :Tr) ? model.Tr : (Tinit+Tnow)/2)::T
dt = T(model.dt)::T
dr = T(model.dr)::T
σmodel = T(model.σmodel)::T
σannealing = T(model.σannealing)::T
λannealing = T(model.λannealing)::T

# Arrays to hold all t and T points (up to npoints=maxpoints)
agepoints = zeros(T, maxpoints)
Tpoints = zeros(T, maxpoints)

# Fill some intermediate points to give the MCMC something to work with
agepoints[1:npoints] .= range(tnow, tinit, length=npoints)
Tpoints[1:npoints] .= Tr # Degrees C

# Calculate number of boundary and unconformity points and allocate buffer for interpolating
agepointbuffer = similar(agepoints, totalpoints)::DenseVector{T}
Tpointbuffer = similar(agepoints, totalpoints)::DenseVector{T}
Expand Down Expand Up @@ -269,10 +277,8 @@
end
export MCMC

function MCMC_varkinetics(data::NamedTuple, model::NamedTuple, npoints::Int, agepoints::DenseVector{T}, Tpoints::DenseVector{T}, boundary::Boundary{T}, constraint::Constraint{T}=Constraint(T), detail::DetailInterval{T}=DetailInterval(T)) where T <: AbstractFloat
# Sanitize inputs
@assert firstindex(agepoints) === 1
@assert firstindex(Tpoints) === 1
function MCMC_varkinetics(data::NamedTuple, model::NamedTuple, boundary::Boundary{T}, constraint::Constraint{T}=Constraint(T), detail::DetailInterval{T}=DetailInterval(T)) where T <: AbstractFloat
# Process inputs
halfwidth = T.(data.halfwidth)::DenseVector{T}
crystAge = T.(data.crystAge)::DenseVector{T}
HeAge = T.(data.HeAge)::DenseVector{T}
Expand All @@ -281,26 +287,36 @@
Th = T.(data.Th)::DenseVector{T}
Sm = (haskey(data, :Sm) ? T.(data.Th) : zeros(T, size(U)))::DenseVector{T}
nsteps = (haskey(model, :nsteps) ? model.nsteps : 10^6)::Int
maxpoints = (haskey(model, :maxpoints) ? model.maxpoints : 50)::Int
minpoints = (haskey(model, :minpoints) ? model.minpoints : 1)::Int
maxpoints = (haskey(model, :maxpoints) ? model.maxpoints : 50)::Int
npoints = (haskey(model, :npoints) ? model.npoints : minpoints)::Int
totalpoints = maxpoints + boundary.npoints + constraint.npoints::Int
simplified = (haskey(model, :simplified) ? model.simplified : false)::Bool
boundarytype = (haskey(model, :boundarytype) ? model.boundarytype : :hard)::Symbol
dynamicjumping = (haskey(model, :dynamicjumping) ? model.dynamicjumping : false)::Bool
dTmax = T(haskey(model, :dTmax) ? model.dTmax : 10.)::T
dTmax_sigma = T(haskey(model, :dTmax_sigma) ? model.dTmax_sigma : 5.)::T
dTmax = T(haskey(model, :dTmax) ? model.dTmax : 10)::T
dTmax_sigma = T(haskey(model, :dTmax_sigma) ? model.dTmax_sigma : 5)::T
agesteps = T.(model.agesteps)::DenseVector{T}
tsteps = T.(model.tsteps)::DenseVector{T}
tinit = T(model.tinit)::T
tnow = zero(T) # It's always time zero today
Tinit = T(model.Tinit)::T
Tnow = T(model.Tnow)::T
Tr = T(haskey(model, :Tr) ? model.Tr : (Tinit+Tnow)/2)::T
dt = T(model.dt)::T
dr = T(model.dr)::T
σmodel = T(model.σmodel)::T
σannealing = T(model.σannealing)::T
λannealing = T(model.λannealing)::T

# Arrays to hold all t and T points (up to npoints=maxpoints)
agepoints = zeros(T, maxpoints)
Tpoints = zeros(T, maxpoints)

# Fill some intermediate points to give the MCMC something to work with
agepoints[1:npoints] .= range(tnow, tinit, length=npoints)
Tpoints[1:npoints] .= Tr # Degrees C

# Calculate number of boundary and unconformity points and allocate buffer for interpolating
agepointbuffer = similar(agepoints, totalpoints)::DenseVector{T}
Tpointbuffer = similar(agepoints, totalpoints)::DenseVector{T}
Expand Down
27 changes: 10 additions & 17 deletions test/testcomplete.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,8 @@ model = (
tinitMax = 4000.0, # Ma -- forbid anything older than this
minpoints = 1, # Minimum allowed number of t-T points
maxpoints = 40, # Maximum allowed number of t-T points
npoints = 5, # Initial number of t-T points
Tr = 250., # Residence temperature of initial proposal
simplified = false, # Prefer simpler tT paths?
boundarytype = :reflecting, # Reflecting boundary conditions
# Model uncertainty is not well known (depends on annealing parameters,
Expand Down Expand Up @@ -66,18 +68,9 @@ unconf = Constraint()

## --- Invert for maximum likelihood t-T path

# This is where the "transdimensional" part comes in
agepoints = Array{Float64}(undef, model.maxpoints+1) # Array of fixed size to hold all optional age points
Tpoints = Array{Float64}(undef, model.maxpoints+1) # Array of fixed size to hold all optional age points
# Fill some intermediate points to give the MCMC something to work with
Tr = 250 # Residence temperature
npoints = 5
agepoints[1:npoints] .= (model.tinit/30,model.tinit/4,model.tinit/2,model.tinit-model.tinit/4,model.tinit-model.tinit/30) # Ma
Tpoints[1:npoints] .= Tr # Degrees C

# Run Markov Chain
@time "Compiling MCMC" MCMC(data, model, npoints, agepoints, Tpoints, boundary, unconf)
@time "Running MCMC" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, npoints, agepoints, Tpoints, boundary, unconf)
@time "Compiling MCMC" MCMC(data, model, boundary, unconf)
@time "Running MCMC" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, boundary, unconf)

@test isa(Tpointdist, AbstractMatrix)
@test maximum(Tpointdist) <= model.Tinit
Expand Down Expand Up @@ -115,7 +108,7 @@ detail = DetailInterval(
agemax = 541, # Oldest end of detail interval
minpoints = 5, # Minimum number of points in detail interval
)
@time "MCMC with Detail interval" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, npoints, agepoints, Tpoints, boundary, unconf, detail)
@time "MCMC with Detail interval" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, boundary, unconf, detail)

@test isa(Tpointdist, AbstractMatrix)
@test maximum(Tpointdist) <= model.Tinit
Expand Down Expand Up @@ -151,7 +144,7 @@ llmean = mean(@view(lldist[model.burnin:end]))
model = (model...,
dynamicjumping=true
)
@time "MCMC with Detail interval & dynamicjumping" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, npoints, agepoints, Tpoints, boundary, unconf, detail)
@time "MCMC with Detail interval & dynamicjumping" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, boundary, unconf, detail)

@test isa(Tpointdist, AbstractMatrix)
@test maximum(Tpointdist) <= model.Tinit
Expand Down Expand Up @@ -186,8 +179,8 @@ llmean = mean(@view(lldist[model.burnin:end]))
## --- As above, but with variable kinetic parameters

# Run Markov Chain
@time "Compiling MCMC_varkinetics" MCMC_varkinetics(data, model, npoints, agepoints, Tpoints, boundary, unconf)
@time "Running MCMC_varkinetics" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, npoints, agepoints, Tpoints, boundary, unconf)
@time "Compiling MCMC_varkinetics" MCMC_varkinetics(data, model, boundary, unconf)
@time "Running MCMC_varkinetics" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, boundary, unconf)

@test isa(Tpointdist, AbstractMatrix)
@test maximum(Tpointdist) <= model.Tinit
Expand Down Expand Up @@ -225,7 +218,7 @@ detail = DetailInterval(
agemax = 541, # Oldest end of detail interval
minpoints = 5, # Minimum number of points in detail interval
)
@time "MCMC_varkinetics with Detail interval" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC_varkinetics(data, model, npoints, agepoints, Tpoints, boundary, unconf, detail)
@time "MCMC_varkinetics with Detail interval" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC_varkinetics(data, model, boundary, unconf, detail)

@test isa(Tpointdist, AbstractMatrix)
@test maximum(Tpointdist) <= model.Tinit
Expand Down Expand Up @@ -261,7 +254,7 @@ llmean = mean(@view(lldist[model.burnin:end]))
model = (model...,
dynamicjumping=true
)
@time "MCMC with Detail interval & dynamicjumping" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, npoints, agepoints, Tpoints, boundary, unconf, detail)
@time "MCMC with Detail interval & dynamicjumping" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, boundary, unconf, detail)

@test isa(Tpointdist, AbstractMatrix)
@test maximum(Tpointdist) <= model.Tinit
Expand Down
19 changes: 6 additions & 13 deletions test/testcompletezrn.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,8 @@ model = (
tinitMax = 4000.0, # Ma -- forbid anything older than this
minpoints = 1, # Minimum allowed number of t-T points
maxpoints = 40, # Maximum allowed number of t-T points
npoints = 5, # Number of initial t-T points
Tr = 250., # Residence temperature in initial proposal
simplified = false, # Prefer simpler tT paths?
boundarytype = :reflecting, # Reflecting boundary conditions
# Model uncertainty is not well known (depends on annealing parameters,
Expand Down Expand Up @@ -65,18 +67,9 @@ constraint = Constraint()

## --- Invert for maximum likelihood t-T path

# This is where the "transdimensional" part comes in
agepoints = Array{Float64}(undef, model.maxpoints+1) # Array of fixed size to hold all optional age points
Tpoints = Array{Float64}(undef, model.maxpoints+1) # Array of fixed size to hold all optional age points
# Fill some intermediate points to give the MCMC something to work with
Tr = 250 # Residence temperature
npoints = 5
agepoints[1:npoints] .= (model.tinit/30,model.tinit/4,model.tinit/2,model.tinit-model.tinit/4,model.tinit-model.tinit/30) # Ma
Tpoints[1:npoints] .= Tr # Degrees C

# Run Markov Chain
@time "Compiling MCMC" MCMC(data, model, npoints, agepoints, Tpoints, boundary, constraint)
@time "Running MCMC" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, npoints, agepoints, Tpoints, boundary, constraint)
@time "Compiling MCMC" MCMC(data, model, boundary, constraint)
@time "Running MCMC" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, boundary, constraint)

@test isa(Tpointdist, AbstractMatrix)
@test maximum(Tpointdist) <= model.Tinit
Expand Down Expand Up @@ -114,7 +107,7 @@ detail = DetailInterval(
agemax = 541, # Oldest end of detail interval
minpoints = 5, # Minimum number of points in detail interval
)
@time "MCMC with Detail interval" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, npoints, agepoints, Tpoints, boundary, constraint, detail)
@time "MCMC with Detail interval" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, boundary, constraint, detail)

@test isa(Tpointdist, AbstractMatrix)
@test maximum(Tpointdist) <= model.Tinit
Expand Down Expand Up @@ -150,7 +143,7 @@ llmean = mean(@view(lldist[model.burnin:end]))
model = (model...,
dynamicjumping=true
)
@time "MCMC with Detail interval & dynamicjumping" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, npoints, agepoints, Tpoints, boundary, constraint, detail)
@time "MCMC with Detail interval & dynamicjumping" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, boundary, constraint, detail)

@test isa(Tpointdist, AbstractMatrix)
@test maximum(Tpointdist) <= model.Tinit
Expand Down

2 comments on commit 9a442f3

@brenhinkeller
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register
Release notes:

  • Move npoints, agepoints, and Tpoints inside MCMC function (breaking)

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/118777

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.10.0 -m "<description of version>" 9a442f3b6881342b6f845e4ec48fd10ac1a009ab
git push origin v0.10.0

Please sign in to comment.