diff --git a/examples/ZrnHeInversionVartCryst.jl b/examples/ZrnHeInversionVartCryst.jl index eccca44..4e54100 100755 --- a/examples/ZrnHeInversionVartCryst.jl +++ b/examples/ZrnHeInversionVartCryst.jl @@ -33,151 +33,169 @@ # Literature samples from Guenthner et al. 2013 (AJS), Minnesota name = "MinnesotaInversion" - data = importdataset("minnesota.csv", ',', importas=:Tuple) + ds = importdataset("minnesota.csv", ',', importas=:Tuple) ## --- Prepare problem - # How many steps of the Markov chain should we run? - nsteps = 500000 - - # How long should we wait for MC to converge (become stationary) - burnin = 250000 - - simplified = false # Prefer simpler tT paths? - - dt = 10 # time step size in Myr - dTmax = 10.0 # Maximum reheating/burial per model timestep (to prevent numerical overflow) - - # Other model parameters - CrystAgeMax = 4000.0 # Ma -- forbid anything older than this - TStart = 400.0 # Temperature (in C) - dr = 1 # Radius step, in microns - - diffusionparams = (; + burnin = 100000 + model = ( + nsteps = 200000, # How many steps of the Markov chain should we run? + burnin = burnin, # 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 = 10.0, # Maximum reheating/burial per model timestep (to prevent numerical overflow) + TInit = 400.0, # Initial model temperature (in C) (i.e., crystallization temperature) + ΔTInit = -50.0, # TInit can vary from TInit to TInit+ΔTinit + TNow = 0.0, # Current surface temperature (in C) + ΔTNow = 10.0, # TNow may vary from TNow to TNow+ΔTNow + tInitMax = 4000.0, # Ma -- forbid anything older than this + maxPoints = 40, # Maximum allowed number of t-T points + simplified = false, # Prefer simpler tT paths? + # Diffusion parameters DzEa = 165.0, # kJ/mol DzD0 = 193188.0, # cm^2/sec DN17Ea = 71.0, # kJ/mol DN17D0 = 0.0034, #6.367E-3 # cm^2/sec + # Model uncertainty is not well known (depends on annealing parameters, + # decay constants, diffusion parameters, etc.), but is certainly non-zero. + # Here we add (in quadrature) a blanket model uncertainty of 25 Ma. + σModel = 25.0, # Ma + σAnnealing = 35.0, # Initial annealing uncertainty [Ma] + λAnnealing = 10 ./ burnin # Annealing decay [1/n] + ) + + + # 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 concentration, in PPM + HeAge = ds.HeAge_Ma_raw, # He age, in Ma + HeAge_sigma = ds.HeAge_Ma_sigma_raw, # He age uncertainty (1-sigma), in Ma + CrystAge = ds.CrystAge_Ma, # Crystallization age, in Ma ) - # Model uncertainty is not well known (depends on annealing parameters, - # decay constants, diffusion parameters, etc.), but is certainly non-zero. - # Here we add (in quadrature) a blanket model uncertainty of 25 Ma. - simannealparams = ( - 25.0, # Model uncertainty [Ma] - 35.0, # Initial uncertainty [Ma] - 10 ./ burnin, # lambda [1/n] + # Sort out crystallization ages and start time + map!(x->max(x, model.tInitMax), data.CrystAge, data.CrystAge) + model = (model..., + tInit = ceil(maximum(data.CrystAge)/model.dt) * model.dt, + ) + + tSteps = Array{Float64,1}(0+model.dt/2 : model.dt : model.tInit-model.dt/2) + ageSteps = reverse(tSteps) + + # Boundary conditions (e.g. 10C at present and 650 C at the time of zircon formation). + boundary = ( + agePoints = Float64[model.TNow, model.tInit], # Ma + TPoints = Float64[model.TNow, model.TInit], # Degrees C + T₀ = Float64[model.TNow, model.TInit], + ΔT = Float64[model.ΔTNow, model.ΔTInit], ) # Default: No unconformity is imposed - unconf_agePoints = Float64[] - unconf_TPoints = Float64[] + unconf = ( + agePoints = Float64[], # Ma + TPoints = Float64[], # Degrees C + ) # # Uncomment this section if you wish to impose an unconformity at any point in the record - # unconf_agePoints = Float64[600.0] # Ma - # unconf_TPoints = Float64[20.0] # Degrees C - # unconf_params = Float64[560,80,0,50] # Age0, DAge, T0, DT - # name = "$(name)Imposed-lowff" - - # Populate local variables from data frame with specified options - Halfwidth = data.Halfwidth_um - Uppm = data.U238_ppm - Thppm = data.Th232_ppm - HeAge = data.HeAge_Ma_raw - HeAge_sigma = data.HeAge_Ma_sigma_raw - CrystAge = data.CrystAge_Ma - - CrystAge[CrystAge .> CrystAgeMax] .= CrystAgeMax - tStart = ceil(maximum(CrystAge)/dt) * dt - - tSteps = Array{Float64,1}(0+dt/2 : dt : tStart-dt/2) - ageSteps = reverse(tSteps) - ntSteps = length(tSteps) # Number of time steps - eU = Uppm+.238*Thppm # Used only for plotting + # # Uniform distributions from Age₀ to Age₀+ΔAge, T₀ to T₀+ΔT, + # unconf = (; + # agePoints = Float64[560.0,], # Ma + # TPoints = Float64[20.0,], # Degrees C + # Age₀ = Float64[500,], + # ΔAge = Float64[80,], + # T₀ = Float64[0,], + # ΔT = Float64[50,], + # ) ## --- Test proscribed t-T paths with Neoproterozoic exhumation step # # Generate T path to test # Tr = 150 # T0 = 30 - # agePoints = Float64[tStart, tStart*29/30, 720, 580, 250, 0] # Age (Ma) - # TPoints = Float64[TStart, Tr+T0, Tr+T0, T0, 70, 10] # Temp. (C) + # agePoints = Float64[model.tInit, model.tInit*29/30, 720, 580, 250, 0] # Age (Ma) + # TPoints = Float64[model.TInit, Tr+T0, Tr+T0, T0, 70, 10] # Temp. (C) # TSteps = linterp1s(agePoints,TPoints,ageSteps) # # # Plot t-T path # plot(ageSteps,TSteps,xflip=true) # # # Calculate model ages - # CalcHeAges = Array{Float64}(undef, size(HeAge)) - # @time pr = DamageAnnealing(dt,tSteps,TSteps) - # zircons = Array{Zircon{Float64}}(undef, length(Halfwidth)) + # calcHeAges = Array{Float64}(undef, size(data.HeAge)) + # @time pr = anneal(model.dt, tSteps, TSteps, :zrdaam) + # zircons = Array{Zircon{Float64}}(undef, length(data.halfwidth)) # @time for i=1:length(zircons) # # Iterate through each grain, calculate the modeled age for each - # first_index = 1 + floor(Int64,(tStart - CrystAge[i])/dt) - # zircons[i] = Zircon(Halfwidth[i], dr, Uppm[i], Thppm[i], dt, ageSteps[first_index:end]) - # CalcHeAges[i] = HeAgeSpherical(zircons[i], @views(TSteps[first_index:end]), @views(pr[first_index:end,first_index:end]), diffusionparams) + # first_index = 1 + floor(Int64,(model.tInit - data.CrystAge[i])/model.dt) + # zircons[i] = Zircon(data.halfwidth[i], model.dr, data.U[i], data.Th[i], model.dt, ageSteps[first_index:end]) + # calcHeAges[i] = HeAgeSpherical(zircons[i], @views(TSteps[first_index:end]), @views(pr[first_index:end,first_index:end]), model) # end # # # Plot Comparison of results - # p2 = plot(eU, CalcHeAges, seriestype=:scatter,label="Model") - # plot!(p2, eU, HeAge, yerror=HeAge_sigma*2, seriestype=:scatter, label="Data") + # eU = data.U+.238*data.Th # Used only for plotting + # p2 = plot(eU, calcHeAges, seriestype=:scatter,label="Model") + # plot!(p2, eU, data.HeAge, yerror=data.HeAge_sigma*2, seriestype=:scatter, label="Data") # xlabel!(p2,"eU"); ylabel!(p2,"Age (Ma)") # display(p2) # # # Check log likelihood - # ll = sum(-(CalcHeAges - HeAge).^2 ./ (2 .* HeAge_sigma.^2)) - # ll = sum(-(CalcHeAges - HeAge).^2 ./ (2 .* AnnealedSigma.^2)) + # ll = sum(-(calcHeAges - data.HeAge).^2 ./ (2 .* data.HeAge_sigma.^2)) + # ll = sum(-(calcHeAges - data.HeAge).^2 ./ (2 .* σₐ.^2)) ## --- Invert for maximum likelihood t-T path - # Boundary conditions (10C at present and 650 C at the time of zircon - # formation). Optional: Apppend required unconformities at base of Cambrian, - # or elsewhere - boundary_agePoints = Float64[0, tStart] # Ma - boundary_TPoints = Float64[0, TStart] # Degrees C + # Add additional vectors for proposed unconformity and boundary points + unconf = (unconf..., + agePointsₚ = similar(unconf.agePoints), + TPointsₚ = similar(unconf.TPoints), + ) + boundary = (boundary..., + agePointsₚ = similar(boundary.agePoints), + TPointsₚ = similar(boundary.TPoints), + ) # This is where the "transdimensional" part comes in - nPoints = 0 - maxPoints = 40 # Maximum allowed number of age points - agePoints = Array{Float64}(undef, maxPoints+1) # Array of fixed size to hold all optional age points - TPoints = Array{Float64}(undef, maxPoints+1) # Array of fixed size to hold all optional age points - + agePoints = Array{Float64}(undef, model.maxPoints+1) # Array of fixed size to hold all optional age points + TPoints = Array{Float64}(undef, model.maxPoints+1) # Array of fixed size to hold all optional age points # Fill some intermediate points to give the MCMC something to work with - Tr = 250 # Residence temperature of initial proposal (value should not matter too much) - for t in [tStart/30,tStart/4,tStart/2,tStart-tStart/4,tStart-tStart/30] - global nPoints += 1 - agePoints[nPoints] = t # Ma - TPoints[nPoints] = Tr # Degrees C - end + Tr = 250 # Residence temperature + nPoints = 5 + agePoints[1:nPoints] .= [model.tInit/30,model.tInit/4,model.tInit/2,model.tInit-model.tInit/4,model.tInit-model.tInit/30] # Ma + TPoints[1:nPoints] .+ Tr # Degrees C - # # (Optional) Start with something close to the expected path - #Tr = 150 - #T0 = 30 - #agePoints[1:4] = Float64[tStart*29/30, 720, 510, 250]) # Age (Ma) - #TPoints[1:4] = Float64[ Tr+T0, Tr+T0, T0, 70]) # Temp. (C) - #nPoints+=4 + # # (Optional) Initialize with something close to the expected path + # Tr = 150 + # T0 = 30 + # nPoints = 4 + # agePoints[1:nPoints] = Float64[model.tInit*29/30, 720, 510, 250]) # Age (Ma) + # TPoints[1:nPoints] = Float64[ Tr+T0, Tr+T0, T0, 70]) # Temp. (C) - function MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints, unconf_TPoints, boundary_agePoints, boundary_TPoints, simannealparams, diffusionparams) + function MCMC_vartcryst(data, model, nPoints, agePoints, TPoints, unconf, boundary) # Calculate model ages for initial proposal - TSteps = linterp1s([view(agePoints, 1:nPoints) ; boundary_agePoints ; unconf_agePoints], - [view(TPoints, 1:nPoints) ; boundary_TPoints ; unconf_TPoints], ageSteps) - CalcHeAges = Array{Float64}(undef, size(HeAge)) - pr = DamageAnnealing(dt,tSteps,TSteps) # Damage annealing history - zircons = Array{Zircon{Float64}}(undef, length(Halfwidth)) + TSteps = linterp1s([view(agePoints, 1:nPoints) ; boundary.agePoints ; unconf.agePoints], + [view(TPoints, 1:nPoints) ; boundary.TPoints ; unconf.TPoints], ageSteps) + calcHeAges = Array{Float64}(undef, size(data.HeAge)) + pr = anneal(model.dt, tSteps, TSteps, :zrdaam) # Damage annealing history + + # Prepare a Mineral object for each analysis + zircons = Array{Zircon{Float64}}(undef, length(data.halfwidth)) for i=1:length(zircons) # Iterate through each grain, calculate the modeled age for each - first_index = 1 + floor(Int64,(tStart - CrystAge[i])/dt) - zircons[i] = Zircon(Halfwidth[i], dr, Uppm[i], Thppm[i], dt, ageSteps[first_index:end]) - CalcHeAges[i] = HeAgeSpherical(zircons[i], @views(TSteps[first_index:end]), @views(pr[first_index:end,first_index:end]), diffusionparams) + first_index = 1 + floor(Int64,(model.tInit - data.CrystAge[i])/model.dt) + zircons[i] = Zircon(data.halfwidth[i], model.dr, data.U[i], data.Th[i], model.dt, ageSteps[first_index:end]) + calcHeAges[i] = HeAgeSpherical(zircons[i], @views(TSteps[first_index:end]), @views(pr[first_index:end,first_index:end]), model) end - AnnealedSigma = simannealsigma.(1, HeAge_sigma; params=simannealparams) - UnAnnealedSigma = simannealsigma.(nsteps, HeAge_sigma; params=simannealparams) + + # Simulated annealing of uncertainty + σₐ = simannealsigma.(1, data.HeAge_sigma; params=model) + σₙ = simannealsigma.(model.nsteps, data.HeAge_sigma; params=model) # Log-likelihood for initial proposal - ll = normpdf_ll(HeAge, AnnealedSigma, CalcHeAges) - if simplified + ll = normpdf_ll(data.HeAge, σₐ, calcHeAges) + if model.simplified ll -= log(nPoints) end @@ -187,39 +205,36 @@ agePointsₚ = similar(agePoints) TPointsₚ = similar(TPoints) TStepsₚ = similar(TSteps) - unconf_agePointsₚ = similar(unconf_agePoints) - unconf_TPointsₚ = similar(unconf_TPoints) - boundary_TPointsₚ = similar(boundary_TPoints) - CalcHeAgesₚ = similar(CalcHeAges) + calcHeAgesₚ = similar(calcHeAges) # Distributions to populate - HeAgeDist = Array{Float64}(undef, length(HeAge), nsteps) - TStepDist = Array{Float64}(undef, ntSteps, nsteps) - llDist = Array{Float64}(undef, nsteps) - nDist = zeros(Int, nsteps) - acceptanceDist = zeros(Bool, nsteps) + HeAgeDist = Array{Float64}(undef, length(data.HeAge), model.nsteps) + TStepDist = Array{Float64}(undef, length(tSteps), model.nsteps) + llDist = Array{Float64}(undef, model.nsteps) + nDist = zeros(Int, model.nsteps) + acceptanceDist = zeros(Bool, model.nsteps) # Standard deviations of Gaussian proposal distributions for temperature and time - t_sigma = tStart/60 - T_sigma = TStart/60 + t_sigma = model.tInit/60 + T_sigma = model.TInit/60 # Proposal probabilities (must sum to 1) - move = 0.64 + move = 0.65 birth = 0.15 death = 0.15 # Should equal birth - boundary = 0.06 + boundarymove = 0.05 maxattempts = 1000 - @showprogress 10 "Running MCMC..." for n=1:nsteps + @showprogress 10 "Running MCMC..." for n=1:model.nsteps # Copy proposal from last accepted solution nPointsₚ = nPoints copyto!(agePointsₚ, agePoints) copyto!(TPointsₚ, TPoints) copyto!(TStepsₚ, TSteps) - copyto!(unconf_agePointsₚ, unconf_agePoints) - copyto!(unconf_TPointsₚ, unconf_TPoints) - copyto!(boundary_TPointsₚ, boundary_TPoints) + copyto!(unconf.agePointsₚ, unconf.agePoints) + copyto!(unconf.TPointsₚ, unconf.TPoints) + copyto!(boundary.TPointsₚ, boundary.TPoints) # Adjust the proposal r = rand() @@ -229,46 +244,46 @@ k = ceil(Int, rand() * nPoints) agePointsₚ[k] += randn() * t_sigma - if agePointsₚ[k] < dt + if agePointsₚ[k] < model.dt # Don't let any point get too close to 0 - agePointsₚ[k] += (dt - agePointsₚ[k]) - elseif agePointsₚ[k] > (tStart - dt) - # Don't let any point get too close to tStart - agePointsₚ[k] -= (agePointsₚ[k] - (tStart - dt)) + agePointsₚ[k] += (model.dt - agePointsₚ[k]) + elseif agePointsₚ[k] > (model.tInit - model.dt) + # Don't let any point get too close to tInit + agePointsₚ[k] -= (agePointsₚ[k] - (model.tInit - model.dt)) end # Move the Temperature of one model point if TPointsₚ[k] < 0 # Don't allow T<0 TPointsₚ[k] = 0 - elseif TPointsₚ[k] > TStart - # Don't allow T>TStart - TPointsₚ[k] = TStart + elseif TPointsₚ[k] > model.TInit + # Don't allow T>TInit + TPointsₚ[k] = model.TInit end # Interpolate proposed t-T path - TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary_agePoints ; unconf_agePointsₚ], - [view(TPointsₚ, 1:nPointsₚ) ; boundary_TPointsₚ ; unconf_TPointsₚ], ageSteps) + TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary.agePoints ; unconf.agePointsₚ], + [view(TPointsₚ, 1:nPointsₚ) ; boundary.TPointsₚ ; unconf.TPointsₚ], ageSteps) # Accept the proposal (and break out of the loop) if it satisfies the maximum reheating rate - maximum(diff(TStepsₚ)) < dTmax && break + maximum(diff(TStepsₚ)) < model.dTmax && break # Copy last accepted solution to re-modify if we don't break copyto!(agePointsₚ, agePoints) copyto!(TPointsₚ, TPoints) end - elseif (r < move+birth) && (nPointsₚ < maxPoints) + elseif (r < move+birth) && (nPointsₚ < model.maxPoints) # Birth: add a new model point nPointsₚ += 1 for i=1:maxattempts # Try maxattempts times to satisfy the reheating rate limit - agePointsₚ[nPointsₚ] = rand()*tStart - TPointsₚ[nPointsₚ] = rand()*TStart + agePointsₚ[nPointsₚ] = rand()*model.tInit + TPointsₚ[nPointsₚ] = rand()*model.TInit # Interpolate proposed t-T path - TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary_agePoints ; unconf_agePointsₚ], - [view(TPointsₚ, 1:nPointsₚ) ; boundary_TPointsₚ ; unconf_TPointsₚ], ageSteps) + TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary.agePoints ; unconf.agePointsₚ], + [view(TPointsₚ, 1:nPointsₚ) ; boundary.TPointsₚ ; unconf.TPointsₚ], ageSteps) # Accept the proposal (and break out of the loop) if it satisfies the maximum reheating rate - maximum(diff(TStepsₚ)) < dTmax && break + maximum(diff(TStepsₚ)) < model.dTmax && break end elseif (r < move+birth+death) && (r >= move+birth) && (nPointsₚ > 1) # Death: remove a model point @@ -279,59 +294,55 @@ TPointsₚ[k] = TPointsₚ[nPoints] # Interpolate proposed t-T path - TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary_agePoints ; unconf_agePointsₚ], - [view(TPointsₚ, 1:nPointsₚ) ; boundary_TPointsₚ ; unconf_TPointsₚ], ageSteps) + TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary.agePoints ; unconf.agePointsₚ], + [view(TPointsₚ, 1:nPointsₚ) ; boundary.TPointsₚ ; unconf.TPointsₚ], ageSteps) # Accept the proposal (and break out of the loop) if it satisfies the maximum reheating rate - maximum(diff(TStepsₚ)) < dTmax && break + maximum(diff(TStepsₚ)) < model.dTmax && break end else # Move boundary conditions for i=1:maxattempts # Try maxattempts times to satisfy the reheating rate limit - r2 = rand() - if r2 < 0.5 - # Allow the present temperature to vary from 0 to 10 degrees C - boundary_TPointsₚ[1] = 0+rand()*10 - else - # Allow the initial temperature to vary from TStart to TStart-50 C - boundary_TPointsₚ[2] = TStart-rand()*50 - end - if length(unconf_agePointsₚ) > 0 - # If there's an imposed unconformity, adjust within parameters - unconf_agePointsₚ = unconf_params[1] + rand()*unconf_params[2] - unconf_TPointsₚ = unconf_params[3] + rand()*unconf_params[4] + # Move the temperatures of the starting and ending boundaries + @. boundary.TPointsₚ = boundary.T₀ + rand()*boundary.ΔT + + # If there's an imposed unconformity, adjust within parameters + if length(unconf.agePoints) > 0 + @. unconf.agePointsₚ = unconf.Age₀ + rand()*unconf.ΔAge + @. unconf.TPointsₚ = unconf.T₀ + rand()*unconf.ΔT end + # Recalculate interpolated proposed t-T path - TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary_agePoints ; unconf_agePointsₚ], - [view(TPointsₚ, 1:nPointsₚ) ; boundary_TPointsₚ ; unconf_TPointsₚ], ageSteps) + TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary.agePoints ; unconf.agePointsₚ], + [view(TPointsₚ, 1:nPointsₚ) ; boundary.TPointsₚ ; unconf.TPointsₚ], ageSteps) # Accept the proposal (and break out of the loop) if it satisfies the maximum reheating rate - maximum(diff(TStepsₚ)) < dTmax && break + maximum(diff(TStepsₚ)) < model.dTmax && break # Copy last accepted solution to re-modify if we don't break - copyto!(unconf_agePointsₚ, unconf_agePoints) - copyto!(unconf_TPointsₚ, unconf_TPoints) - copyto!(boundary_TPointsₚ, boundary_TPoints) + copyto!(unconf.agePointsₚ, unconf.agePoints) + copyto!(unconf.TPointsₚ, unconf.TPoints) + copyto!(boundary.TPointsₚ, boundary.TPoints) end end # Recalculate interpolated proposed t-T path - TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary_agePoints ; unconf_agePointsₚ], - [view(TPointsₚ, 1:nPointsₚ) ; boundary_TPointsₚ ; unconf_TPointsₚ], ageSteps) + TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary.agePoints ; unconf.agePointsₚ], + [view(TPointsₚ, 1:nPointsₚ) ; boundary.TPointsₚ ; unconf.TPointsₚ], ageSteps) # Calculate model ages for each grain - DamageAnnealing!(pr, dt, tSteps, TStepsₚ) + anneal!(pr, model.dt, tSteps, TStepsₚ, :zrdaam) for i=1:length(zircons) - first_index = 1 + floor(Int64,(tStart - CrystAge[i])/dt) - CalcHeAgesₚ[i] = HeAgeSpherical(zircons[i], @views(TStepsₚ[first_index:end]), @views(pr[first_index:end,first_index:end]), diffusionparams) + first_index = 1 + floor(Int64,(model.tInit - data.CrystAge[i])/model.dt) + calcHeAgesₚ[i] = HeAgeSpherical(zircons[i], @views(TStepsₚ[first_index:end]), @views(pr[first_index:end,first_index:end]), model) end # Calculate log likelihood of proposal - AnnealedSigma .= simannealsigma.(n, HeAge_sigma; params=simannealparams) - llₚ = normpdf_ll(HeAge, AnnealedSigma, CalcHeAgesₚ) - llₗ = normpdf_ll(HeAge, AnnealedSigma, CalcHeAges) # Recalulate last one too with new AnnealedSigma - if simplified # slightly penalize more complex t-T paths + σₐ .= simannealsigma.(n, data.HeAge_sigma; params=model) + llₚ = normpdf_ll(data.HeAge, σₐ, calcHeAgesₚ) + llₗ = normpdf_ll(data.HeAge, σₐ, calcHeAges) # Recalulate last one too with new σₐ + if model.simplified # slightly penalize more complex t-T paths llₚ -= log(nPointsₚ) llₗ -= log(nPoints) end @@ -345,10 +356,10 @@ nPoints = nPointsₚ copyto!(agePoints, agePointsₚ) copyto!(TPoints, TPointsₚ) - copyto!(unconf_agePoints, unconf_agePointsₚ) - copyto!(unconf_TPoints, unconf_TPointsₚ) - copyto!(boundary_TPoints, boundary_TPointsₚ) - copyto!(CalcHeAges, CalcHeAgesₚ) + copyto!(unconf.agePoints, unconf.agePointsₚ) + copyto!(unconf.TPoints, unconf.TPointsₚ) + copyto!(boundary.TPoints, boundary.TPointsₚ) + copyto!(calcHeAges, calcHeAgesₚ) # These are saved for ouput, but not critical to the function of the MCMC loop copyto!(TSteps, TStepsₚ) @@ -356,9 +367,9 @@ end # Record results for analysis and troubleshooting - llDist[n] = normpdf_ll(HeAge, UnAnnealedSigma, CalcHeAges) # Recalculated to constant baseline + llDist[n] = normpdf_ll(data.HeAge, σₙ, calcHeAges) # Recalculated to constant baseline nDist[n] = nPoints # Distribution of # of points - HeAgeDist[:,n] = CalcHeAges # Distribution of He ages + HeAgeDist[:,n] = calcHeAges # Distribution of He ages # This is the actual output we want -- the distribution of t-T paths (t path is always identical) TStepDist[:,n] = TSteps # Distribution of T paths @@ -368,32 +379,32 @@ end # Run Markov Chain - @time (TStepDist, HeAgeDist, nDist, llDist, acceptanceDist) = MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints, unconf_TPoints, boundary_agePoints, boundary_TPoints, simannealparams, diffusionparams) + @time (TStepDist, HeAgeDist, nDist, llDist, acceptanceDist) = MCMC_vartcryst(data, model, nPoints, agePoints, TPoints, unconf, boundary) # # Save results using JLD - @save string(name, ".jld") ageSteps tSteps TStepDist burnin nsteps TStart tStart + @save string(name, ".jld") ageSteps tSteps TStepDist model - # Plot sample age-eU correlations - h = scatter(eU,HeAgeDist[:,burnin:burnin+50], label="") - plot!(h, eU, HeAge, yerror=HeAge_sigma, seriestype=:scatter, label="Data") +## --- Plot sample age-eU correlations + + eU = data.U+.238*data.Th # Used only for plotting + h = scatter(eU,HeAgeDist[:,model.burnin:model.burnin+50], label="") + plot!(h, eU, data.HeAge, yerror=data.HeAge_sigma, seriestype=:scatter, label="Data") xlabel!(h,"eU (ppm)"); ylabel!(h,"Age (Ma)") savefig(h,string(name,"Age-eU.pdf")) - # - ## --- Create image of paths # Resize the post-burnin part of the stationary distribution - TStepDistResized = Array{Float64}(undef, 2001, size(TStepDist,2)-burnin) - xq = collect(range(0,tStart,length=2001)) - for i=1:size(TStepDist,2)-burnin - TStepDistResized[:,i] = linterp1s(tSteps,TStepDist[:,i+burnin],xq) + TStepDistResized = Array{Float64}(undef, 2001, size(TStepDist,2)-model.burnin) + xq = collect(range(0,model.tInit,length=2001)) + for i=1:size(TStepDist,2)-model.burnin + TStepDistResized[:,i] = linterp1s(tSteps,TStepDist[:,i+model.burnin],xq) end # Calculate composite image - tTimage = zeros(ceil(Int, TStart)*2, size(TStepDistResized,1)) - yq = collect(0:0.5:TStart) + tTimage = zeros(ceil(Int, TInit)*2, size(TStepDistResized,1)) + yq = collect(0:0.5:TInit) for i=1:size(TStepDistResized,1) hist = fit(Histogram,TStepDistResized[i,:],yq,closed=:right) tTimage[:,i] = hist.weights diff --git a/src/Thermochron.jl b/src/Thermochron.jl index 325ca0f..e6da0b6 100644 --- a/src/Thermochron.jl +++ b/src/Thermochron.jl @@ -6,7 +6,7 @@ module Thermochron using ProgressMeter: @showprogress include("minerals.jl") - include("ZrnHe.jl") + include("zirconhelium.jl") include("inversion.jl") end diff --git a/src/inversion.jl b/src/inversion.jl index a720bd2..6693b0a 100644 --- a/src/inversion.jl +++ b/src/inversion.jl @@ -15,9 +15,8 @@ sigma = sqrt(sigma_analytical^2 + sigma_annealing^2) """ - function simannealsigma(n, sigma_analytical; params=(25.0, 35.0, 10/10^5)) - σₘ, σᵢ, λ = params - sigma_annealing = σᵢ*exp(-λ*n) + σₘ - return sqrt(sigma_analytical^2 + sigma_annealing^2) + function simannealsigma(n, sigma_analytical; params=(σModel=25.0, σAnnealing=35.0, λAnnealing=10/10^5)) + sigma_combined = params.σAnnealing * exp(-params.λAnnealing*n) + params.σModel + return sqrt(sigma_analytical^2 + sigma_combined^2) end export simannealsigma diff --git a/src/ZrnHe.jl b/src/zirconhelium.jl similarity index 92% rename from src/ZrnHe.jl rename to src/zirconhelium.jl index a5434ed..d611021 100644 --- a/src/ZrnHe.jl +++ b/src/zirconhelium.jl @@ -1,25 +1,26 @@ """ ```julia -ρᵣ = DamageAnnealing(dt::Number, tSteps::Vector, TSteps::Matrix) +ρᵣ = anneal(dt::Number, tSteps::Vector, TSteps::Matrix, [model::Symbol]) ``` Zircon damage annealing model as in Guenthner et al. 2013 (AJS) """ -function DamageAnnealing(dt::Number, tSteps::DenseVector, TSteps::DenseVector) +function anneal(dt::Number, tSteps::DenseVector, TSteps::DenseVector, model=:zrdaam) # Allocate matrix to hold reduced track lengths for all previous timesteps ρᵣ = zeros(length(tSteps),length(tSteps)) # In=-place version - DamageAnnealing!(ρᵣ, dt, tSteps, TSteps) + anneal!(ρᵣ, dt, tSteps, TSteps, model) end -export DamageAnnealing +export anneal """ ```julia -DamageAnnealing!(ρᵣ::Matrix, dt::Number, tSteps::Vector, TSteps::Vector) +anneal!(ρᵣ::Matrix, dt::Number, tSteps::Vector, TSteps::Vector, [model::Symbol]) ``` -In-place version of `DamageAnnealing` +In-place version of `anneal` """ -function DamageAnnealing!(ρᵣ::DenseMatrix, dt::Number, tSteps::DenseVector, TSteps::DenseVector) +anneal!(ρᵣ::DenseMatrix, dt::Number, tSteps::DenseVector, TSteps::DenseVector, model::Symbol) = anneal!(ρᵣ, dt, tSteps, TSteps, Val(model)) +function anneal!(ρᵣ::DenseMatrix, dt::Number, tSteps::DenseVector, TSteps::DenseVector, ::Val{:zrdaam}) # Annealing model constants B=-0.05721 C0=6.24534 @@ -63,7 +64,7 @@ function DamageAnnealing!(ρᵣ::DenseMatrix, dt::Number, tSteps::DenseVector, T return ρᵣ end -export DamageAnnealing! +export anneal! """ diff --git a/test/runtests.jl b/test/runtests.jl index 2bc3995..8be438b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,6 +2,6 @@ using Thermochron using LinearAlgebra using Test -@testset "Zircon helium" begin include("testZrnHe.jl") end +@testset "Zircon helium" begin include("testzirconhelium.jl") end @testset "Inversion" begin include("testinversion.jl") end @testset "Integrated Examples" begin include("testexamples.jl") end diff --git a/test/testexamples.jl b/test/testexamples.jl index ab4313c..d41734b 100644 --- a/test/testexamples.jl +++ b/test/testexamples.jl @@ -5,117 +5,125 @@ using StatGeochem # Read in data from file using StatGeochem datapath = joinpath("..", "examples", "minnesota.csv") -data = importdataset(datapath, ',', importas=:Tuple) - -# How many steps of the Markov chain should we run? -nsteps = 200 - -# How long should we wait for MC to converge (become stationary) -burnin = 100 +ds = importdataset(datapath, ',', importas=:Tuple); ## --- Prepare problem -simplified = false # Prefer simpler tT paths? - -dt = 10 # time step size in Myr -dTmax = 10.0 # Maximum reheating/burial per model timestep (to prevent numerical overflow) - -# Other model parameters -CrystAgeMax = 4000.0 # Ma -- forbid anything older than this -TStart = 400.0 # Temperature (in C) -dr = 1 # Radius step, in microns - -diffusionparams = (; +burnin = 100 +model = ( + nsteps = 200, # How many steps of the Markov chain should we run? + burnin = burnin, # 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 = 10.0, # Maximum reheating/burial per model timestep (to prevent numerical overflow) + TInit = 400.0, # Initial model temperature (in C) (i.e., crystallization temperature) + ΔTInit = -50.0, # TInit can vary from TInit to TInit+ΔTinit + TNow = 0.0, # Current surface temperature (in C) + ΔTNow = 10.0, # TNow may vary from TNow to TNow+ΔTNow + tInitMax = 4000.0, # Ma -- forbid anything older than this + maxPoints = 40, # Maximum allowed number of t-T points + simplified = false, # Prefer simpler tT paths? + # Diffusion parameters DzEa = 165.0, # kJ/mol DzD0 = 193188.0, # cm^2/sec DN17Ea = 71.0, # kJ/mol DN17D0 = 0.0034, #6.367E-3 # cm^2/sec + # Model uncertainty is not well known (depends on annealing parameters, + # decay constants, diffusion parameters, etc.), but is certainly non-zero. + # Here we add (in quadrature) a blanket model uncertainty of 25 Ma. + σModel = 25.0, # Ma + σAnnealing = 35.0, # Initial annealing uncertainty [Ma] + λAnnealing = 10 ./ burnin # Annealing decay [1/n] ) -# Model uncertainty is not well known (depends on annealing parameters, -# decay constants, diffusion parameters, etc.), but is certainly non-zero. -# Here we add (in quadrature) a blanket model uncertainty of 25 Ma. -simannealparams = ( - 25.0, # Model uncertainty [Ma] - 35.0, # Initial uncertainty [Ma] - 10 ./ burnin, # lambda [1/n] + +# 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 concentration, in PPM + HeAge = ds.HeAge_Ma_raw, # He age, in Ma + HeAge_sigma = ds.HeAge_Ma_sigma_raw, # He age uncertainty (1-sigma), in Ma + CrystAge = ds.CrystAge_Ma, # Crystallization age, in Ma ) -# Default: No unconformity is imposed -unconf_agePoints = Float64[] -unconf_TPoints = Float64[] +# Sort out crystallization ages and start time +map!(x->max(x, model.tInitMax), data.CrystAge, data.CrystAge) +model = (model..., + tInit = ceil(maximum(data.CrystAge)/model.dt) * model.dt, +) -# # Uncomment this section if you wish to impose an unconformity at any point in the record -# unconf_agePoints = Float64[600.0] # Ma -# unconf_TPoints = Float64[20.0] # Degrees C -# unconf_params = Float64[560,80,0,50] # Age0, DAge, T0, DT -# name = "$(name)Imposed-lowff" - -# Populate local variables from data frame with specified options -Halfwidth = data.Halfwidth_um -Uppm = data.U238_ppm -Thppm = data.Th232_ppm -HeAge = data.HeAge_Ma_raw -HeAge_sigma = data.HeAge_Ma_sigma_raw -CrystAge = data.CrystAge_Ma - -CrystAge[CrystAge .> CrystAgeMax] .= CrystAgeMax -tStart = ceil(maximum(CrystAge)/dt) * dt - -tSteps = Array{Float64,1}(0+dt/2 : dt : tStart-dt/2) +tSteps = Array{Float64,1}(0+model.dt/2 : model.dt : model.tInit-model.dt/2) ageSteps = reverse(tSteps) -ntSteps = length(tSteps) # Number of time steps -eU = Uppm+.238*Thppm # Used only for plotting +# Boundary conditions (e.g. 10C at present and 650 C at the time of zircon formation). +boundary = ( + agePoints = Float64[model.TNow, model.tInit], # Ma + TPoints = Float64[model.TNow, model.TInit], # Degrees C + T₀ = Float64[model.TNow, model.TInit], + ΔT = Float64[model.ΔTNow, model.ΔTInit], +) +# Default: No unconformity is imposed +unconf = ( + agePoints = Float64[], # Ma + TPoints = Float64[], # Degrees C +) + +# # Uncomment this section if you wish to impose an unconformity at any point in the record +# # Uniform distributions from Age₀ to Age₀+ΔAge, T₀ to T₀+ΔT, +# unconf = (; +# agePoints = Float64[560.0,], # Ma +# TPoints = Float64[20.0,], # Degrees C +# Age₀ = Float64[500,], +# ΔAge = Float64[80,], +# T₀ = Float64[0,], +# ΔT = Float64[50,], +# ) ## --- Invert for maximum likelihood t-T path -# Boundary conditions (10C at present and 650 C at the time of zircon -# formation). Optional: Apppend required unconformities at base of Cambrian, -# or elsewhere -boundary_agePoints = Float64[0, tStart] # Ma -boundary_TPoints = Float64[0, TStart] # Degrees C +# Add additional vectors for proposed unconformity and boundary points +unconf = (unconf..., + agePointsₚ = similar(unconf.agePoints), + TPointsₚ = similar(unconf.TPoints), +) +boundary = (boundary..., + agePointsₚ = similar(boundary.agePoints), + TPointsₚ = similar(boundary.TPoints), +) # This is where the "transdimensional" part comes in -nPoints = 0 -maxPoints = 40 # Maximum allowed number of age points -agePoints = Array{Float64}(undef, maxPoints+1) # Array of fixed size to hold all optional age points -TPoints = Array{Float64}(undef, maxPoints+1) # Array of fixed size to hold all optional age points - +agePoints = Array{Float64}(undef, model.maxPoints+1) # Array of fixed size to hold all optional age points +TPoints = Array{Float64}(undef, model.maxPoints+1) # Array of fixed size to hold all optional age points # Fill some intermediate points to give the MCMC something to work with -Tr = 250 # Residence temperature of initial proposal (value should not matter too much) -for t in [tStart/30,tStart/4,tStart/2,tStart-tStart/4,tStart-tStart/30] - global nPoints += 1 - agePoints[nPoints] = t # Ma - TPoints[nPoints] = Tr # Degrees C -end - -# # (Optional) Start with something close to the expected path -#Tr = 150 -#T0 = 30 -#agePoints[1:4] = Float64[tStart*29/30, 720, 510, 250]) # Age (Ma) -#TPoints[1:4] = Float64[ Tr+T0, Tr+T0, T0, 70]) # Temp. (C) -#nPoints+=4 +Tr = 250 # Residence temperature +nPoints = 5 +agePoints[1:nPoints] .= [model.tInit/30,model.tInit/4,model.tInit/2,model.tInit-model.tInit/4,model.tInit-model.tInit/30] # Ma +TPoints[1:nPoints] .+ Tr # Degrees C -function MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints, unconf_TPoints, boundary_agePoints, boundary_TPoints, simannealparams, diffusionparams) +function MCMC_vartcryst(data, model, nPoints, agePoints, TPoints, unconf, boundary) # Calculate model ages for initial proposal - TSteps = linterp1s([view(agePoints, 1:nPoints) ; boundary_agePoints ; unconf_agePoints], - [view(TPoints, 1:nPoints) ; boundary_TPoints ; unconf_TPoints], ageSteps) - CalcHeAges = Array{Float64}(undef, size(HeAge)) - pr = DamageAnnealing(dt,tSteps,TSteps) # Damage annealing history - zircons = Array{Zircon{Float64}}(undef, length(Halfwidth)) + TSteps = linterp1s([view(agePoints, 1:nPoints) ; boundary.agePoints ; unconf.agePoints], + [view(TPoints, 1:nPoints) ; boundary.TPoints ; unconf.TPoints], ageSteps) + calcHeAges = Array{Float64}(undef, size(data.HeAge)) + pr = anneal(model.dt, tSteps, TSteps, :zrdaam) # Damage annealing history + + # Prepare a Mineral object for each analysis + zircons = Array{Zircon{Float64}}(undef, length(data.halfwidth)) for i=1:length(zircons) # Iterate through each grain, calculate the modeled age for each - first_index = 1 + floor(Int64,(tStart - CrystAge[i])/dt) - zircons[i] = Zircon(Halfwidth[i], dr, Uppm[i], Thppm[i], dt, ageSteps[first_index:end]) - CalcHeAges[i] = HeAgeSpherical(zircons[i], @views(TSteps[first_index:end]), @views(pr[first_index:end,first_index:end]), diffusionparams) + first_index = 1 + floor(Int64,(model.tInit - data.CrystAge[i])/model.dt) + zircons[i] = Zircon(data.halfwidth[i], model.dr, data.U[i], data.Th[i], model.dt, ageSteps[first_index:end]) + calcHeAges[i] = HeAgeSpherical(zircons[i], @views(TSteps[first_index:end]), @views(pr[first_index:end,first_index:end]), model) end - AnnealedSigma = simannealsigma.(1, HeAge_sigma; params=simannealparams) - UnAnnealedSigma = simannealsigma.(nsteps, HeAge_sigma; params=simannealparams) + + # Simulated annealing of uncertainty + σₐ = simannealsigma.(1, data.HeAge_sigma; params=model) + σₙ = simannealsigma.(model.nsteps, data.HeAge_sigma; params=model) # Log-likelihood for initial proposal - ll = normpdf_ll(HeAge, AnnealedSigma, CalcHeAges) - if simplified + ll = normpdf_ll(data.HeAge, σₐ, calcHeAges) + if model.simplified ll -= log(nPoints) end @@ -125,39 +133,36 @@ function MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints agePointsₚ = similar(agePoints) TPointsₚ = similar(TPoints) TStepsₚ = similar(TSteps) - unconf_agePointsₚ = similar(unconf_agePoints) - unconf_TPointsₚ = similar(unconf_TPoints) - boundary_TPointsₚ = similar(boundary_TPoints) - CalcHeAgesₚ = similar(CalcHeAges) + calcHeAgesₚ = similar(calcHeAges) # Distributions to populate - HeAgeDist = Array{Float64}(undef, length(HeAge), nsteps) - TStepDist = Array{Float64}(undef, ntSteps, nsteps) - llDist = Array{Float64}(undef, nsteps) - nDist = zeros(Int, nsteps) - acceptanceDist = zeros(Bool, nsteps) + HeAgeDist = Array{Float64}(undef, length(data.HeAge), model.nsteps) + TStepDist = Array{Float64}(undef, length(tSteps), model.nsteps) + llDist = Array{Float64}(undef, model.nsteps) + nDist = zeros(Int, model.nsteps) + acceptanceDist = zeros(Bool, model.nsteps) # Standard deviations of Gaussian proposal distributions for temperature and time - t_sigma = tStart/60 - T_sigma = TStart/60 + t_sigma = model.tInit/60 + T_sigma = model.TInit/60 # Proposal probabilities (must sum to 1) - move = 0.64 + move = 0.65 birth = 0.15 death = 0.15 # Should equal birth - boundary = 0.06 + boundarymove = 0.05 maxattempts = 1000 - @showprogress 10 "Running MCMC..." for n=1:nsteps + @showprogress 10 "Running MCMC..." for n=1:model.nsteps # Copy proposal from last accepted solution nPointsₚ = nPoints copyto!(agePointsₚ, agePoints) copyto!(TPointsₚ, TPoints) copyto!(TStepsₚ, TSteps) - copyto!(unconf_agePointsₚ, unconf_agePoints) - copyto!(unconf_TPointsₚ, unconf_TPoints) - copyto!(boundary_TPointsₚ, boundary_TPoints) + copyto!(unconf.agePointsₚ, unconf.agePoints) + copyto!(unconf.TPointsₚ, unconf.TPoints) + copyto!(boundary.TPointsₚ, boundary.TPoints) # Adjust the proposal r = rand() @@ -167,46 +172,46 @@ function MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints k = ceil(Int, rand() * nPoints) agePointsₚ[k] += randn() * t_sigma - if agePointsₚ[k] < dt + if agePointsₚ[k] < model.dt # Don't let any point get too close to 0 - agePointsₚ[k] += (dt - agePointsₚ[k]) - elseif agePointsₚ[k] > (tStart - dt) - # Don't let any point get too close to tStart - agePointsₚ[k] -= (agePointsₚ[k] - (tStart - dt)) + agePointsₚ[k] += (model.dt - agePointsₚ[k]) + elseif agePointsₚ[k] > (model.tInit - model.dt) + # Don't let any point get too close to tInit + agePointsₚ[k] -= (agePointsₚ[k] - (model.tInit - model.dt)) end # Move the Temperature of one model point if TPointsₚ[k] < 0 # Don't allow T<0 TPointsₚ[k] = 0 - elseif TPointsₚ[k] > TStart - # Don't allow T>TStart - TPointsₚ[k] = TStart + elseif TPointsₚ[k] > model.TInit + # Don't allow T>TInit + TPointsₚ[k] = model.TInit end # Interpolate proposed t-T path - TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary_agePoints ; unconf_agePointsₚ], - [view(TPointsₚ, 1:nPointsₚ) ; boundary_TPointsₚ ; unconf_TPointsₚ], ageSteps) + TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary.agePoints ; unconf.agePointsₚ], + [view(TPointsₚ, 1:nPointsₚ) ; boundary.TPointsₚ ; unconf.TPointsₚ], ageSteps) # Accept the proposal (and break out of the loop) if it satisfies the maximum reheating rate - maximum(diff(TStepsₚ)) < dTmax && break + maximum(diff(TStepsₚ)) < model.dTmax && break # Copy last accepted solution to re-modify if we don't break copyto!(agePointsₚ, agePoints) copyto!(TPointsₚ, TPoints) end - elseif (r < move+birth) && (nPointsₚ < maxPoints) + elseif (r < move+birth) && (nPointsₚ < model.maxPoints) # Birth: add a new model point nPointsₚ += 1 for i=1:maxattempts # Try maxattempts times to satisfy the reheating rate limit - agePointsₚ[nPointsₚ] = rand()*tStart - TPointsₚ[nPointsₚ] = rand()*TStart + agePointsₚ[nPointsₚ] = rand()*model.tInit + TPointsₚ[nPointsₚ] = rand()*model.TInit # Interpolate proposed t-T path - TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary_agePoints ; unconf_agePointsₚ], - [view(TPointsₚ, 1:nPointsₚ) ; boundary_TPointsₚ ; unconf_TPointsₚ], ageSteps) + TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary.agePoints ; unconf.agePointsₚ], + [view(TPointsₚ, 1:nPointsₚ) ; boundary.TPointsₚ ; unconf.TPointsₚ], ageSteps) # Accept the proposal (and break out of the loop) if it satisfies the maximum reheating rate - maximum(diff(TStepsₚ)) < dTmax && break + maximum(diff(TStepsₚ)) < model.dTmax && break end elseif (r < move+birth+death) && (r >= move+birth) && (nPointsₚ > 1) # Death: remove a model point @@ -217,59 +222,55 @@ function MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints TPointsₚ[k] = TPointsₚ[nPoints] # Interpolate proposed t-T path - TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary_agePoints ; unconf_agePointsₚ], - [view(TPointsₚ, 1:nPointsₚ) ; boundary_TPointsₚ ; unconf_TPointsₚ], ageSteps) + TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary.agePoints ; unconf.agePointsₚ], + [view(TPointsₚ, 1:nPointsₚ) ; boundary.TPointsₚ ; unconf.TPointsₚ], ageSteps) # Accept the proposal (and break out of the loop) if it satisfies the maximum reheating rate - maximum(diff(TStepsₚ)) < dTmax && break + maximum(diff(TStepsₚ)) < model.dTmax && break end else # Move boundary conditions for i=1:maxattempts # Try maxattempts times to satisfy the reheating rate limit - r2 = rand() - if r2 < 0.5 - # Allow the present temperature to vary from 0 to 10 degrees C - boundary_TPointsₚ[1] = 0+rand()*10 - else - # Allow the initial temperature to vary from TStart to TStart-50 C - boundary_TPointsₚ[2] = TStart-rand()*50 - end - if length(unconf_agePointsₚ) > 0 - # If there's an imposed unconformity, adjust within parameters - unconf_agePointsₚ = unconf_params[1] + rand()*unconf_params[2] - unconf_TPointsₚ = unconf_params[3] + rand()*unconf_params[4] + # Move the temperatures of the starting and ending boundaries + @. boundary.TPointsₚ = boundary.T₀ + rand()*boundary.ΔT + + # If there's an imposed unconformity, adjust within parameters + if length(unconf.agePoints) > 0 + @. unconf.agePointsₚ = unconf.Age₀ + rand()*unconf.ΔAge + @. unconf.TPointsₚ = unconf.T₀ + rand()*unconf.ΔT end + # Recalculate interpolated proposed t-T path - TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary_agePoints ; unconf_agePointsₚ], - [view(TPointsₚ, 1:nPointsₚ) ; boundary_TPointsₚ ; unconf_TPointsₚ], ageSteps) + TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary.agePoints ; unconf.agePointsₚ], + [view(TPointsₚ, 1:nPointsₚ) ; boundary.TPointsₚ ; unconf.TPointsₚ], ageSteps) # Accept the proposal (and break out of the loop) if it satisfies the maximum reheating rate - maximum(diff(TStepsₚ)) < dTmax && break + maximum(diff(TStepsₚ)) < model.dTmax && break # Copy last accepted solution to re-modify if we don't break - copyto!(unconf_agePointsₚ, unconf_agePoints) - copyto!(unconf_TPointsₚ, unconf_TPoints) - copyto!(boundary_TPointsₚ, boundary_TPoints) + copyto!(unconf.agePointsₚ, unconf.agePoints) + copyto!(unconf.TPointsₚ, unconf.TPoints) + copyto!(boundary.TPointsₚ, boundary.TPoints) end end # Recalculate interpolated proposed t-T path - TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary_agePoints ; unconf_agePointsₚ], - [view(TPointsₚ, 1:nPointsₚ) ; boundary_TPointsₚ ; unconf_TPointsₚ], ageSteps) + TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary.agePoints ; unconf.agePointsₚ], + [view(TPointsₚ, 1:nPointsₚ) ; boundary.TPointsₚ ; unconf.TPointsₚ], ageSteps) # Calculate model ages for each grain - DamageAnnealing!(pr, dt, tSteps, TStepsₚ) + anneal!(pr, model.dt, tSteps, TStepsₚ, :zrdaam) for i=1:length(zircons) - first_index = 1 + floor(Int64,(tStart - CrystAge[i])/dt) - CalcHeAgesₚ[i] = HeAgeSpherical(zircons[i], @views(TStepsₚ[first_index:end]), @views(pr[first_index:end,first_index:end]), diffusionparams) + first_index = 1 + floor(Int64,(model.tInit - data.CrystAge[i])/model.dt) + calcHeAgesₚ[i] = HeAgeSpherical(zircons[i], @views(TStepsₚ[first_index:end]), @views(pr[first_index:end,first_index:end]), model) end # Calculate log likelihood of proposal - AnnealedSigma .= simannealsigma.(n, HeAge_sigma; params=simannealparams) - llₚ = normpdf_ll(HeAge, AnnealedSigma, CalcHeAgesₚ) - llₗ = normpdf_ll(HeAge, AnnealedSigma, CalcHeAges) # Recalulate last one too with new AnnealedSigma - if simplified # slightly penalize more complex t-T paths + σₐ .= simannealsigma.(n, data.HeAge_sigma; params=model) + llₚ = normpdf_ll(data.HeAge, σₐ, calcHeAgesₚ) + llₗ = normpdf_ll(data.HeAge, σₐ, calcHeAges) # Recalulate last one too with new σₐ + if model.simplified # slightly penalize more complex t-T paths llₚ -= log(nPointsₚ) llₗ -= log(nPoints) end @@ -283,10 +284,10 @@ function MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints nPoints = nPointsₚ copyto!(agePoints, agePointsₚ) copyto!(TPoints, TPointsₚ) - copyto!(unconf_agePoints, unconf_agePointsₚ) - copyto!(unconf_TPoints, unconf_TPointsₚ) - copyto!(boundary_TPoints, boundary_TPointsₚ) - copyto!(CalcHeAges, CalcHeAgesₚ) + copyto!(unconf.agePoints, unconf.agePointsₚ) + copyto!(unconf.TPoints, unconf.TPointsₚ) + copyto!(boundary.TPoints, boundary.TPointsₚ) + copyto!(calcHeAges, calcHeAgesₚ) # These are saved for ouput, but not critical to the function of the MCMC loop copyto!(TSteps, TStepsₚ) @@ -294,9 +295,9 @@ function MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints end # Record results for analysis and troubleshooting - llDist[n] = normpdf_ll(HeAge, UnAnnealedSigma, CalcHeAges) # Recalculated to constant baseline + llDist[n] = normpdf_ll(data.HeAge, σₙ, calcHeAges) # Recalculated to constant baseline nDist[n] = nPoints # Distribution of # of points - HeAgeDist[:,n] = CalcHeAges # Distribution of He ages + HeAgeDist[:,n] = calcHeAges # Distribution of He ages # This is the actual output we want -- the distribution of t-T paths (t path is always identical) TStepDist[:,n] = TSteps # Distribution of T paths @@ -305,9 +306,10 @@ function MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints return (TStepDist, HeAgeDist, nDist, llDist, acceptanceDist) end + # Run Markov Chain -@time (TStepDist, HeAgeDist, nDist, llDist, acceptanceDist) = MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints, unconf_TPoints, boundary_agePoints, boundary_TPoints, simannealparams, diffusionparams) -@time (TStepDist, HeAgeDist, nDist, llDist, acceptanceDist) = MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints, unconf_TPoints, boundary_agePoints, boundary_TPoints, simannealparams, diffusionparams) +@time (TStepDist, HeAgeDist, nDist, llDist, acceptanceDist) = MCMC_vartcryst(data, model, nPoints, agePoints, TPoints, unconf, boundary) +@time (TStepDist, HeAgeDist, nDist, llDist, acceptanceDist) = MCMC_vartcryst(data, model, nPoints, agePoints, TPoints, unconf, boundary) @test isa(TStepDist, AbstractMatrix) diff --git a/test/testinversion.jl b/test/testinversion.jl index a6cd9dd..61525a6 100644 --- a/test/testinversion.jl +++ b/test/testinversion.jl @@ -1,8 +1,8 @@ simannealparams = ( - 10., # Model uncertainty [Ma] - 100., # Initial uncertainty [Ma] - 10 / 10^5, # lambda [1/n] + σModel = 10., # Model uncertainty [Ma] + σAnnealing = 100., # Initial uncertainty [Ma] + λAnnealing = 10 / 10^5, # lambda [1/n] ) @test simannealsigma(1, 10, params=simannealparams) ≈ 110.44365174144839 diff --git a/test/testZrnHe.jl b/test/testzirconhelium.jl similarity index 99% rename from test/testZrnHe.jl rename to test/testzirconhelium.jl index 50b6781..9e39045 100644 --- a/test/testZrnHe.jl +++ b/test/testzirconhelium.jl @@ -33,7 +33,7 @@ pr_known = [0.000205838 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.00018167 0.000345328 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000175903 0.000306842 0.000575964 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000174306 0.00029782 0.000515098 0.000955418 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000173858 0.000295383 0.000501104 0.000859803 0.00157683 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000173735 0.000294718 0.00049742 0.000838273 0.00142761 0.00259 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000173702 0.000294541 0.000496445 0.000832761 0.00139476 0.00235863 0.00423469 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000173694 0.000294495 0.000496194 0.000831348 0.00138659 0.00230891 0.00387837 0.00689242 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000173692 0.000294484 0.000496131 0.000830996 0.00138456 0.00229693 0.00380376 0.00634758 0.0111652 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496117 0.000830911 0.00138408 0.00229407 0.00378637 0.00623665 0.0103389 0.0179917 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496113 0.000830892 0.00138397 0.00229341 0.00378237 0.00621171 0.0101757 0.0167508 0.0288076 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830888 0.00138394 0.00229326 0.00378149 0.0062062 0.0101404 0.0165134 0.0269675 0.0457439 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229323 0.0037813 0.00620503 0.0101329 0.0164642 0.0266275 0.0430615 0.0718099 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378126 0.0062048 0.0101314 0.0164542 0.02656 0.0425841 0.0679927 0.110913 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378126 0.00620476 0.0101311 0.0164523 0.026547 0.0424937 0.0673402 0.105664 0.167408 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378125 0.00620475 0.010131 0.0164519 0.0265447 0.0424772 0.0672227 0.104805 0.160534 0.244771 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378125 0.00620475 0.010131 0.0164519 0.0265443 0.0424744 0.0672025 0.104658 0.15946 0.236352 0.34328 0 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378125 0.00620475 0.010131 0.0164519 0.0265442 0.0424739 0.0671992 0.104635 0.159287 0.235101 0.333839 0.457702 0 0 0 0 0 0 0 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378125 0.00620475 0.010131 0.0164519 0.0265442 0.0424739 0.0671987 0.104631 0.159261 0.234913 0.332512 0.448193 0.577183 0 0 0 0 0 0 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378125 0.00620475 0.010131 0.0164519 0.0265442 0.0424739 0.0671987 0.104631 0.159257 0.234886 0.332326 0.446936 0.568685 0.688768 0 0 0 0 0 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378125 0.00620475 0.010131 0.0164519 0.0265442 0.0424739 0.0671987 0.104631 0.159257 0.234883 0.332302 0.446773 0.567636 0.682035 0.782538 0 0 0 0 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378125 0.00620475 0.010131 0.0164519 0.0265442 0.0424739 0.0671987 0.104631 0.159257 0.234883 0.332299 0.446753 0.56751 0.681264 0.777758 0.85451 0 0 0 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378125 0.00620475 0.010131 0.0164519 0.0265442 0.0424739 0.0671987 0.104631 0.159257 0.234883 0.332299 0.446751 0.567497 0.681181 0.777256 0.851417 0.905927 0 0 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378125 0.00620475 0.010131 0.0164519 0.0265442 0.0424739 0.0671987 0.104631 0.159257 0.234883 0.332299 0.446751 0.567496 0.681172 0.777207 0.851122 0.904069 0.940745 0 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378125 0.00620475 0.010131 0.0164519 0.0265442 0.0424739 0.0671987 0.104631 0.159257 0.234883 0.332299 0.446751 0.567495 0.681172 0.777203 0.851097 0.903911 0.939694 0.96343 0 0 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378125 0.00620475 0.010131 0.0164519 0.0265442 0.0424739 0.0671987 0.104631 0.159257 0.234883 0.332299 0.446751 0.567495 0.681172 0.777203 0.851095 0.903899 0.939615 0.962863 0.977805 0 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378125 0.00620475 0.010131 0.0164519 0.0265442 0.0424739 0.0671987 0.104631 0.159257 0.234883 0.332299 0.446751 0.567495 0.681172 0.777203 0.851095 0.903898 0.939609 0.962826 0.977511 0.986726 0 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378125 0.00620475 0.010131 0.0164519 0.0265442 0.0424739 0.0671987 0.104631 0.159257 0.234883 0.332299 0.446751 0.567495 0.681172 0.777203 0.851095 0.903898 0.939609 0.962824 0.977495 0.98658 0.992174 0 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378125 0.00620475 0.010131 0.0164519 0.0265442 0.0424739 0.0671987 0.104631 0.159257 0.234883 0.332299 0.446751 0.567495 0.681172 0.777203 0.851095 0.903898 0.939609 0.962824 0.977494 0.986574 0.992104 0.995453 0.0; 0.000173691 0.000294481 0.000496112 0.000830887 0.00138394 0.00229322 0.00378125 0.00620475 0.010131 0.0164519 0.0265442 0.0424739 0.0671987 0.104631 0.159257 0.234883 0.332299 0.446751 0.567495 0.681172 0.777203 0.851095 0.903898 0.939609 0.962824 0.977494 0.986573 0.992102 0.995421 0.9974] - pr = DamageAnnealing(dt,tSteps,TSteps) + pr = anneal(dt, tSteps, TSteps, :zrdaam) @test size(pr) == (30,30) @test round.(pr, sigdigits=6) ≈ pr_known