diff --git a/Project.toml b/Project.toml index 3212a19..2bd9b86 100644 --- a/Project.toml +++ b/Project.toml @@ -4,8 +4,11 @@ authors = ["Graham Harper Edwards", "Sarah U. Neuhaus"] version = "1.1.0-DEV" [deps] +CleanHistograms = "77dce33c-ab31-460d-b478-44a825891bef" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Requires = "ae029012-a4dd-5104-9daa-d747884805df" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [compat] DelimitedFiles = "1" diff --git a/example.ipynb b/example.ipynb index e444435..26eca26 100644 --- a/example.ipynb +++ b/example.ipynb @@ -26,15 +26,8 @@ ], "source": [ "p = Proposal(5320., 1e-4,1e-3,3.5, 4.2, 1000, deepbonney()...)\n", - "jumpsize = Proposal(20., .1e-4, .1e-3, .1, .1,10, 10, 1.)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ + "jumpsize = Proposal(20., .1e-4, .1e-3, .1, .1,10, 10, 1.)\n", + "\n", "chains = porewatermetropolis(p, jumpsize, andrill2a(); burnin=100, chainsteps=100, k=Constants(), seawater=mcmurdosound(), climate=LR04(), onlychloride=true)" ] } diff --git a/src/CorePore.jl b/src/CorePore.jl index 31554d7..eb53a05 100644 --- a/src/CorePore.jl +++ b/src/CorePore.jl @@ -1,5 +1,5 @@ module CorePore -import DelimitedFiles +import DelimitedFiles, Requires, Statistics using Random export Water, mcmurdoshelf, mcmurdosound, deepbonney, CoreData, andrill2a, andrill1b, PorewaterProperty, SedimentColumn, LR04, Constants, Proposal, update @@ -11,10 +11,15 @@ include("diffuse.jl") export porewaterhistory, porewaterhistory! include("history.jl") -export loglikelihood, normpdf +export loglikelihood, normpdf, means, medians include("statistics.jl") export porewatermetropolis include("metropolis.jl") +function __init__() + Requires.@require CairoMakie="13f3f980-e62b-5c42-98c6-ff1f3baf88f0" include("plot.jl") + Requires.@require GLMakie="e9467ef8-e4e7-5192-8a1a-b1aee30e663a" include("plot.jl") + Requires.@require WGLMakie="276b4fcb-3e11-5398-bf8b-a0c2d153d008" include("plot.jl") +end end diff --git a/src/parameters.jl b/src/parameters.jl index 585c31c..0a6ac97 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -334,7 +334,7 @@ Struct containing porewater parameters (described below). All inputs must be of | `frz2mlt` | Benthic δ¹⁸O threshold for subglacial melting | ‰ | 4.2 | | `flr` | depth of diffusion-dominated porewater column | m | 1000. | | `basalCl` | chloridity at base of diffusion-dominated column | g/kg | 119.4 | -| `basalO` | δ¹⁸O at base of diffusion-dominated column | g/kg | -25.2 | +| `basalO` | δ¹⁸O at base of diffusion-dominated column | ‰ | -25.2 | see also: [`update`](@ref) diff --git a/src/plot.jl b/src/plot.jl new file mode 100644 index 0000000..623d553 --- /dev/null +++ b/src/plot.jl @@ -0,0 +1,38 @@ +Requires.@require CairoMakie="13f3f980-e62b-5c42-98c6-ff1f3baf88f0" @eval using .CairoMakie +Requires.@require GLMakie="e9467ef8-e4e7-5192-8a1a-b1aee30e663a" @eval using .WGLMakie +Requires.@require WGLMakie="276b4fcb-3e11-5398-bf8b-a0c2d153d008" @eval using .WGLMakie + +import CleanHistograms +export histograms, defaultlabels + +defaultlabels() = (onset = "Onset date (ka)", dfrz = "Freezing rate (m/yr)", dmlt = "Melting rate (m/yr)", sea2frz = "Basal freezing (‰, δ¹⁸O)", frz2mlt = "Basal melting (‰, δ¹⁸O)", flr = "Column depth (m)", basalCl = "Deep chloridity (g/kg)", basalO = "Deep δ¹⁸O (‰)") + + +function histograms(x::NamedTuple; f=Makie.Figure(size=(600,600)), panels::NamedTuple = defaultlabels(), cols::Int=0, bins::Int=32, darkmode::Bool=false) + n = length(panels) + + ks = keys(panels) + + pltclr = ifelse(darkmode,:white,:black) + + cols = ifelse(cols < 1, ceil(Int,sqrt(n)), cols) + rows = ceil(Int,n/cols) + + @inbounds for i = 1:rows, j=1:cols + ij = cols * (i-1) + j + if ij < n + k = ks[ij] + @assert k ∈ keys(x) "Provided panel $k is not a key in chains provided in `x`" + + xk = x[k] + + ax = Makie.Axis(f[i,j], xlabel=panels[k], + bottomspinecolor=pltclr,xtickcolor=pltclr,xticklabelcolor=pltclr, xlabelcolor=pltclr,backgroundcolor=ifelse(darkmode,:transparent,:white), + xgridvisible=false,ygridvisible=false,yticklabelsvisible=false,yticksvisible=false,rightspinevisible=false,leftspinevisible=false,topspinevisible=false) + h = CleanHistograms.cleanhist(xk,bins=bins) + band!(ax,h.x,h.y,zero(h.y), color=(pltclr,0.1)) + lines!(ax,h.x,h.y, color=pltclr, linewidth=2,) + end + end + f +end diff --git a/src/statistics.jl b/src/statistics.jl index d0d55b7..ecf33dd 100644 --- a/src/statistics.jl +++ b/src/statistics.jl @@ -58,4 +58,42 @@ function loglikelihood(zo::T, muo::T, sigo::T, zm::AbstractRange{Float64}, m::T) ll += normll(muo[i], sigo[i], xmi) end ll +end + + +""" + + means(x; std::Integer) + +Calculate the means of each field in NamedTuple `x`. Optionally provide an integer number of standard deviations to calculate, returned as the form `key` = (m=μ, s=σ). Use with the returned chains of [`porewatermetropolis`](@ref). + +""" +function means(x::NamedTuple; std::Int=0) + k = keys(x) + if std > 0 + NamedTuple{k}( (m=Statistics.mean(x[i]), s=std*Statistics.std(x[i])) for i in k) + else + NamedTuple{k}(Statistics.mean(x[i]) for i in k) + end +end + + + + +""" + + medians(x; ci::Float64) + +Calculate the medians of each field in NamedTuple `x`. Optionally provide a `ci` ∈ [0,1] of standard deviations to calculate, returned as the form `key` = (m=μ, s=σ). Use with the returned chains of [`porewatermetropolis`](@ref). + +""" +function medians(x; ci::Float64=0.) + @assert 0 ≤ ci ≤ 1 + k = keys(x) + if ci>0. + upper, lower = (0.5 + 0.5ci), 0.5ci + NamedTuple{k}( (m=Statistics.median(x[i]), l=Statistics.quantile(x[i],lower), u=Statistics.quantile(x[i],upper)) for i in k) + else + NamedTuple{k}(Statistics.median(x[i]) for i in k) + end end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 35a4e3c..f27a2d5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,8 +12,6 @@ macro silence(block) end end -mean(x) = sum(x)/length(x) - @testset "parameters" begin include("param.jl") end @testset "diffusion" begin include("diff.jl") end