Fix bug in new hiatus-moving routine

brenhinkeller · May 27, 2024 · 1b76661 · 1b76661
1 parent bd39e4d
commit 1b76661
Show file tree

Hide file tree

Showing 3 changed files with 17 additions and 25 deletions.
diff --git a/src/StratMetropolis.jl b/src/StratMetropolis.jl
@@ -715,15 +715,11 @@
                 for _ in 1:(length(closest_hiatus_unique)÷2)
                     h = rand(closest_hiatus_unique)
                     if model_agesₚ[h-1] == model_agesₚ[h]
-                        n = findclosestunequal(model_agesₚ, h)
-                        if n < h
-                            @inbounds for i = n:h-1
-                                model_agesₚ[i] = model_agesₚ[n]
-                            end
-                        elseif n > h
-                            @inbounds for i = h:n
-                                model_agesₚ[i] = model_agesₚ[n]
-                            end
+                        iₙ = findclosestunequal(model_agesₚ, h)
+                        if iₙ < h
+                            model_agesₚ[iₙ:h-1] .= model_agesₚ[iₙ]
+                        elseif iₙ > h
+                            model_agesₚ[h:iₙ] .= model_agesₚ[iₙ]
                         end
                     end
                 end
@@ -813,15 +809,11 @@
                 for _ in 1:(length(closest_hiatus_unique)÷2)
                     h = rand(closest_hiatus_unique)
                     if model_agesₚ[h-1] == model_agesₚ[h]
-                        n = findclosestunequal(model_agesₚ, h)
-                        if n < h
-                            @inbounds for i = n:h-1
-                                model_agesₚ[i] = model_agesₚ[n]
-                            end
-                        elseif n > h
-                            @inbounds for i = h:n
-                                model_agesₚ[i] = model_agesₚ[n]
-                            end
+                        iₙ = findclosestunequal(model_agesₚ, h)
+                        if iₙ < h
+                            model_agesₚ[iₙ:h-1] .= model_agesₚ[iₙ]
+                        elseif iₙ > h
+                            model_agesₚ[h:iₙ] .= model_agesₚ[iₙ]
                         end
                     end
                 end

diff --git a/test/testRadiocarbon.jl b/test/testRadiocarbon.jl
@@ -149,8 +149,8 @@ hiatus.Duration_sigma = [  30.5,   20.0 ]
 # Test that results match expectation, within some tolerance
 @test mdl.Age isa Vector{Float64}
 @test mdl.Age ≈ [8327.67, 8283.87, 8099.66, 8057.74, 8023.23, 7990.66, 7957.71, 7944.6, 7931.19, 7823.06, 7810.29, 7797.89] atol=20
-@test mdl.Age_025CI ≈ [8202.92, 8138.36, 7974.42, 7960.03, 7896.94, 7874.23, 7858.76, 7842.18, 7827.89, 7717.17, 7709.2, 7703.73] atol=40
-@test mdl.Age_975CI ≈ [8405.45, 8391.43, 8225.43, 8163.38, 8151.17, 8123.45, 8030.52, 8023.59, 8016.78, 7913.76, 7904.32, 7877.73] atol=40
+@test mdl.Age_025CI ≈ [8210.55, 8152.11, 7983.06, 7968.49, 7913.95, 7885.02, 7868.98, 7852.67, 7838.42, 7730.2, 7719.48, 7710.82] atol=40
+@test mdl.Age_975CI ≈ [8405.89, 8391.62, 8225.61, 8163.41, 8151.31, 8123.92, 8030.35, 8023.65, 8016.86, 7916.1, 7905.95, 7886.48] atol=40
 # Test that all age-depth models are in stratigraphic order
 @test all([issorted(x, rev=true) for x in eachcol(agedist)])
 @test size(hiatusdist) == (nHiatuses, config.nsteps)
diff --git a/test/testStratOnly.jl b/test/testStratOnly.jl
@@ -16,7 +16,7 @@ config.resolution = 5 # Same units as sample height. Smaller is slower!
 config.bounding = 0.5 # how far away do we place runaway bounds, as a fraction of total section height. Larger is slower.
 (bottom, top) = extrema(smpl.Height)
 npoints_approx = round(Int,length(bottom:config.resolution:top) * (1 + 2*config.bounding))
-config.nsteps = 1000000 # Number of steps to run in distribution MCMC
+config.nsteps = 100000 # Number of steps to run in distribution MCMC
 config.burnin = 100000*npoints_approx # Number to discard
 config.sieve = round(Int,npoints_approx) # Record one out of every nsieve steps
 
@@ -26,7 +26,7 @@ config.sieve = round(Int,npoints_approx) # Record one out of every nsieve steps
 # Test that results match expectation, within some tolerance
 @test mdl.Age isa Vector{Float64}
 @test mdl.Age ≈ [751.87, 742.76, 733.65, 724.51, 720.15, 715.96, 711.8, 709.15, 706.6, 704.04, 701.42, 698.67] atol=1
-@test mdl.Age_025CI ≈ [742.48, 723.81, 718.7, 715.73, 707.42, 703.73, 701.34, 698.2, 696.35, 695.04, 694.0, 693.05] atol=2
+@test mdl.Age_025CI ≈ [742.48, 723.81, 718.7, 715.73, 707.42, 703.73, 701.34, 698.2, 696.35, 695.04, 694.0, 693.05] atol=3
 @test mdl.Age_975CI ≈ [761.17, 757.98, 752.64, 733.57, 731.21, 728.07, 722.24, 720.59, 718.6, 716.11, 712.48, 704.28] atol=3
 # Test that all age-depth models are in stratigraphic order
 @test all([issorted(x, rev=true) for x in eachcol(agedist)])
@@ -46,9 +46,9 @@ hiatus.Duration_sigma = [   3.1,    2.0 ]
 
 # Test that results match expectation, within some tolerance
 @test mdl.Age isa Vector{Float64}
-@test mdl.Age ≈ [752.76, 747.88, 728.92, 724.14, 720.8, 717.62, 714.44, 712.95, 711.41, 700.82, 699.4, 697.99] atol=3
-@test mdl.Age_025CI ≈ [743.41, 734.96, 717.54, 715.41, 709.67, 707.05, 705.26, 703.72, 702.49, 693.27, 692.65, 692.06] atol=4
-@test mdl.Age_975CI ≈ [761.77, 759.23, 741.37, 732.59, 730.64, 728.05, 723.52, 722.34, 720.9, 709.03, 706.78, 703.46] atol=4
+@test mdl.Age ≈ [752.76, 747.88, 728.92, 724.14, 720.8, 717.62, 714.44, 712.95, 711.41, 700.82, 699.4, 697.99] atol=1
+@test mdl.Age_025CI ≈ [743.41, 734.96, 717.54, 715.41, 709.67, 707.05, 705.26, 703.72, 702.49, 693.27, 692.65, 692.06] atol=3
+@test mdl.Age_975CI ≈ [761.77, 759.23, 741.37, 732.59, 730.64, 728.05, 723.52, 722.34, 720.9, 709.03, 706.78, 703.46] atol=3
 # Test that all age-depth models are in stratigraphic order
 @test all([issorted(x, rev=true) for x in eachcol(agedist)])
 @test size(hiatusdist) == (nHiatuses, config.nsteps)