add basalCl and basalO to Proposal etc

src/history.jl

@@ -1,18 +1,18 @@
 """
 
 ```julia
-porewaterhistory!(sc::SedimentColumn, p::Proposal, k::Constants, climhist::NamedTuple, seawater::Water, basalwater::Water ka_dt::Int)
+porewaterhistory!(sc::SedimentColumn, p::Proposal, k::Constants, climhist::NamedTuple, seawater::Water, ka_dt::Int)
 ```
 
 In-place version of [`porewaterhistory`](@ref), which takes every input as an arg (rather than some defaults as kwargs). It also requires you to provide `ka_dt` -- the number of diffusion timesteps in each thousand-year climate timestep.
 
 see also: [`porewaterhistory`](@ref)
 
 """
-function porewaterhistory!(sc::SedimentColumn, p::Proposal, k::Constants, climhist::ClimateHistory, sw::Water, bw::Water, ka_dt::Int)
+function porewaterhistory!(sc::SedimentColumn, p::Proposal, k::Constants, climhist::ClimateHistory, sw::Water, ka_dt::Int)
     #sc.Cl.p .= sc.Cl.o .= sw.Cl
     #sc.O.p .= sc.O.o .= sw.O
-    equilibratecolumn!(sc,sw,bw,k.z,k.depth)
+    equilibratecolumn!(sc,sw,water(p.basalCl,p.basalO),k.z,k.depth)
 
     isd = searchsortedfirst(climhist.t, p.onset, rev=true)
     #isd = ifelse(isd<climhist.n, isd, climhist.n)
@@ -37,7 +37,7 @@ end
 """
 
 ```julia
-porewaterhistory(proposals [; k=Constants(), climatehistory=LR04(), seawater=mcmurdosound(), basalwater=deepbonney()])
+porewaterhistory(proposals [; k=Constants(), climatehistory=LR04(), seawater=mcmurdosound()])
 ```
 
 Calculate the porewater advection-diffusion history of chlorinity and O-isotope-traced water in a sediment column described by properties in `k` (::[`Constants`](@ref)) over a given [`ClimateHistory`](@ref) ([`LR04`](@ref) by default) and coretop `seawater` compositions.
@@ -49,10 +49,10 @@ See [`diffuseadvectcolumn!`](@ref) for the underlying diffusion-advection transp
 see also: [`porewaterhistory!`](@ref), [`Proposal`](@ref), [`Constants`](@ref), [`LR04`](@ref), [`water`](@ref)
 
 """
-function porewaterhistory(p::Proposal; k::Constants=Constants(), climatehistory::ClimateHistory=LR04(), seawater::Water=mcmurdosound(), basalwater::Water=deepbonney())
+function porewaterhistory(p::Proposal; k::Constants=Constants(), climatehistory::ClimateHistory=LR04(), seawater::Water=mcmurdosound())
 
     sc = SedimentColumn(k.nz,seawater...)
-    porewaterhistory!(sc, p, k, climatehistory, seawater, basalwater, dt_climatetimestep(climatehistory.t,k.dt))
+    porewaterhistory!(sc, p, k, climatehistory, seawater, dt_climatetimestep(climatehistory.t,k.dt))
 
     (; Cl = sc.Cl.p, d18O = sc.O.p)
 end

src/metropolis.jl

