diff --git a/src/metropolis.jl b/src/metropolis.jl index 9d7e5e9..eed75de 100644 --- a/src/metropolis.jl +++ b/src/metropolis.jl @@ -14,21 +14,24 @@ end """ - strictpriors(p::Proposal, record_max_age::Number, climatelimits::Tuple{Number,Number}) + PorewaterDiffusion.strictpriors(p::Proposal, record_max_age::Number, climatelimits::Tuple{Number,Number}, k::Constants) Evalute strict constraints on priors that will automatically reject a proposal with... - Onset date beyond the climate record timespan (in ka) (`record_max_age`) - Nonphysical subglacial thresholds -- melting at lower benthic δ¹⁸O than freezing or values exceeding the record extrema (`climatelimits`). +- Freezing rate is ≤ 0.4dt/dz (to prevent an error in a log-calculation in [`PorewaterDiffusion.boundaryconditions`](@ref)). """ -function strictpriors(p::Proposal, record_max_age::Number, climatelimits::Tuple{Number,Number}) +function strictpriors(p::Proposal, record_max_age::Number, climatelimits::Tuple{Number,Number}, k::Constants) x=true - x *= p.onset <= record_max_age - x *= 0 < p.onset + x &= p.onset <= record_max_age + x &= 0 < p.onset - x *= p.sea2frz < p.frz2mlt - x *= climatelimits[1] < p.sea2frz - x *= p.frz2mlt < climatelimits[2] + x &= p.sea2frz < p.frz2mlt + x &= climatelimits[1] < p.sea2frz + x &= p.frz2mlt < climatelimits[2] + + x &= p.dfrz <= 0.4k.dz/k.dt x end @@ -97,7 +100,7 @@ function porewatermetropolis(p::Proposal, jumpsigma::Proposal, prior::CoreData; ϕ, jumpname, jump = proposaljump(p, jumpsigma, f=explore, rng=rng) - if strictpriors(ϕ, record_max_age, climate_limits) + if strictpriors(ϕ, record_max_age, climate_limits, k) porewaterhistory!(sc, ϕ, k, climate, seawater, ka_dt) @@ -127,7 +130,7 @@ function porewatermetropolis(p::Proposal, jumpsigma::Proposal, prior::CoreData; @inbounds for i = 1:chainsteps ϕ, jumpname, jump = proposaljump(p, jumpsigma, f=explore, rng=rng) - if strictpriors(ϕ, record_max_age, climate_limits) + if strictpriors(ϕ, record_max_age, climate_limits, k) porewaterhistory!(sc, ϕ, k, climate, seawater, ka_dt) diff --git a/test/metro.jl b/test/metro.jl index 13b5794..640e286 100644 --- a/test/metro.jl +++ b/test/metro.jl @@ -1,10 +1,13 @@ # strictpriors -spargs = (1000., (1.,4.)) -@test PorewaterDiffusion.strictpriors(Proposal(500,1,1,2,3), spargs...) -@test !PorewaterDiffusion.strictpriors(Proposal(1100,1,1,2,3), spargs...) -@test !PorewaterDiffusion.strictpriors(Proposal(500,1,1,0,3), spargs...) -@test !PorewaterDiffusion.strictpriors(Proposal(500,1,1,2,5), spargs...) +spargs = (1000., (1.,4.), Constants()) +@test PorewaterDiffusion.strictpriors(Proposal(500,.1,1,2,3), spargs...) # passing test +@test !PorewaterDiffusion.strictpriors(Proposal(1100,.1,1,2,3), spargs...) # onset too old +@test !PorewaterDiffusion.strictpriors(Proposal(0,.1,1,2,3), spargs...) # onset too young +@test !PorewaterDiffusion.strictpriors(Proposal(500,.1,1,0,3), spargs...) # sea2frz too low +@test !PorewaterDiffusion.strictpriors(Proposal(500,.1,1,2,5), spargs...) # frz2mlt too high +@test !PorewaterDiffusion.strictpriors(Proposal(500,.1,1,3,2), spargs...) # sez2frz < frz2mlt +@test !PorewaterDiffusion.strictpriors(Proposal(500,.3,1,2,3), spargs...) # p.dfrz too high. # proposaljump