Refactor pb-loss-aware eruption age estimation, add options

.gitignore

@@ -1,4 +1,5 @@
 .DS_Store
 Manifest.toml
 /dev
-.vscode
+.vscode
+*.pdf

Project.toml

@@ -1,15 +1,19 @@
 name = "Isoplot"
 uuid = "5adc30d5-9ddf-423c-bb15-ece697bec3ab"
 authors = ["C. Brenhin Keller <[email protected]>"]
-version = "0.3.7"
+version = "0.3.8"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 NaNStatistics = "b946abbf-3ea7-4610-9019-9858bfdeaf2d"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [weakdeps]
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
@@ -21,9 +25,13 @@ PlotsExt = "Plots"
 
 [compat]
 Distributions = "0.25"
+LogExpFunctions = "0.3.29"
 LoopVectorization = "0.12"
 Measurements = "2"
 NaNStatistics = "0.6"
 Plots = "1"
+Reexport = "1.2.2"
+Rotations = "1.7"
 SpecialFunctions = "2"
+StaticArrays = "1.9"
 julia = "1.9"

examples/demo.jl

@@ -1,4 +1,4 @@
-using Isoplot, Plots, VectorizedStatistics
+using Isoplot, Plots
 
 cd(@__DIR__)
 
@@ -33,7 +33,7 @@ display(hdl)
 
 # Plot in rankorder plot
 age_06_38 = last.(age.(analyses))
-rankorder_plot = rankorder(val.(age_06_38), σ1.(age_06_38), ylabel="Age (Ma)")
+rankorder_plot = rankorder(Isoplot.val.(age_06_38), 2*Isoplot.err.(age_06_38), ylabel="Age (Ma)")
 savefig(rankorder_plot, "rank_order.pdf")
 display(rankorder_plot)
 
@@ -61,17 +61,16 @@ display(h)
 
 ## --- Show eruption age relative to distribution of upper intercepts
 
-uis = upperintercept.(vmean(t0dist), analyses)
-h = rankorder(Isoplot.val.(uis), Isoplot.err.(uis))
+uis = upperintercept.(nanmean(t0dist), analyses)
+h = rankorder(Isoplot.val.(uis), 2*Isoplot.err.(uis))
 plot!(h,1:length(uis),fill(terupt.lower,length(uis)),fillto=terupt.upper,color=:blue,fillalpha=0.5,linealpha=0, label="Model ($terupt Ma, 95% CI)")
 plot!(h,1:length(uis),fill(terupt.mean,length(uis)),linecolor=:black,linestyle=:dot,label="",legend=:topleft,fg_color_legend=:white,framestyle=:box)
 
-
-
 ## --- Histogram of distribution of time of Pb-loss
 
 h = histogram(t0dist, xlabel="Age [Ma]", ylabel="Probability Density", normalize=true, label="Time of Pb-loss", color=:darkblue, alpha=0.65, linealpha=0.1, framestyle=:box)
-xlims!(h, 0, last(xlims()))
-ylims!(h, 0, last(ylims()))
+plot!(h, xlims=(0,last(xlims())), ylims=(0, last(ylims())))
 savefig(h, "PbLoss.pdf")
 display(h)
+
+## --- End of File

examples/pbloss.jl

@@ -1,4 +1,4 @@
-using Isoplot, Plots, VectorizedStatistics
+using Isoplot, Plots
 
 # Example U-Pb dataset (MacLennan et al. 2020)
 #       207/235  1σ abs   206/236     1σ abs     correlation
@@ -84,4 +84,4 @@ display(teruptold)
 # h = plot(ccfiltered,cc, layout=(2,1), size=(500,660), left_margin=(4,:mm))
 # savefig(h, "filteredunfiltered.pdf")
 
-## ---
+## --- End of File

ext/PlotsExt.jl

@@ -104,7 +104,7 @@ module PlotsExt
         xoffset = (xl[2]-xl[1])/200
         yoffset = (yl[2]-yl[1])/100
         t = (xl[1]+8*xoffset) .< r75tticks .< xl[2]
-        ticklabels = Plots.text.(string.(round.(Int, tticks[t])), 10, :right)
+        ticklabels = Plots.text.(string.(round.(tticks[t], digits=1)), 10, :right)
         Plots.annotate!(hdl, r75tticks[t].-xoffset,r68tticks[t].+yoffset,ticklabels)
 
         return hdl

src/Isoplot.jl

@@ -1,10 +1,14 @@
 module Isoplot
 
-    using NaNStatistics
+    using Reexport
+    @reexport using NaNStatistics
+    using Distributions
     using LoopVectorization: @turbo
+    using LogExpFunctions
     using LinearAlgebra
-    using Distributions
     using Measurements
+    using StaticArrays
+    using Rotations
 
 
     # A type alias for array-ish types

src/concordia.jl

@@ -61,9 +61,9 @@ function upperintercept(tₗₗ::Number, d::UPbAnalysis{T}, nresamplings::Intege
     # Get ratios
     r75₀, r68₀ = d.μ
     # Return early if our lead loss time is too old or anything is NaN'd
-    tₗₗ < log(r68₀+1)/λ238U.val || return fill!(uis, T(NaN))
-    tₗₗ < log(r75₀+1)/λ235U.val || return fill!(uis, T(NaN))
-
+    tₗₗ < log(r68₀+1)/λ238U.val || return fill(T(NaN), nresamplings)
+    tₗₗ < log(r75₀+1)/λ235U.val || return fill(T(NaN), nresamplings)
+    
     # Calculate isotopic ratios of our time of Pb-loss
     r75ₗₗ = exp(λ235U.val*tₗₗ) - 1
     r68ₗₗ = exp(λ238U.val*tₗₗ) - 1