Skip to content

Commit 46992e5

Browse files
brenhinkellerbrenhinkeller
committed
Allow switching method of handling one-sided constraints, using precalculated mean and std of BilinearExponential 
1 parent 93c71ad commit 46992e5

File tree

4 files changed

+154
-74
lines changed

4 files changed

+154
-74
lines changed

‎src/Objects.jl

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@
2121
inputSigmaLevel::Int
2222
Age_Unit::String
2323
Height_Unit::String
24+
Sidedness_Method::Symbol
2425
end
2526

2627
function ChronAgeData(nSamples::Integer)
@@ -43,6 +44,7 @@
4344
2, # i.e., are the data files 1-sigma or 2-sigma
4445
"Ma",
4546
"m",
47+
:cdf,
4648
)
4749
end
4850

@@ -55,6 +57,7 @@
5557
Chronometer::NTuple{N, Symbol}
5658
Age_Unit::String
5759
Height_Unit::String
60+
Sidedness_Method::Symbol
5861
end
5962

6063
function GeneralAgeData(nSamples::Integer)
@@ -67,9 +70,27 @@
6770
ntuple(i->:Chronometer, nSamples), # Age Types (e.g., :UPb or :ArAr)
6871
"Ma",
6972
"m",
73+
:cdf,
7074
)
7175
end
7276

77+
# A set of types to allow dispatch to different age sidedness methods
78+
abstract type Sidedness end
79+
struct CDFSidedness{T<:Real} <: Sidedness
80+
s::Vector{T}
81+
end
82+
struct FastSidedness{T<:Real} <: Sidedness
83+
s::Vector{T}
84+
end
85+
Base.getindex(A::Sidedness, args...) = getindex(A.s, args...)
86+
Base.setindex!(A::Sidedness, args...) = setindex!(A.s, args...)
87+
Base.eachindex(A::Sidedness) = eachindex(A.s)
88+
Base.size(A::Sidedness, args...) = size(A.s, args...)
89+
Base.length(A::Sidedness) = length(A.s)
90+
Base.ndims(A::Sidedness) = ndims(A.s)
91+
Base.ndims(::Type{<:Sidedness}) = 1
92+
Base.copyto!(A::Sidedness, args...) = copyto!(A.s, args...)
93+
7394
# One-sigma systematic uncertainty
7495
mutable struct SystematicUncertainty
7596
UPb::Float64

‎src/StratMetropolis.jl

Lines changed: 48 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -68,9 +68,16 @@
6868
(a,b) = hcat(fill!(similar(Height), 1), Height) \ Age
6969
model_ages = a .+ b .* collect(model_heights)
7070

71+
# Select sidedness method
72+
sidedness = if smpl.Sidedness_Method === :fast || all(iszero, smpl.Age_Sidedness)
73+
FastSidedness(Age_Sidedness)
74+
else
75+
CDFSidedness(Age_Sidedness)
76+
end
77+
7178
# Run the Markov chain
7279
ages = Normal.(Age, Age_sigma)
73-
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve, Chronometer, systematic)
80+
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve, Chronometer, systematic)
7481

7582
# Crop the result
7683
agedist = agedist[active_height_t,:]
@@ -124,9 +131,16 @@
124131
(a,b) = hcat(fill!(similar(Height), 1), Height) \ Age
125132
model_ages = a .+ b .* collect(model_heights)
126133

134+
# Select sidedness method
135+
sidedness = if smpl.Sidedness_Method === :fast || all(iszero, smpl.Age_Sidedness)
136+
FastSidedness(Age_Sidedness)
137+
else
138+
CDFSidedness(Age_Sidedness)
139+
end
140+
127141
# Run the Markov chain
128142
ages = Normal.(Age, Age_sigma)
129-
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
143+
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
130144

131145
# Crop the result
132146
agedist = agedist[active_height_t,:]
@@ -210,9 +224,16 @@
210224
(a,b) = hcat(fill!(similar(Height), 1), Height) \ Age
211225
model_ages = a .+ b .* collect(model_heights)
212226

227+
# Select sidedness method
228+
sidedness = if smpl.Sidedness_Method === :fast || all(iszero, smpl.Age_Sidedness)
229+
FastSidedness(Age_Sidedness)
230+
else
231+
CDFSidedness(Age_Sidedness)
232+
end
233+
213234
# Run the Markov chain
214235
ages = BilinearExponential.(eachcol(p))
215-
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve, Chronometer, systematic)
236+
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve, Chronometer, systematic)
216237

217238
# Crop the result
218239
agedist = agedist[active_height_t,:]
@@ -270,9 +291,16 @@
270291
(a,b) = hcat(fill!(similar(Height), 1), Height) \ Age
271292
model_ages = a .+ b .* collect(model_heights)
272293

294+
# Select sidedness method
295+
sidedness = if smpl.Sidedness_Method === :fast || all(iszero, smpl.Age_Sidedness)
296+
FastSidedness(Age_Sidedness)
297+
else
298+
CDFSidedness(Age_Sidedness)
299+
end
300+
273301
# Run the Markov chain
274302
ages = BilinearExponential.(eachcol(p))
275-
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
303+
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
276304

277305
# Crop the result
278306
agedist = agedist[active_height_t,:]
@@ -355,9 +383,16 @@
355383
(a,b) = hcat(fill!(similar(Height), 1), Height) \ Age
356384
model_ages = a .+ b .* collect(model_heights)
357385

386+
# Select sidedness method
387+
sidedness = if smpl.Sidedness_Method === :fast || all(iszero, smpl.Age_Sidedness)
388+
FastSidedness(Age_Sidedness)
389+
else
390+
CDFSidedness(Age_Sidedness)
391+
end
392+
358393
# Run the Markov chain
359394
ages = Radiocarbon.(Age, Age_sigma, (collect(c) for c in eachcol(p)))
360-
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
395+
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
361396

362397
# Crop the result
363398
agedist = agedist[active_height_t,:]
@@ -419,9 +454,16 @@
419454
(a,b) = hcat(fill!(similar(Height), 1), Height) \ Age
420455
model_ages = a .+ b .* collect(model_heights)
421456

457+
# Select sidedness method
458+
sidedness = if smpl.Sidedness_Method === :fast || all(iszero, smpl.Age_Sidedness)
459+
FastSidedness(Age_Sidedness)
460+
else
461+
CDFSidedness(Age_Sidedness)
462+
end
463+
422464
# Run the Markov chain
423465
ages = Radiocarbon.(Age, Age_sigma, (collect(c) for c in eachcol(p)))
424-
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
466+
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
425467

426468
# Crop the result
427469
agedist = agedist[active_height_t,:]

0 commit comments

Comments
 (0)