add histogram plotting and mean/med calc

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"

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)"
    ]
   }

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

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)
 

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

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

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