Skip to content

Commit

Permalink
Incorporate reheating rate limit into log likelihood
Browse files Browse the repository at this point in the history
  • Loading branch information
@brenhinkeller
brenhinkeller committed Apr 21, 2024
1 parent b0a8d19 commit 75e92a6
Show file tree
Hide file tree
Showing 2 changed files with 50 additions and 63 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.8.0"
version = "0.9.0"

[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Expand Down
111 changes: 49 additions & 62 deletions src/inversion.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,8 @@
totalpoints = maxpoints + boundary.npoints + unconf.npoints::Int
simplified = (haskey(model, :simplified) ? model.simplified : false)::Bool
dynamicjumping = (haskey(model, :dynamicjumping) ? model.dynamicjumping : false)::Bool
dTmax = T(model.dTmax)::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
Expand Down Expand Up @@ -128,13 +129,10 @@
σ = sqrt.(HeAge_sigma.^2 .+ σmodel^2)

# Log-likelihood for initial proposal
ll = normpdf_ll(HeAge, σₐ, calcHeAges)
if simplified
ll -= log(npoints)
end
llna = llnaₚ = diff_ll(Tsteps, dTmax, dTmax_sigma) + (simplified ? -log(npoints) : zero(T))
ll = llₚ = normpdf_ll(HeAge, σₐ, calcHeAges) + llna

# Variables to hold proposals
llₚ = ll
npointsₚ = npoints
agepointsₚ = copy(agepoints)::DenseVector{T}
Tpointsₚ = copy(Tpoints)::DenseVector{T}
Expand Down Expand Up @@ -235,7 +233,7 @@
temperatures = collectto!(Tpointbuffer, view(Tpointsₚ, 1:npointsₚ), boundary.Tpointsₚ, unconf.Tpointsₚ)::StridedVector{T}
linterp1s!(Tsteps, knot_index, ages, temperatures, agesteps)

(maxdiff(Tsteps) > dTmax) && @goto restart
# (maxdiff(Tsteps) > dTmax) && @goto restart

# Calculate model ages for each grain
if any(tzr)
Expand All @@ -258,13 +256,11 @@
end

# Calculate log likelihood of proposal
llnaₚ = diff_ll(Tsteps, dTmax, dTmax_sigma)
simplified && (llnaₚ += -log(npointsₚ))
σₐ .= simannealsigma.(n, HeAge_sigma, σmodel, σannealing, λannealing)
llₚ = normpdf_ll(HeAge, σₐ, calcHeAgesₚ)
llₗ = normpdf_ll(HeAge, σₐ, calcHeAges) # Recalulate last one too with new σₐ
if simplified # slightly penalize more complex t-T paths
llₚ -= log(npointsₚ)
llₗ -= log(npoints)
end
llₚ = normpdf_ll(HeAge, σₐ, calcHeAgesₚ) + llnaₚ
llₗ = normpdf_ll(HeAge, σₐ, calcHeAges) + llna # Recalulate last one too with new σₐ

# Accept or reject proposal based on likelihood
if log(rand()) < (llₚ - llₗ)
Expand All @@ -281,6 +277,7 @@

# Update the currently accepted proposal
ll = llₚ
llna = llnaₚ
npoints = npointsₚ
copyto!(agepoints, 1, agepointsₚ, 1, npoints)
copyto!(Tpoints, 1, Tpointsₚ, 1, npoints)
Expand All @@ -294,7 +291,7 @@
end

# Record results for analysis and troubleshooting
lldist[n] = normpdf_ll(HeAge, σ, calcHeAges) # Recalculated to constant baseline
lldist[n] = llna + normpdf_ll(HeAge, σ, calcHeAges) # Recalculated to constant baseline
ndist[n] = npoints # distribution of # of points
σⱼtdist[n] = σⱼt
σⱼTdist[n] = σⱼT
Expand Down Expand Up @@ -327,7 +324,8 @@
totalpoints = maxpoints + boundary.npoints + unconf.npoints::Int
simplified = (haskey(model, :simplified) ? model.simplified : false)::Bool
dynamicjumping = (haskey(model, :dynamicjumping) ? model.dynamicjumping : false)::Bool
dTmax = T(model.dTmax)::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
Expand Down Expand Up @@ -401,13 +399,10 @@
σ = sqrt.(HeAge_sigma.^2 .+ σmodel^2)

# Log-likelihood for initial proposal
ll = normpdf_ll(HeAge, σₐ, calcHeAges)
if simplified
ll -= log(npoints)
end
llna = llnaₚ = diff_ll(Tsteps, dTmax, dTmax_sigma) + (simplified ? -log(npoints) : zero(T))
ll = llₚ = normpdf_ll(HeAge, σₐ, calcHeAges) + llna

# Variables to hold proposals
llₚ = ll
npointsₚ = npoints
agepointsₚ = copy(agepoints)::DenseVector{T}
Tpointsₚ = copy(Tpoints)::DenseVector{T}
Expand Down Expand Up @@ -509,7 +504,7 @@
temperatures = collectto!(Tpointbuffer, view(Tpointsₚ, 1:npointsₚ), boundary.Tpointsₚ, unconf.Tpointsₚ)::StridedVector{T}
linterp1s!(Tsteps, knot_index, ages, temperatures, agesteps)

(maxdiff(Tsteps) > dTmax) && @goto restart
# (maxdiff(Tsteps) > dTmax) && @goto restart
(pointsininterval(agepointsₚ, npointsₚ, detail.agemin, detail.agemax, dt) < enoughpoints) && @goto restart

# Calculate model ages for each grain
Expand All @@ -533,13 +528,11 @@
end

# Calculate log likelihood of proposal
llnaₚ = diff_ll(Tsteps, dTmax, dTmax_sigma)
simplified && (llnaₚ += -log(npointsₚ))
σₐ .= simannealsigma.(n, HeAge_sigma, σmodel, σannealing, λannealing)
llₚ = normpdf_ll(HeAge, σₐ, calcHeAgesₚ)
llₗ = normpdf_ll(HeAge, σₐ, calcHeAges) # Recalulate last one too with new σₐ
if simplified # slightly penalize more complex t-T paths
llₚ -= log(npointsₚ)
llₗ -= log(npoints)
end
llₚ = normpdf_ll(HeAge, σₐ, calcHeAgesₚ) + llnaₚ
llₗ = normpdf_ll(HeAge, σₐ, calcHeAges) + llna # Recalulate last one too with new σₐ

# Accept or reject proposal based on likelihood
if log(rand()) < (llₚ - llₗ)
Expand All @@ -556,6 +549,7 @@

# Update the currently accepted proposal
ll = llₚ
llna = llnaₚ
npoints = npointsₚ
copyto!(agepoints, 1, agepointsₚ, 1, npoints)
copyto!(Tpoints, 1, Tpointsₚ, 1, npoints)
Expand All @@ -569,7 +563,7 @@
end

# Record results for analysis and troubleshooting
lldist[n] = normpdf_ll(HeAge, σ, calcHeAges) # Recalculated to constant baseline
lldist[n] = llna + normpdf_ll(HeAge, σ, calcHeAges) # Recalculated to constant baseline
ndist[n] = npoints # distribution of # of points
σⱼtdist[n] = σⱼt
σⱼTdist[n] = σⱼT
Expand Down Expand Up @@ -604,7 +598,8 @@
totalpoints = maxpoints + boundary.npoints + unconf.npoints::Int
simplified = (haskey(model, :simplified) ? model.simplified : false)::Bool
dynamicjumping = (haskey(model, :dynamicjumping) ? model.dynamicjumping : false)::Bool
dTmax = T(model.dTmax)::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
Expand Down Expand Up @@ -678,13 +673,10 @@
σ = sqrt.(HeAge_sigma.^2 .+ σmodel^2)

# Log-likelihood for initial proposal
ll = normpdf_ll(HeAge, σₐ, calcHeAges)
if simplified
ll -= log(npoints)
end
llna = llnaₚ = diff_ll(Tsteps, dTmax, dTmax_sigma) + loglikelihood(admₚ, adm₀) + loglikelihood(zdmₚ, zdm₀) + (simplified ? -log(npoints) : zero(T))
ll = llₚ = normpdf_ll(HeAge, σₐ, calcHeAges) + llna

# Variables to hold proposals
llₚ = ll
npointsₚ = npoints
agepointsₚ = copy(agepoints)::DenseVector{T}
Tpointsₚ = copy(Tpoints)::DenseVector{T}
Expand Down Expand Up @@ -814,7 +806,7 @@
temperatures = collectto!(Tpointbuffer, view(Tpointsₚ, 1:npointsₚ), boundary.Tpointsₚ, unconf.Tpointsₚ)::StridedVector{T}
linterp1s!(Tsteps, knot_index, ages, temperatures, agesteps)

(maxdiff(Tsteps) > dTmax) && @goto restart
# (maxdiff(Tsteps) > dTmax) && @goto restart
(pointsininterval(agepointsₚ, npointsₚ, detail.agemin, detail.agemax, dt) < enoughpoints) && @goto restart

# Calculate model ages for each grain
Expand All @@ -838,15 +830,12 @@
end

# Calculate log likelihood of proposal
llnaₚ = diff_ll(Tsteps, dTmax, dTmax_sigma)
llnaₚ += loglikelihood(admₚ, adm₀) + loglikelihood(zdmₚ, zdm₀)
simplified && (llnaₚ += -log(npointsₚ))
σₐ .= simannealsigma.(n, HeAge_sigma, σmodel, σannealing, λannealing)
llₚ = normpdf_ll(HeAge, σₐ, calcHeAgesₚ)
llₗ = normpdf_ll(HeAge, σₐ, calcHeAges) # Recalulate last one too with new σₐ
if simplified # slightly penalize more complex t-T paths
llₚ -= log(npointsₚ)
llₗ -= log(npoints)
end
llₚ += loglikelihood(admₚ, adm₀) + loglikelihood(zdmₚ, zdm₀)
llₗ += loglikelihood(adm, adm₀) + loglikelihood(zdm, zdm₀)
llₚ = normpdf_ll(HeAge, σₐ, calcHeAgesₚ) + llnaₚ
llₗ = normpdf_ll(HeAge, σₐ, calcHeAges) + llna # Recalulate last one too with new σₐ

# Accept or reject proposal based on likelihood
if log(rand()) < (llₚ - llₗ)
Expand All @@ -862,9 +851,10 @@
end

# Update the currently accepted proposal
ll = llₚ
adm = admₚ
zdm = zdmₚ
ll = llₚ
llna = llnaₚ
npoints = npointsₚ
copyto!(agepoints, 1, agepointsₚ, 1, npoints)
copyto!(Tpoints, 1, Tpointsₚ, 1, npoints)
Expand All @@ -878,7 +868,7 @@
end

# Record results for analysis and troubleshooting
lldist[n] = normpdf_ll(HeAge, σ, calcHeAges) # Recalculated to constant baseline
lldist[n] = llna + normpdf_ll(HeAge, σ, calcHeAges) # Recalculated to constant baseline
ndist[n] = npoints # distribution of # of points
σⱼtdist[n] = σⱼt
σⱼTdist[n] = σⱼT
Expand Down Expand Up @@ -911,7 +901,8 @@
totalpoints = maxpoints + boundary.npoints + unconf.npoints::Int
simplified = (haskey(model, :simplified) ? model.simplified : false)::Bool
dynamicjumping = (haskey(model, :dynamicjumping) ? model.dynamicjumping : false)::Bool
dTmax = T(model.dTmax)::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
Expand Down Expand Up @@ -985,13 +976,10 @@
σ = sqrt.(HeAge_sigma.^2 .+ σmodel^2)

# Log-likelihood for initial proposal
ll = normpdf_ll(HeAge, σₐ, calcHeAges)
if simplified
ll -= log(npoints)
end
llna = llnaₚ = diff_ll(Tsteps, dTmax, dTmax_sigma) + loglikelihood(admₚ, adm₀) + loglikelihood(zdmₚ, zdm₀) + (simplified ? -log(npoints) : zero(T))
ll = llₚ = normpdf_ll(HeAge, σₐ, calcHeAges) + llna

# Variables to hold proposals
llₚ = ll
npointsₚ = npoints
agepointsₚ = copy(agepoints)::DenseVector{T}
Tpointsₚ = copy(Tpoints)::DenseVector{T}
Expand Down Expand Up @@ -1120,7 +1108,7 @@
temperatures = collectto!(Tpointbuffer, view(Tpointsₚ, 1:npointsₚ), boundary.Tpointsₚ, unconf.Tpointsₚ)::StridedVector{T}
linterp1s!(Tsteps, knot_index, ages, temperatures, agesteps)

(maxdiff(Tsteps) > dTmax) && @goto restart
# (maxdiff(Tsteps) > dTmax) && @goto restart

# Calculate model ages for each grain
if any(tzr)
Expand All @@ -1143,15 +1131,13 @@
end

# Calculate log likelihood of proposal
llnaₚ = diff_ll(Tsteps, dTmax, dTmax_sigma)
llnaₚ += loglikelihood(admₚ, adm₀) + loglikelihood(zdmₚ, zdm₀)
simplified && (llnaₚ += -log(npointsₚ))
σₐ .= simannealsigma.(n, HeAge_sigma, σmodel, σannealing, λannealing)
llₚ = normpdf_ll(HeAge, σₐ, calcHeAgesₚ)
llₗ = normpdf_ll(HeAge, σₐ, calcHeAges) # Recalulate last one too with new σₐ
if simplified # slightly penalize more complex t-T paths
llₚ -= log(npointsₚ)
llₗ -= log(npoints)
end
llₚ += loglikelihood(admₚ, adm₀) + loglikelihood(zdmₚ, zdm₀)
llₗ += loglikelihood(adm, adm₀) + loglikelihood(zdm, zdm₀)
llₚ = normpdf_ll(HeAge, σₐ, calcHeAgesₚ) + llnaₚ
llₗ = normpdf_ll(HeAge, σₐ, calcHeAges) + llna # Recalulate last one too with new σₐ


# Accept or reject proposal based on likelihood
if log(rand()) < (llₚ - llₗ)
Expand All @@ -1167,9 +1153,10 @@
end

# Update the currently accepted proposal
ll = llₚ
adm = admₚ
zdm = zdmₚ
ll = llₚ
llna = llnaₚ
npoints = npointsₚ
copyto!(agepoints, agepointsₚ)
copyto!(Tpoints, Tpointsₚ)
Expand All @@ -1183,7 +1170,7 @@
end

# Record results for analysis and troubleshooting
lldist[n] = normpdf_ll(HeAge, σ, calcHeAges) # Recalculated to constant baseline
lldist[n] = llna + normpdf_ll(HeAge, σ, calcHeAges) # Recalculated to constant baseline
ndist[n] = npoints # distribution of # of points
σⱼtdist[n] = σⱼt
σⱼTdist[n] = σⱼT
Expand Down

2 comments on commit 75e92a6

@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:

  • Incorporate reheating rate limit into log likelihood

@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/105404

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.9.0 -m "<description of version>" 75e92a61a953fa5c5789780fa7cdf5d15bed1395
git push origin v0.9.0

Please sign in to comment.