Skip to content

Commit

Permalink
Make σⱼt and σⱼT point-specific
Browse files Browse the repository at this point in the history
  • Loading branch information
@brenhinkeller
brenhinkeller committed Nov 7, 2024
1 parent b4c040f commit be07dac
Show file tree
Hide file tree
Showing 4 changed files with 87 additions and 57 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.10.0"
version = "0.10.1"

[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Expand Down
62 changes: 34 additions & 28 deletions src/inversion.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,8 +135,8 @@

# Standard deviations of Gaussian proposal ("jumping") distributions
# for temperature and time
σⱼt = tinit/60
σⱼT = Tinit/60
σⱼt = fill((tinit-tnow)/60, maxpoints)
σⱼT = fill((Tinit-Tnow)/60, maxpoints)
k = 1 # Index of chosen t-T point

# Proposal probabilities (must sum to 1)
Expand All @@ -161,23 +161,26 @@
copyto!(constraint.Tpointsₚ, constraint.Tpoints)
copyto!(boundary.Tpointsₚ, boundary.Tpoints)

# Adjust the proposal
# Randomly choose an option and point (if applicable) to adjust
r = rand()
k = ceil(Int, rand()*npoints)

# Adjust the proposal
if r < p_move
# Move one t-T point
k = ceil(Int, rand() * npoints)
movepoint!(agepointsₚ, Tpointsₚ, k, tnow+dt, tinit-dt, Tnow, Tinit, σⱼt, σⱼT, boundarytype)
movepoint!(agepointsₚ, Tpointsₚ, k, tnow+dt, tinit-dt, Tnow, Tinit, σⱼt[k], σⱼT[k], boundarytype)

elseif (r < p_move+p_birth) && (npoints < maxpoints)
# Birth: add a new model point
npointsₚ = npoints + 1
agepointsₚ[npointsₚ] = 0 + rand()*(tinit-0)
Tpointsₚ[npointsₚ] = Tnow + rand()*(Tinit-Tnow)
k = npointsₚ = npoints + 1
agepointsₚ[k] = tnow + rand()*(tinit-tnow)
Tpointsₚ[k] = Tnow + rand()*(Tinit-Tnow)
σⱼt[k] = (tinit-Tnow)/60
σⱼT[k] = (Tinit-Tnow)/60

elseif (r < p_move+p_birth+p_death) && (r >= p_move+p_birth) && (npoints > max(minpoints, detail.minpoints))
# Death: remove a model point
npointsₚ = npoints - 1 # Delete last point in array from proposal
k = ceil(Int, rand()*npoints) # Choose point to delete
npointsₚ = npoints - 1
agepointsₚ[k] = agepointsₚ[npoints]
Tpointsₚ[k] = Tpointsₚ[npoints]

Expand Down Expand Up @@ -237,10 +240,10 @@
# Update jumping distribution based on size of current accepted p_move
if dynamicjumping && r < p_move
if agepointsₚ[k] != agepoints[k]
σⱼt = max(ℯ * abs(agepointsₚ[k] - agepoints[k]), dt)
σⱼt[k] = ℯ * abs(agepointsₚ[k] - agepoints[k])
end
if Tpointsₚ[k] != Tpoints[k]
σⱼT = max(ℯ * abs(Tpointsₚ[k] - Tpoints[k]), (Tinit-Tnow)/100)
σⱼT[k] = ℯ * abs(Tpointsₚ[k] - Tpoints[k])
end
end

Expand All @@ -262,8 +265,8 @@
# Record results for analysis and troubleshooting
lldist[n] = llna + normpdf_ll(HeAge, σ, calcHeAges) # Recalculated to constant baseline
ndist[n] = npoints # distribution of # of points
σⱼtdist[n] = σⱼt
σⱼTdist[n] = σⱼT
σⱼtdist[n] = σⱼt[k]
σⱼTdist[n] = σⱼT[k]
HeAgedist[:,n] .= calcHeAges # distribution of He ages

# This is the actual output we want -- the distribution of t-T paths (t path is always identical)
Expand Down Expand Up @@ -400,8 +403,8 @@

# Standard deviations of Gaussian proposal ("jumping") distributions
# for temperature and time
σⱼt = tinit/60
σⱼT = Tinit/60
σⱼt = fill((tinit-tnow)/60, maxpoints)
σⱼT = fill((Tinit-Tnow)/60, maxpoints)
k = 1 # Index of chosen t-T point

# Proposal probabilities (must sum to 1)
Expand Down Expand Up @@ -429,23 +432,26 @@
copyto!(constraint.Tpointsₚ, constraint.Tpoints)
copyto!(boundary.Tpointsₚ, boundary.Tpoints)

# Adjust the proposal
# Randomly choose an option and point (if applicable) to adjust
r = rand()
k = ceil(Int, rand()*npoints)

# Adjust the proposal
if r < p_move
# Move one t-T point
k = ceil(Int, rand() * npoints)
movepoint!(agepointsₚ, Tpointsₚ, k, tnow+dt, tinit-dt, Tnow, Tinit, σⱼt, σⱼT, boundarytype)
movepoint!(agepointsₚ, Tpointsₚ, k, tnow+dt, tinit-dt, Tnow, Tinit, σⱼt[k], σⱼT[k], boundarytype)

elseif (r < p_move+p_birth) && (npoints < maxpoints)
# Birth: add a new model point
npointsₚ = npoints + 1
agepointsₚ[npointsₚ] = 0 + rand()*(tinit-0)
Tpointsₚ[npointsₚ] = Tnow + rand()*(Tinit-Tnow)
k = npointsₚ = npoints + 1
agepointsₚ[k] = tnow + rand()*(tinit-tnow)
Tpointsₚ[k] = Tnow + rand()*(Tinit-Tnow)
σⱼt[k] = (tinit-Tnow)/60
σⱼT[k] = (Tinit-Tnow)/60

elseif (r < p_move+p_birth+p_death) && (r >= p_move+p_birth) && (npoints > max(minpoints, detail.minpoints))
# Death: remove a model point
npointsₚ = npoints - 1 # Delete last point in array from proposal
k = ceil(Int, rand()*npoints) # Choose point to delete
npointsₚ = npoints - 1
agepointsₚ[k] = agepointsₚ[npoints]
Tpointsₚ[k] = Tpointsₚ[npoints]

Expand Down Expand Up @@ -513,10 +519,10 @@
# Update jumping distribution based on size of current accepted p_move
if dynamicjumping && r < p_move
if agepointsₚ[k] != agepoints[k]
σⱼt = max(ℯ * abs(agepointsₚ[k] - agepoints[k]), dt)
σⱼt[k] = ℯ * abs(agepointsₚ[k] - agepoints[k])
end
if Tpointsₚ[k] != Tpoints[k]
σⱼT = max(ℯ * abs(Tpointsₚ[k] - Tpoints[k]), (Tinit-Tnow)/100)
σⱼT[k] = ℯ * abs(Tpointsₚ[k] - Tpoints[k])
end
end

Expand All @@ -540,8 +546,8 @@
# Record results for analysis and troubleshooting
lldist[n] = llna + normpdf_ll(HeAge, σ, calcHeAges) # Recalculated to constant baseline
ndist[n] = npoints # distribution of # of points
σⱼtdist[n] = σⱼt
σⱼTdist[n] = σⱼT
σⱼtdist[n] = σⱼt[k]
σⱼTdist[n] = σⱼT[k]
HeAgedist[:,n] .= calcHeAges # distribution of He ages

# This is the actual output we want -- the distribution of t-T paths (t path is always identical)
Expand Down
67 changes: 41 additions & 26 deletions test/testcomplete.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ data = (
)

# Sort out crystallization ages and start time
map!(x->max(x, model.tinitMax), data.crystAge, data.crystAge)
map!(x->min(x, model.tinitMax), data.crystAge, data.crystAge)
tinit = ceil(maximum(data.crystAge)/model.dt) * model.dt
model = (model...,
tinit = tinit,
Expand Down Expand Up @@ -99,16 +99,18 @@ llmean = mean(@view(lldist[model.burnin:end]))
@test maximum(ndist) <= model.maxpoints
@info "Mean npoints: $(mean(ndist[model.burnin:end]))"


@test mean(σⱼtdist[model.burnin:end]) ≈ model.tinit/60
@info "Mean σⱼₜ: $(mean(σⱼtdist[model.burnin:end]))"

@test mean(σⱼTdist[model.burnin:end]) ≈ model.Tinit/60
@info "Mean σⱼT: $(mean(σⱼTdist[model.burnin:end]))"

## ---
detail = DetailInterval(
agemin = 0, # Youngest end of detail interval
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, boundary, unconf, detail)
## --- As above, but with variable kinetic parameters

# Run Markov Chain
@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 All @@ -125,7 +127,7 @@ abserr = abs(sum(nanmean(HeAgedist[:,model.burnin:end], dims=2) - data.HeAge)/le

@test isa(lldist, AbstractVector)
llmean = mean(@view(lldist[model.burnin:end]))
@test -200 < llmean < 0
@test -300 < llmean < 0
@info "Mean ll: $llmean"

@test isa(acceptancedist, AbstractVector{Bool})
Expand All @@ -137,14 +139,20 @@ llmean = mean(@view(lldist[model.burnin:end]))
@test maximum(ndist) <= model.maxpoints
@info "Mean npoints: $(mean(ndist[model.burnin:end]))"


@test mean(σⱼtdist[model.burnin:end]) ≈ model.tinit/60
@info "Mean σⱼₜ: $(mean(σⱼtdist[model.burnin:end]))"

@test mean(σⱼTdist[model.burnin:end]) ≈ model.Tinit/60
@info "Mean σⱼT: $(mean(σⱼTdist[model.burnin:end]))"

## ---
model = (model...,
dynamicjumping=true
detail = DetailInterval(
agemin = 0, # Youngest end of detail interval
agemax = 541, # Oldest end of detail interval
minpoints = 5, # Minimum number of points in detail interval
)
@time "MCMC with Detail interval & dynamicjumping" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, 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 @@ -173,14 +181,15 @@ llmean = mean(@view(lldist[model.burnin:end]))
@test maximum(ndist) <= model.maxpoints
@info "Mean npoints: $(mean(ndist[model.burnin:end]))"

@test mean(σⱼtdist[model.burnin:end]) ≈ model.tinit/60
@info "Mean σⱼₜ: $(mean(σⱼtdist[model.burnin:end]))"

@test mean(σⱼTdist[model.burnin:end]) ≈ model.Tinit/60
@info "Mean σⱼT: $(mean(σⱼTdist[model.burnin:end]))"

## --- As above, but with variable kinetic parameters
## ---

# Run Markov Chain
@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)
@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 All @@ -197,7 +206,7 @@ abserr = abs(sum(nanmean(HeAgedist[:,model.burnin:end], dims=2) - data.HeAge)/le

@test isa(lldist, AbstractVector)
llmean = mean(@view(lldist[model.burnin:end]))
@test -300 < llmean < 0
@test -200 < llmean < 0
@info "Mean ll: $llmean"

@test isa(acceptancedist, AbstractVector{Bool})
Expand All @@ -209,16 +218,18 @@ llmean = mean(@view(lldist[model.burnin:end]))
@test maximum(ndist) <= model.maxpoints
@info "Mean npoints: $(mean(ndist[model.burnin:end]))"


@test mean(σⱼtdist[model.burnin:end]) ≈ model.tinit/60
@info "Mean σⱼₜ: $(mean(σⱼtdist[model.burnin:end]))"

@test mean(σⱼTdist[model.burnin:end]) ≈ model.Tinit/60
@info "Mean σⱼT: $(mean(σⱼTdist[model.burnin:end]))"

## ---
detail = DetailInterval(
agemin = 0, # Youngest end of detail interval
agemax = 541, # Oldest end of detail interval
minpoints = 5, # Minimum number of points in detail interval
model = (model...,
dynamicjumping=true
)
@time "MCMC_varkinetics with Detail interval" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC_varkinetics(data, model, 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 @@ -247,14 +258,15 @@ llmean = mean(@view(lldist[model.burnin:end]))
@test maximum(ndist) <= model.maxpoints
@info "Mean npoints: $(mean(ndist[model.burnin:end]))"

@test model.dt < mean(σⱼtdist[model.burnin:end]) < model.tinit
@info "Mean σⱼₜ: $(mean(σⱼtdist[model.burnin:end]))"

@test 0 < mean(σⱼTdist[model.burnin:end]) < model.Tinit
@info "Mean σⱼT: $(mean(σⱼTdist[model.burnin:end]))"

## ---
model = (model...,
dynamicjumping=true
)
@time "MCMC with Detail interval & dynamicjumping" (tpointdist, Tpointdist, ndist, HeAgedist, lldist, acceptancedist, σⱼtdist, σⱼTdist) = MCMC(data, model, boundary, unconf, detail)

@time "MCMC_varkinetics with Detail interval & dynamicjumping" (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 @@ -283,5 +295,8 @@ llmean = mean(@view(lldist[model.burnin:end]))
@test maximum(ndist) <= model.maxpoints
@info "Mean npoints: $(mean(ndist[model.burnin:end]))"

@test model.dt < mean(σⱼtdist[model.burnin:end]) < model.tinit
@info "Mean σⱼₜ: $(mean(σⱼtdist[model.burnin:end]))"

@test 0 < mean(σⱼTdist[model.burnin:end]) < model.Tinit
@info "Mean σⱼT: $(mean(σⱼTdist[model.burnin:end]))"
13 changes: 11 additions & 2 deletions test/testcompletezrn.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ data = (
)

# Sort out crystallization ages and start time
map!(x->max(x, model.tinitMax), data.crystAge, data.crystAge)
map!(x->min(x, model.tinitMax), data.crystAge, data.crystAge)
tinit = ceil(maximum(data.crystAge)/model.dt) * model.dt
model = (model...,
tinit = tinit,
Expand Down Expand Up @@ -98,7 +98,10 @@ llmean = mean(@view(lldist[model.burnin:end]))
@test maximum(ndist) <= model.maxpoints
@info "Mean npoints: $(mean(ndist[model.burnin:end]))"

@test mean(σⱼtdist[model.burnin:end]) ≈ model.tinit/60
@info "Mean σⱼₜ: $(mean(σⱼtdist[model.burnin:end]))"

@test mean(σⱼTdist[model.burnin:end]) ≈ model.Tinit/60
@info "Mean σⱼT: $(mean(σⱼTdist[model.burnin:end]))"

## ---
Expand Down Expand Up @@ -136,7 +139,10 @@ llmean = mean(@view(lldist[model.burnin:end]))
@test maximum(ndist) <= model.maxpoints
@info "Mean npoints: $(mean(ndist[model.burnin:end]))"

@test mean(σⱼtdist[model.burnin:end]) ≈ model.tinit/60
@info "Mean σⱼₜ: $(mean(σⱼtdist[model.burnin:end]))"

@test mean(σⱼTdist[model.burnin:end]) ≈ model.Tinit/60
@info "Mean σⱼT: $(mean(σⱼTdist[model.burnin:end]))"

## ---
Expand Down Expand Up @@ -172,5 +178,8 @@ llmean = mean(@view(lldist[model.burnin:end]))
@test maximum(ndist) <= model.maxpoints
@info "Mean npoints: $(mean(ndist[model.burnin:end]))"

@test model.dt < mean(σⱼtdist[model.burnin:end]) < model.tinit
@info "Mean σⱼₜ: $(mean(σⱼtdist[model.burnin:end]))"
@info "Mean σⱼT: $(mean(σⱼTdist[model.burnin:end]))"

@test 0 < mean(σⱼTdist[model.burnin:end]) < model.Tinit
@info "Mean σⱼT: $(mean(σⱼTdist[model.burnin:end]))"

2 comments on commit be07dac

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

  • Make σⱼt and σⱼT point-specific

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

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.1 -m "<description of version>" be07daceb39c90e978a72d4e4ba6cc37b528fa5e
git push origin v0.10.1

Please sign in to comment.