Skip to content

Commit

Permalink
make dfrz/dmlt normal jumps again (not lognormal)
Browse files Browse the repository at this point in the history
  • Loading branch information
grahamedwards committed Mar 1, 2024
1 parent 692377d commit 36521c4
Show file tree
Hide file tree
Showing 2 changed files with 12 additions and 18 deletions.
28 changes: 11 additions & 17 deletions src/metropolis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,22 +19,18 @@ end
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)).
- Annual melting rate must be no more than that observed at the Thwaites grounding line (<10 m/yr, [Davis+ 2023](https://www.nature.com/articles/s41586-022-05586-0)).
- Freezing rate is non-zero and ≤ 0.4dt/dz (to prevent an error in a log-calculation in [`PorewaterDiffusion.boundaryconditions`](@ref)).
- Annual melting rate is non-zero and must be no more than that observed at the Thwaites grounding line (<10 m/yr, [Davis+ 2023](https://www.nature.com/articles/s41586-022-05586-0)).
"""
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=true

x &= p.sea2frz < p.frz2mlt
x &= climatelimits[1] < p.sea2frz
x &= p.frz2mlt < climatelimits[2]

x &= p.dfrz <= 0.4k.dz/k.dt
x &= p.dmlt <= 10.
x &= 0 < p.onset <= record_max_age
x &= climatelimits[1] < p.sea2frz < p.frz2mlt < climatelimits[2]
x &= 0 < p.dfrz <= 0.4k.dz/k.dt
x &= 0 < p.dmlt <= 10.

x
end
Expand All @@ -47,18 +43,16 @@ end
Add a random jump to a randomly selected field of `p` with a corresponding normal jumping distribution defined by the corresponding field in `σ`. The possible fields may be specified by providing a Tuple of Symbols `f`, and a specific RNG seed may be provided.
Note that `:dmlt` and `:dfrz` are drawn from a lognormal jumping distribution (where `σ` is in log-space.)
"""
"""#Note that `:dmlt` and `:dfrz` are drawn from a lognormal jumping distribution (where `σ` is in log-space.)
function proposaljump(p::Proposal, j::Proposal; rng::AbstractRNG=Xoshiro(), f::Tuple{Vararg{Symbol}}=fieldnames(Proposal))

jumpname = rand(rng,f)
logdist = (jumpname == :dmlt) | (jumpname == :dfrz)
#logdist = (jumpname == :dmlt) | (jumpname == :dfrz)
jump = getproposal(j,jumpname) * randn(rng)
x = getproposal(p,jumpname)
x = ifelse(logdist, log(x),x)
#x = ifelse(logdist, log(x),x)
x += jump
x = ifelse(logdist, exp(x),x)
#x = ifelse(logdist, exp(x),x)
(update(p, jumpname, x) , jumpname , abs(jump) )
end

Expand Down
2 changes: 1 addition & 1 deletion test/metro.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ spargs = (1000., (1.,4.), Constants())


# proposaljump
@test PorewaterDiffusion.proposaljump(Proposal(ones(5)...), Proposal(ones(5)...), rng=StableRNG(2580))[1].dfrz 1.0949056480096304
@test PorewaterDiffusion.proposaljump(Proposal(ones(5)...), Proposal(ones(5)...), rng=StableRNG(2580))[1].dfrz 1.090668193354693 # if lognormal -> 1.0949056480096304

@test PorewaterDiffusion.proposaljump(Proposal(ones(5)...), Proposal(ones(5)...), rng=StableRNG(10))[1].onset 0.23064139237391934

Expand Down

0 comments on commit 36521c4

Please sign in to comment.