@@ -65,7 +65,7 @@ porewatermetropolis...
 ```
 Not tested, yet...
 """
-function porewatermetropolis(p::Proposal, jumpsigma::Proposal, prior::CoreData; burnin::Int=0, chainsteps::Int=100, k::Constants=Constants(), seawater::Water=mcmurdosound(), basalwater::Water=deepbonney(), explore::Tuple{Vararg{Symbol}}=fieldnames(Proposal), climate::ClimateHistory=LR04(), rng::AbstractRNG=Random.Xoshiro())
+function porewatermetropolis(p::Proposal, jumpsigma::Proposal, prior::CoreData; burnin::Int=0, chainsteps::Int=100, k::Constants=Constants(), seawater::Water=mcmurdosound(), explore::Tuple{Vararg{Symbol}}=fieldnames(Proposal), climate::ClimateHistory=LR04(), rng::AbstractRNG=Random.Xoshiro())
 
     scalejump=2.4
 
@@ -81,7 +81,7 @@ function porewatermetropolis(p::Proposal, jumpsigma::Proposal, prior::CoreData;
     sc = SedimentColumn(k.nz,seawater...)
     ka_dt = PorewaterDiffusion.dt_climatetimestep(climate.t,k.dt)
 
-    porewaterhistory!(sc, ϕ, k, climate, seawater, basalwater, ka_dt)
+    porewaterhistory!(sc, ϕ, k, climate, seawater, ka_dt)
 
     llCl, llO = loglikelihood(prior.z,prior.Cl.mu,prior.Cl.sig,k.z,sc.Cl.p), loglikelihood(prior.z,prior.O.mu,prior.O.sig,k.z,sc.O.p)
     ll = llCl + llO
@@ -99,7 +99,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, k)
 
-            porewaterhistory!(sc, ϕ, k, climate, seawater, basalwater, ka_dt)
+            porewaterhistory!(sc, ϕ, k, climate, seawater, ka_dt)
 
             llCl, llO = loglikelihood(prior.z,prior.Cl.mu,prior.Cl.sig,k.z,sc.Cl.p), loglikelihood(prior.z,prior.O.mu,prior.O.sig,k.z,sc.O.p)
             llϕ = llCl + llO
@@ -129,7 +129,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, k)
 
-            porewaterhistory!(sc, ϕ, k, climate, seawater, basalwater, ka_dt)
+            porewaterhistory!(sc, ϕ, k, climate, seawater, ka_dt)
 
             llCl, llO = loglikelihood(prior.z,prior.Cl.mu,prior.Cl.sig,k.z,sc.Cl.p), loglikelihood(prior.z,prior.O.mu,prior.O.sig,k.z,sc.O.p)
             llϕ = llCl + llO

src/parameters.jl

@@ -334,6 +334,8 @@ struct Proposal
     dmlt::Float64
     sea2frz::Float64
     frz2mlt::Float64
+    basalCl::Float64
+    basalO::Float64
 end
 
 
@@ -348,13 +350,15 @@ function update(x::Proposal, f::Symbol,v::Number)
     @assert f ∈ fieldnames(Proposal)
     v = float(v)
 
-    o = ifelse(f==:onset, v, x.onset)
-    df = ifelse(f==:dfrz, v, x.dfrz)
-    dm = ifelse(f==:dmlt, v, x.dmlt)
-    s = ifelse(f==:sea2frz, v, x.sea2frz)
-    f = ifelse(f==:frz2mlt, v, x.frz2mlt)
-
-    Proposal(o,df,dm,s,f)
+    Proposal(
+            ifelse(f==:onset, v, x.onset),
+            ifelse(f==:dfrz, v, x.dfrz),
+            ifelse(f==:dmlt, v, x.dmlt),
+            ifelse(f==:sea2frz, v, x.sea2frz),
+            ifelse(f==:frz2mlt, v, x.frz2mlt),
+            ifelse(f==:basalCl, v, x.basalCl),
+            ifelse(f==:basalO, v, x.basalO)
+    )
 end
 
 
@@ -376,5 +380,7 @@ function getproposal(x::Proposal, f::Symbol)
     y = ifelse(f==:dmlt, x.dmlt, y)
     y = ifelse(f==:sea2frz, x.sea2frz, y)
     y = ifelse(f==:frz2mlt, x.frz2mlt, y)
+    y = ifelse(f==:basalCl, x.basalCl, y)
+    y = ifelse(f==:basalO, x.basalO, y)
     y
 end

test/hist.jl

@@ -1,16 +1,17 @@
 @test PorewaterDiffusion.dt_climatetimestep(100:-1:0,10) == 100
 
-ch, k, sw, p = LR04(), Constants(), mcmurdosound(), Proposal(5320., 4e-5,1e-2,3.5, 4.2)
+ch, k, sw = LR04(), Constants(), mcmurdosound()
+p = Proposal(5320., 4e-5,1e-2,3.5, 4.2, sw...)
 sc = SedimentColumn(k.nz,sw...)
-porewaterhistory!(sc,p, k, ch,sw, sw, PorewaterDiffusion.dt_climatetimestep(ch.t,k.dt))
+porewaterhistory!(sc,p, k, ch,sw, PorewaterDiffusion.dt_climatetimestep(ch.t,k.dt))
 
 @test 48 < sc.Cl.p[end-1] < 52
 @test 19 < sc.Cl.p[2] < 20
 
 @test -2.0 < sc.O.p[end-1] < -1.9
 @test -2.6 < sc.O.p[2] < -2.5
 
-pwh2 = porewaterhistory(p,k=k,climatehistory=ch,seawater=sw,basalwater=sw)
+pwh2 = porewaterhistory(p,k=k,climatehistory=ch,seawater=sw)
 
 @test 48 < pwh2.Cl[end-1] < 52
 @test 19 < pwh2.Cl[2] < 20

test/metro.jl

@@ -1,20 +1,20 @@
 
 # strictpriors
 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.
-@test !PorewaterDiffusion.strictpriors(Proposal(500,.1,20.,2,3), spargs...) # p.dfrz too high.
+@test PorewaterDiffusion.strictpriors(Proposal(500,.1,1,2,3,1,1), spargs...) # passing test
+@test !PorewaterDiffusion.strictpriors(Proposal(1100,.1,1,2,3,1,1), spargs...) # onset too old
+@test !PorewaterDiffusion.strictpriors(Proposal(0,.1,1,2,3,1,1), spargs...) # onset too young
+@test !PorewaterDiffusion.strictpriors(Proposal(500,.1,1,0,3,1,1), spargs...) # sea2frz too low
+@test !PorewaterDiffusion.strictpriors(Proposal(500,.1,1,2,5,1,1), spargs...) # frz2mlt too high
+@test !PorewaterDiffusion.strictpriors(Proposal(500,.1,1,3,2,1,1), spargs...) # sez2frz < frz2mlt 
+@test !PorewaterDiffusion.strictpriors(Proposal(500,.3,1,2,3,1,1), spargs...) # p.dfrz too high.
+@test !PorewaterDiffusion.strictpriors(Proposal(500,.1,20.,2,3,1,1), spargs...) # p.dfrz too high.
 
 
 # proposaljump
-@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(7)...), Proposal(ones(7)...), 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
+@test PorewaterDiffusion.proposaljump(Proposal(ones(7)...), Proposal(ones(7)...), rng=StableRNG(10))[1].onset ≈ 0.23064139237391934
 
 # stopwatch
 @test "0% |■■□□□□□□□□| 100%  ||  step: 27 / 100  ||  time: 0.0 m" == PorewaterDiffusion.stopwatch(27,100,time())

test/param.jl

@@ -35,10 +35,11 @@ ktest = Constants(k=0.1, dz=5, dt=10, depth=2000)
 @test last(ktest.k2w) ≈ 0.0030695870568939743
 
 
-proposaltest = Proposal(1,1,1,1,1)
-@test update(proposaltest, :onset, 2.) == Proposal(2,1,1,1,1)
-@test update(proposaltest, :dfrz, 2.) == Proposal(1,2,1,1,1)
-@test update(proposaltest, :dmlt, 2.) == Proposal(1,1,2,1,1)
-@test update(proposaltest, :sea2frz, 2.) == Proposal(1,1,1,2,1)
-@test update(proposaltest, :frz2mlt, 2.) == Proposal(1,1,1,1,2)
-
+proposaltest = Proposal(1,1,1,1,1,1,1)
+@test update(proposaltest, :onset, 2.) == Proposal(2,1,1,1,1,1,1)
+@test update(proposaltest, :dfrz, 2.) == Proposal(1,2,1,1,1,1,1)
+@test update(proposaltest, :dmlt, 2.) == Proposal(1,1,2,1,1,1,1)
+@test update(proposaltest, :sea2frz, 2.) == Proposal(1,1,1,2,1,1,1)
+@test update(proposaltest, :frz2mlt, 2.) == Proposal(1,1,1,1,2,1,1)
+@test update(proposaltest, :basalCl, 2.) == Proposal(1,1,1,1,1,2,1)
+@test update(proposaltest, :basalO, 2.) == Proposal(1,1,1,1,1,1,2)