Additional tests on hiatus output

src/StratMetropolis.jl

@@ -712,7 +712,8 @@
                 end
 
                 # Move discontinuities to the nearest hiatus, if one not already present
-                for h in rand(closest_hiatus_unique, length(closest_hiatus_unique))
+                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
@@ -809,7 +810,8 @@
                 end
 
                 # Move discontinuities to the nearest hiatus, if one not already present
-                for h in rand(closest_hiatus_unique, length(closest_hiatus_unique))
+                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

test/testRadiocarbon.jl

@@ -59,7 +59,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 ≈ [8322.4, 8237.9, 8156.4, 8072.6, 8026.7, 7982.9, 7939.3, 7915.2, 7892.4, 7869.8, 7846.7, 7822.2] atol=20
+@test mdl.Age ≈ [8322.52, 8237.68, 8155.65, 8072.77, 8026.99, 7983.31, 7939.54, 7915.44, 7892.49, 7869.92, 7846.76, 7822.28] atol=20
 @test mdl.Age_025CI ≈ [8205.9, 8045.4, 7999.0, 7975.6, 7894.3, 7860.6, 7839.4, 7802.5, 7780.8, 7762.1, 7747.9, 7736.0] atol=35
 @test mdl.Age_975CI ≈ [8405.0, 8386.1, 8347.1, 8163.5, 8151.0, 8123.9, 8029.5, 8016.3, 8003.0, 7985.9, 7958.2, 7916.3] atol=35
 # Test that all age-depth models are in stratigraphic order
@@ -117,7 +117,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
 (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
 # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
@@ -148,8 +148,9 @@ hiatus.Duration_sigma = [  30.5,   20.0 ]
 
 # Test that results match expectation, within some tolerance
 @test mdl.Age isa Vector{Float64}
-@test mdl.Age ≈ [8328.02, 8284.17, 8099.46, 8057.85, 8023.46, 7991.04, 7958.17, 7945.3, 7932.21, 7821.86, 7809.53, 7797.44] atol=20
+@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 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)

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 = 100000 # Number of steps to run in distribution MCMC
+config.nsteps = 1000000 # 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
 
@@ -47,8 +47,11 @@ 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.73, 747.82, 728.84, 724.13, 720.83, 717.67, 714.49, 713.02, 711.48, 700.71, 699.33, 697.95] 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 ≈ [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 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)
+@test mean(hiatusdist, dims=2) ≈ [10.580012942504894; 18.96167245288326;;] atol=2
+@test mean(lldist) ≈ -6.149720449018242 atol=6