Skip to content

Commit

Permalink
Update examples
Browse files Browse the repository at this point in the history
  • Loading branch information
@brenhinkeller
brenhinkeller committed Apr 13, 2024
1 parent 88050ea commit 5a61cdc
Show file tree
Hide file tree
Showing 3 changed files with 128 additions and 77 deletions.
88 changes: 70 additions & 18 deletions examples/ZrnHeInversionVartCryst.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,22 +35,77 @@
ds = importdataset("minnesota.csv", ',', importas=:Tuple)

# Populate data NamedTuple from imported dataset
data = (
halfwidth = ds.Halfwidth_um, # Crystal half-width, in microns
U = ds.U238_ppm, # U concentration, in PPM
Th = ds.Th232_ppm, # Th-232 concentration, in PPM
Sm = ds.Sm147_ppm, # Sm-147 concentration, in PPM (optional)
HeAge = ds.HeAge_Ma_raw, # He age, in Ma
HeAge_sigma = ds.HeAge_Ma_sigma_10pct, # He age uncertainty (1-sigma), in Ma
crystAge = ds.CrystAge_Ma, # Crystallization age, in Ma
mineral = ds.Mineral
data = (;
halfwidth = copy(ds.Halfwidth_um), # Crystal half-width, in microns
U = copy(ds.U238_ppm), # U concentration, in PPM
Th = copy(ds.Th232_ppm), # Th-232 concentration, in PPM
Sm = copy(ds.Sm147_ppm), # Sm-147 concentration, in PPM (optional)
HeAge = copy(ds.HeAge_Ma_raw), # He age, in Ma
HeAge_sigma = copy(ds.HeAge_Ma_sigma_raw), # He age uncertainty (1-sigma), in Ma
crystAge = copy(ds.CrystAge_Ma), # Crystallization age, in Ma
mineral = copy(ds.Mineral), # i.e., "zircon" or "apatite"
);

ta = containsi.(data.mineral, "apatite")
tz = containsi.(data.mineral, "zircon")
eU = data.U+0.238*data.Th # Used for plotting, etc.

## --- Empirical uncertainty estimation

age =copy(data.HeAge)
age_sigma = copy(ds.HeAge_Ma_sigma_raw)
age_sigma_empirical = copy(data.HeAge_sigma)
min_rel_uncert = 10/100 # 10% minmum relative age uncertainty

if any(ta)
# Standard deviation of a Gaussian kernel in eU space, representing the
# range of eU over which zircons with similar eU should have similar ages
σeU = 10

# Calculate errors
for i ∈ findall(ta)
W = normpdf.(eU[i], σeU, eU[ta])
σ_external = nanstd(age[ta], W) # Weighted standard deviation
σ_internal = age_sigma[i]
age_sigma_empirical[i] = sqrt(σ_external^2 + σ_internal^2)
age_sigma_empirical[i] = max(age_sigma_empirical[i], min_rel_uncert*age[i])
end

h = plot(xlabel="eU", ylabel="Age", framestyle=:box, title="apatite")
plot!(eU[ta], age[ta], yerror=age_sigma_empirical[ta], seriestype=:scatter, c=:black, msc=:black, label="empirical")
plot!(eU[ta], age[ta], yerror=age_sigma[ta], seriestype=:scatter, c=mineralcolors["apatite"], msc=mineralcolors["apatite"], label="internal")
display(h)
end

if any(tz)
# Stzndard deviation of a Gaussian kernel in eU space, representing the
# range of eU over which zircons with similar eU should have similar ages
σeU = 100

# Calculate errors
for i ∈ findall(tz)
W = normpdf.(eU[i], σeU, eU[tz])
σ_external = nanstd(age[tz], W) # Weighted stzndard deviation
σ_internal = age_sigma[i]
age_sigma_empirical[i] = sqrt(σ_external^2 + σ_internal^2)
age_sigma_empirical[i] = max(age_sigma_empirical[i], min_rel_uncert*age[i])
end

h = plot(xlabel="eU", ylabel="Age", framestyle=:box, title="zircon")
plot!(eU[tz], age[tz], yerror=age_sigma_empirical[tz], seriestype=:scatter, c=:black, msc=:black, label="empirical")
plot!(eU[tz], age[tz], yerror=age_sigma[tz], seriestype=:scatter, c=mineralcolors["zircon"], msc=mineralcolors["zircon"], label="internal")
display(h)
end

## --- Prepare problem

# Use empirical ages for zircon, 10 % for apatite
data.HeAge_sigma[tz] .= age_sigma_empirical[tz]
data.HeAge_sigma[ta] .= ds.HeAge_Ma_sigma_10pct[ta]

model = (
nsteps = 20_000, # How many steps of the Markov chain should we run?
burnin = 10_000, # How long should we wait for MC to converge (become stationary)
nsteps = 100_000, # How many steps of the Markov chain should we run?
burnin = 50_000, # How long should we wait for MC to converge (become stationary)
dr = 1.0, # Radius step, in microns
dt = 10.0, # Time step size in Myr
dTmax = 25.0, # Maximum reheating/burial per model timestep. If too high, may cause numerical problems in Crank-Nicholson solve
Expand Down Expand Up @@ -174,8 +229,6 @@

## --- Plot model ages vs observed ages in age-eU space (zircon)

eU = data.U+.238*data.Th # Used only for plotting
tz = containsi.(data.mineral, "zircon")
if any(tz)
h = scatter(eU[tz], data.HeAge[tz],
yerror=2*data.HeAge_sigma[tz],
Expand All @@ -192,16 +245,15 @@
scatter!(h, eU[tz], m,
yerror=(m-l, u-m),
label="Model + 95%CI",
color=:blue,
msc=:blue,
color=mineralcolors["zircon"],
msc=mineralcolors["zircon"],
)
savefig(h, name*"_zircon_Age-eU.pdf")
display(h)
end

## --- Plot model ages vs observed ages in age-eU space (apatite)

ta = containsi.(data.mineral, "apatite")
if any(ta)
h = scatter(eU[ta], data.HeAge[ta],
yerror=2*data.HeAge_sigma[ta],
Expand All @@ -218,8 +270,8 @@
scatter!(h, eU[ta], m,
yerror=(m-l, u-m),
label="Model + 95%CI",
color=:blue,
msc=:blue,
color=mineralcolors["apatite"],
msc=mineralcolors["apatite"],
)
savefig(h, name*"_apatite_Age-eU.pdf")
display(h)
Expand Down
69 changes: 34 additions & 35 deletions examples/minnesota.csv
View file
Original file line number Diff line number Diff line change
@@ -1,35 +1,34 @@
Grain_Name,Mineral,CrystAge_Ma,Halfwidth_um,HeAge_Ma_raw,HeAge_Ma_sigma_raw,HeAge_Ma_sigma_5pct,HeAge_Ma_sigma_10pct,HeAge_Ma_sigma_emp,U238_ppm,Th232_ppm,Sm147_ppm,HeAge_reference
04RF1zA,zircon,3500,64,770,15,38.5,77,100.4,84,67,0,Guenthner et al. 2013
04SC1zB,zircon,1779,71,659,51,32.95,65.9,114.6,415,111,0,Guenthner et al. 2013
04SC1zA,zircon,1779,47,649,13,32.45,64.9,107,422,110,0,Guenthner et al. 2013
04GF1zB,zircon,3500,63,638,12.6,31.9,63.8,69.3,254,104,0,Guenthner et al. 2013
04GF1zA,zircon,3500,56,620,18.4,31,62,71.3,247,88,0,Guenthner et al. 2013
04RF1zB,zircon,3500,72,557,10,27.85,55.7,82.7,155,177,0,Guenthner et al. 2013
04MT1zA,zircon,3500,67,545,12,27.25,54.5,80.2,362,89,0,Guenthner et al. 2013
04R1zB,zircon,1783,51,500,10,25,50,137.8,580,174,0,Guenthner et al. 2013
04MT1zB,zircon,3500,68,493,43,24.65,49.3,86.2,346,78,0,Guenthner et al. 2013
04R1zA,zircon,1783,42,357,7,17.85,35.7,116.4,656,167,0,Guenthner et al. 2013
04MT1zF,zircon,3500,54,329,7,16.45,32.9,147.7,567,75,0,Miltich 2005
04EQ1zA,zircon,3500,41,253,5.2,12.65,25.3,101.7,781,203,0,Guenthner et al. 2013
04EQ1zF,zircon,3500,43,241,4.6,12.05,24.1,112.4,647,284,0,Guenthner et al. 2013
04SH1zB,zircon,2603,56,225,4.6,11.25,22.5,95.9,883,184,0,Guenthner et al. 2013
04EQ1zB,zircon,3500,42,225,6,11.25,22.5,102.8,748,218,0,Guenthner et al. 2013
04MT1zD,zircon,3500,80,217,5,10.85,21.7,104.4,739,141,0,Miltich 2005
04MT1zC,zircon,3500,45,193,4,9.65,19.3,88.9,980,252,0,Miltich 2005
04MT1zE,zircon,3500,49,190,4,9.5,19,89.5,974,197,0,Miltich 2005
04SG1zB,zircon,3500,55,72,2,3.6,36,36,1738,1171,0,Guenthner et al. 2013
04EQ1zC,zircon,3500,36,57,1.1,2.85,28.5,93.6,894,246,0,Guenthner et al. 2013
04EQ1zE,zircon,3500,35,42,1.1,2.1,21,74.8,1107,351,0,Guenthner et al. 2013
04EQ1zD,zircon,3500,31,29,0.8,1.45,14.5,97.6,866,166,0,Guenthner et al. 2013
04SH1zA,zircon,2603,50,11,0.2,0.55,5.5,88.1,1011,228,0,Guenthner et al. 2013
04GF1a,apatite,3500,29,234,9.4,11.7,23.4,,33,24.4,746.39,Miltich 2005
04GF1AB,apatite,3500,57,98,2,4.9,9.8,,6.6,1.8,0,Miltich 2005
04GF1AC,apatite,3500,59,233,5,11.65,23.3,,10.3,6.3,0,Miltich 2005
04SG1AB,apatite,3500,51,339,6,16.95,33.9,,107.1,91,0,Miltich 2005
04SH1AB,apatite,2604,43,378,9.3,18.9,37.8,,19.4,25.6,155.41,Miltich 2005
04SH1AC,apatite,2604,52,258,5,12.9,25.8,,16.9,30.2,0,Miltich 2005
04MT1a,apatite,3500,44,158,3.3,7.9,15.8,,53.3,36.7,240.99,Miltich 2005
04MT1AB,apatite,3500,53,269,8.1,13.45,26.9,,11,5.2,144.26,Miltich 2005
04MT1AC,apatite,3500,79,313,5.8,15.65,31.3,,25.9,16.8,111.11,Miltich 2005
04EQ1AB,apatite,3500,76,309,6,15.45,30.9,,21.9,8.7,0,Miltich 2005
04EQ1AC,apatite,3500,95,392,8,19.6,39.2,,46,6.4,0,Miltich 2005
Grain_Name,Mineral,CrystAge_Ma,Halfwidth_um,HeAge_Ma_raw,HeAge_Ma_sigma_raw,HeAge_Ma_sigma_5pct,HeAge_Ma_sigma_10pct,U238_ppm,Th232_ppm,Sm147_ppm,HeAge_reference
04RF1zA,zircon,3500,64,770,15,38.5,77,84,67,0,Guenthner et al. 2013
04SC1zB,zircon,1779,71,659,51,32.95,65.9,415,111,0,Guenthner et al. 2013
04SC1zA,zircon,1779,47,649,13,32.45,64.9,422,110,0,Guenthner et al. 2013
04GF1zB,zircon,3500,63,638,12.6,31.9,63.8,254,104,0,Guenthner et al. 2013
04GF1zA,zircon,3500,56,620,18.4,31,62,247,88,0,Guenthner et al. 2013
04RF1zB,zircon,3500,72,557,10,27.85,55.7,155,177,0,Guenthner et al. 2013
04MT1zA,zircon,3500,67,545,12,27.25,54.5,362,89,0,Guenthner et al. 2013
04R1zB,zircon,1783,51,500,10,25,50,580,174,0,Guenthner et al. 2013
04MT1zB,zircon,3500,68,493,43,24.65,49.3,346,78,0,Guenthner et al. 2013
04R1zA,zircon,1783,42,357,7,17.85,35.7,656,167,0,Guenthner et al. 2013
04MT1zF,zircon,3500,54,329,7,16.45,32.9,567,75,0,Miltich 2005
04EQ1zA,zircon,3500,41,253,5.2,12.65,25.3,781,203,0,Guenthner et al. 2013
04EQ1zF,zircon,3500,43,241,4.6,12.05,24.1,647,284,0,Guenthner et al. 2013
04SH1zB,zircon,2603,56,225,4.6,11.25,22.5,883,184,0,Guenthner et al. 2013
04EQ1zB,zircon,3500,42,225,6,11.25,22.5,748,218,0,Guenthner et al. 2013
04MT1zD,zircon,3500,80,217,5,10.85,21.7,739,141,0,Miltich 2005
04MT1zC,zircon,3500,45,193,4,9.65,19.3,980,252,0,Miltich 2005
04MT1zE,zircon,3500,49,190,4,9.5,19,974,197,0,Miltich 2005
04EQ1zC,zircon,3500,36,57,1.1,2.85,28.5,894,246,0,Guenthner et al. 2013
04EQ1zE,zircon,3500,35,42,1.1,2.1,21,1107,351,0,Guenthner et al. 2013
04EQ1zD,zircon,3500,31,29,0.8,1.45,14.5,866,166,0,Guenthner et al. 2013
04SH1zA,zircon,2603,50,11,0.2,0.55,5.5,1011,228,0,Guenthner et al. 2013
04GF1a,apatite,3500,29,234,9.4,11.7,23.4,33,24.4,746.39,Miltich 2005
04GF1AB,apatite,3500,57,98,2,4.9,9.8,6.6,1.8,0,Miltich 2005
04GF1AC,apatite,3500,59,233,5,11.65,23.3,10.3,6.3,0,Miltich 2005
04SG1AB,apatite,3500,51,339,6,16.95,33.9,107.1,91,0,Miltich 2005
04SH1AB,apatite,2604,43,378,9.3,18.9,37.8,19.4,25.6,155.41,Miltich 2005
04SH1AC,apatite,2604,52,258,5,12.9,25.8,16.9,30.2,0,Miltich 2005
04MT1a,apatite,3500,44,158,3.3,7.9,15.8,53.3,36.7,240.99,Miltich 2005
04MT1AB,apatite,3500,53,269,8.1,13.45,26.9,11,5.2,144.26,Miltich 2005
04MT1AC,apatite,3500,79,313,5.8,15.65,31.3,25.9,16.8,111.11,Miltich 2005
04EQ1AB,apatite,3500,76,309,6,15.45,30.9,21.9,8.7,0,Miltich 2005
04EQ1AC,apatite,3500,95,392,8,19.6,39.2,46,6.4,0,Miltich 2005
48 changes: 24 additions & 24 deletions examples/minnesotazrn.csv
View file
Original file line number Diff line number Diff line change
@@ -1,24 +1,24 @@
Grain_Name,Mineral,CrystAge_Ma,Halfwidth_um,HeAge_Ma_raw,HeAge_Ma_sigma_raw,HeAge_Ma_sigma_5pct,HeAge_Ma_sigma_10pct,HeAge_Ma_sigma_emp,U238_ppm,Th232_ppm,HeAge_reference
04RF1zA,zircon,3500,64,770,15,38.5,77,100.4,84,67,Guenthner et al. 2013
04SC1zB,zircon,1779,71,659,51,32.95,65.9,114.6,415,111,Guenthner et al. 2013
04SC1zA,zircon,1779,47,649,13,32.45,64.9,107,422,110,Guenthner et al. 2013
04GF1zB,zircon,3500,63,638,12.6,31.9,63.8,69.3,254,104,Guenthner et al. 2013
04GF1zA,zircon,3500,56,620,18.4,31,62,71.3,247,88,Guenthner et al. 2013
04RF1zB,zircon,3500,72,557,10,27.85,55.7,82.7,155,177,Guenthner et al. 2013
04MT1zA,zircon,3500,67,545,12,27.25,54.5,80.2,362,89,Guenthner et al. 2013
04R1zB,zircon,1783,51,500,10,25,50,137.8,580,174,Guenthner et al. 2013
04MT1zB,zircon,3500,68,493,43,24.65,49.3,86.2,346,78,Guenthner et al. 2013
04R1zA,zircon,1783,42,357,7,17.85,35.7,116.4,656,167,Guenthner et al. 2013
04MT1zF,zircon,3500,54,329,7,16.45,32.9,147.7,567,75,Miltich 2005
04EQ1zA,zircon,3500,41,253,5.2,12.65,25.3,101.7,781,203,Guenthner et al. 2013
04EQ1zF,zircon,3500,43,241,4.6,12.05,24.1,112.4,647,284,Guenthner et al. 2013
04SH1zB,zircon,2603,56,225,4.6,11.25,22.5,95.9,883,184,Guenthner et al. 2013
04EQ1zB,zircon,3500,42,225,6,11.25,22.5,102.8,748,218,Guenthner et al. 2013
04MT1zD,zircon,3500,80,217,5,10.85,21.7,104.4,739,141,Miltich 2005
04MT1zC,zircon,3500,45,193,4,9.65,19.3,88.9,980,252,Miltich 2005
04MT1zE,zircon,3500,49,190,4,9.5,19,89.5,974,197,Miltich 2005
04SG1zB,zircon,3500,55,72,2,3.6,36,36,1738,1171,Guenthner et al. 2013
04EQ1zC,zircon,3500,36,57,1.1,2.85,28.5,93.6,894,246,Guenthner et al. 2013
04EQ1zE,zircon,3500,35,42,1.1,2.1,21,74.8,1107,351,Guenthner et al. 2013
04EQ1zD,zircon,3500,31,29,0.8,1.45,14.5,97.6,866,166,Guenthner et al. 2013
04SH1zA,zircon,2603,50,11,0.2,0.55,5.5,88.1,1011,228,Guenthner et al. 2013
Grain_Name,Mineral,CrystAge_Ma,Halfwidth_um,HeAge_Ma_raw,HeAge_Ma_sigma_raw,HeAge_Ma_sigma_5pct,HeAge_Ma_sigma_10pct,U238_ppm,Th232_ppm,HeAge_reference
04RF1zA,zircon,3500,64,770,15,38.5,77,84,67,Guenthner et al. 2013
04SC1zB,zircon,1779,71,659,51,32.95,65.9,415,111,Guenthner et al. 2013
04SC1zA,zircon,1779,47,649,13,32.45,64.9,422,110,Guenthner et al. 2013
04GF1zB,zircon,3500,63,638,12.6,31.9,63.8,254,104,Guenthner et al. 2013
04GF1zA,zircon,3500,56,620,18.4,31,62,247,88,Guenthner et al. 2013
04RF1zB,zircon,3500,72,557,10,27.85,55.7,155,177,Guenthner et al. 2013
04MT1zA,zircon,3500,67,545,12,27.25,54.5,362,89,Guenthner et al. 2013
04R1zB,zircon,1783,51,500,10,25,50,580,174,Guenthner et al. 2013
04MT1zB,zircon,3500,68,493,43,24.65,49.3,346,78,Guenthner et al. 2013
04R1zA,zircon,1783,42,357,7,17.85,35.7,656,167,Guenthner et al. 2013
04MT1zF,zircon,3500,54,329,7,16.45,32.9,567,75,Miltich 2005
04EQ1zA,zircon,3500,41,253,5.2,12.65,25.3,781,203,Guenthner et al. 2013
04EQ1zF,zircon,3500,43,241,4.6,12.05,24.1,647,284,Guenthner et al. 2013
04SH1zB,zircon,2603,56,225,4.6,11.25,22.5,883,184,Guenthner et al. 2013
04EQ1zB,zircon,3500,42,225,6,11.25,22.5,748,218,Guenthner et al. 2013
04MT1zD,zircon,3500,80,217,5,10.85,21.7,739,141,Miltich 2005
04MT1zC,zircon,3500,45,193,4,9.65,19.3,980,252,Miltich 2005
04MT1zE,zircon,3500,49,190,4,9.5,19,974,197,Miltich 2005
04SG1zB,zircon,3500,55,72,2,3.6,36,1738,1171,Guenthner et al. 2013
04EQ1zC,zircon,3500,36,57,1.1,2.85,28.5,894,246,Guenthner et al. 2013
04EQ1zE,zircon,3500,35,42,1.1,2.1,21,1107,351,Guenthner et al. 2013
04EQ1zD,zircon,3500,31,29,0.8,1.45,14.5,866,166,Guenthner et al. 2013
04SH1zA,zircon,2603,50,11,0.2,0.55,5.5,1011,228,Guenthner et al. 2013

0 comments on commit 5a61cdc

Please sign in to comment.