diff --git a/Project.toml b/Project.toml
index 560eb7b5..8ad8a317 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,13 +4,16 @@ authors = ["CLIMA contributors <clima-software@caltech.edu>"]
 version = "0.6.0"
 
 [deps]
+AbstractGPs = "99985d1d-32ba-4be9-9821-2ec096f28918"
 AbstractMCMC = "80f14c24-f653-4e6a-9b94-39d6b0f70001"
 AdvancedMH = "5b7e9947-ddc0-4b3f-9b55-0d8042f74170"
 Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 EnsembleKalmanProcesses = "aa8a2aa5-91d8-4396-bcef-d4f2ec43552d"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9"
+KernelFunctions = "ec8451be-7e33-11e9-00cf-bbf324bd1392"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
@@ -19,25 +22,30 @@ ProgressBars = "49802e3a-d2f1-5c88-81d8-b72133a6f568"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RandomFeatures = "36c3bae2-c0c3-419d-b3b4-eebadd35c5e5"
+ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 ScikitLearn = "3646fa90-6ef7-5e7e-9f22-8aca16db6324"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
+AbstractGPs = "0.5.21"
 AbstractMCMC = "3.3, 4, 5"
 AdvancedMH = "0.6, 0.7, 0.8"
 Conda = "1.7"
 Distributions = "0.24, 0.25"
 DocStringExtensions = "0.8, 0.9"
 EnsembleKalmanProcesses = "2"
+ForwardDiff = "0.10.38"
 GaussianProcesses = "0.12"
+KernelFunctions = "0.10.64"
 MCMCChains = "4.14, 5, 6"
 Printf = "1"
 ProgressBars = "1"
 PyCall = "1.93"
 Random = "1"
 RandomFeatures = "0.3"
+ReverseDiff = "1.15.3"
 ScikitLearn = "0.6, 0.7"
 StableRNGs = "1"
 Statistics = "1"
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index f8c0a082..47f74e31 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -26,6 +26,12 @@ weakdeps = ["ChainRulesCore", "Test"]
     AbstractFFTsChainRulesCoreExt = "ChainRulesCore"
     AbstractFFTsTestExt = "Test"
 
+[[deps.AbstractGPs]]
+deps = ["ChainRulesCore", "Distributions", "FillArrays", "IrrationalConstants", "KernelFunctions", "LinearAlgebra", "PDMats", "Random", "RecipesBase", "Reexport", "Statistics", "StatsBase", "Test"]
+git-tree-sha1 = "010e6b5e06a9a6b1be6b146a7f29d5848490cc64"
+uuid = "99985d1d-32ba-4be9-9821-2ec096f28918"
+version = "0.5.21"
+
 [[deps.AbstractMCMC]]
 deps = ["BangBang", "ConsoleProgressMonitor", "Distributed", "LogDensityProblems", "Logging", "LoggingExtras", "ProgressLogging", "Random", "StatsBase", "TerminalLoggers", "Transducers"]
 git-tree-sha1 = "87e63dcb990029346b091b170252f3c416568afc"
@@ -71,9 +77,9 @@ uuid = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"
 version = "1.1.3"
 
 [[deps.ArgCheck]]
-git-tree-sha1 = "a3a402a35a2f7e0b87828ccabbd5ebfbebe356b4"
+git-tree-sha1 = "680b3b8759bd4c54052ada14e52355ab69e07876"
 uuid = "dce04be8-c92d-5529-be00-80e4d2c0e197"
-version = "2.3.0"
+version = "2.4.0"
 
 [[deps.ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
@@ -93,9 +99,9 @@ version = "3.5.1+1"
 
 [[deps.ArrayInterface]]
 deps = ["Adapt", "LinearAlgebra"]
-git-tree-sha1 = "d60a1922358aa203019b7857a2c8c37329b8736c"
+git-tree-sha1 = "017fcb757f8e921fb44ee063a7aafe5f89b86dd1"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.17.0"
+version = "7.18.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
@@ -169,10 +175,10 @@ uuid = "9718e550-a3fa-408a-8086-8db961cd8217"
 version = "0.1.1"
 
 [[deps.BenchmarkTools]]
-deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
-git-tree-sha1 = "f1dff6729bc61f4d49e140da1af55dcd1ac97b2f"
+deps = ["Compat", "JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
+git-tree-sha1 = "e38fbc49a620f5d0b660d7f543db1009fe0f8336"
 uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-version = "1.5.0"
+version = "1.6.0"
 
 [[deps.BitTwiddlingConvenienceFunctions]]
 deps = ["Static"]
@@ -184,7 +190,7 @@ version = "0.1.6"
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "8873e196c2eb87962a2048b3b8e08946535864a1"
 uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
-version = "1.0.8+2"
+version = "1.0.8+4"
 
 [[deps.CPUSummary]]
 deps = ["CpuId", "IfElse", "PrecompileTools", "Static"]
@@ -199,16 +205,16 @@ uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
 version = "1.18.2+1"
 
 [[deps.CalibrateEmulateSample]]
-deps = ["AbstractMCMC", "AdvancedMH", "Conda", "Distributions", "DocStringExtensions", "EnsembleKalmanProcesses", "GaussianProcesses", "LinearAlgebra", "MCMCChains", "Pkg", "Printf", "ProgressBars", "PyCall", "Random", "RandomFeatures", "ScikitLearn", "StableRNGs", "Statistics", "StatsBase"]
+deps = ["AbstractGPs", "AbstractMCMC", "AdvancedMH", "Conda", "Distributions", "DocStringExtensions", "EnsembleKalmanProcesses", "ForwardDiff", "GaussianProcesses", "KernelFunctions", "LinearAlgebra", "MCMCChains", "Pkg", "Printf", "ProgressBars", "PyCall", "Random", "RandomFeatures", "ReverseDiff", "ScikitLearn", "StableRNGs", "Statistics", "StatsBase"]
 path = ".."
 uuid = "95e48a1f-0bec-4818-9538-3db4340308e3"
-version = "0.5.3"
+version = "0.6.0"
 
 [[deps.ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra"]
-git-tree-sha1 = "3e4b134270b372f2ed4d4d0e936aabaefc1802bc"
+git-tree-sha1 = "1713c74e00545bfe14605d2a2be1712de8fbcb58"
 uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "1.25.0"
+version = "1.25.1"
 weakdeps = ["SparseArrays"]
 
     [deps.ChainRulesCore.extensions]
@@ -222,9 +228,9 @@ version = "0.1.13"
 
 [[deps.CodecBzip2]]
 deps = ["Bzip2_jll", "TranscodingStreams"]
-git-tree-sha1 = "e7c529cc31bb85b97631b922fa2e6baf246f5905"
+git-tree-sha1 = "84990fa864b7f2b4901901ca12736e45ee79068c"
 uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"
-version = "0.8.4"
+version = "0.8.5"
 
 [[deps.CodecZlib]]
 deps = ["TranscodingStreams", "Zlib_jll"]
@@ -292,9 +298,9 @@ weakdeps = ["IntervalSets", "LinearAlgebra", "StaticArrays"]
 
 [[deps.Convex]]
 deps = ["AbstractTrees", "BenchmarkTools", "LDLFactorizations", "LinearAlgebra", "MathOptInterface", "OrderedCollections", "SparseArrays", "Test"]
-git-tree-sha1 = "dac1878b4996fa56292d2c3bd28f2498b980bb93"
+git-tree-sha1 = "dec769959be3af9ba94970b1f14b31c196b0fb9e"
 uuid = "f65535da-76fb-5f13-bab9-19810c17039a"
-version = "0.16.3"
+version = "0.16.4"
 
 [[deps.CpuId]]
 deps = ["Markdown"]
@@ -369,9 +375,9 @@ version = "1.11.0"
 
 [[deps.Distributions]]
 deps = ["AliasTables", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"]
-git-tree-sha1 = "3101c32aab536e7a27b1763c0797dba151b899ad"
+git-tree-sha1 = "03aa5d44647eaec98e1920635cdfed5d5560a8b9"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
-version = "0.25.113"
+version = "0.25.117"
 
     [deps.Distributions.extensions]
     DistributionsChainRulesCoreExt = "ChainRulesCore"
@@ -391,9 +397,9 @@ version = "0.9.3"
 
 [[deps.Documenter]]
 deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "CodecZlib", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "TOML", "Test", "Unicode"]
-git-tree-sha1 = "5a1ee886566f2fa9318df1273d8b778b9d42712d"
+git-tree-sha1 = "d0ea2c044963ed6f37703cead7e29f70cba13d7e"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "1.7.0"
+version = "1.8.0"
 
 [[deps.Downloads]]
 deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
@@ -420,9 +426,9 @@ version = "2.0.1"
 
 [[deps.Expat_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "1c6317308b9dc757616f0b5cb379db10494443a7"
+git-tree-sha1 = "e51db81749b0777b2147fbe7b783ee79045b8e99"
 uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
-version = "2.6.2+0"
+version = "2.6.4+3"
 
 [[deps.FFMPEG]]
 deps = ["FFMPEG_jll"]
@@ -438,15 +444,15 @@ version = "4.4.4+1"
 
 [[deps.FFTW]]
 deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
-git-tree-sha1 = "4820348781ae578893311153d69049a93d05f39d"
+git-tree-sha1 = "7de7c78d681078f027389e067864a8d53bd7c3c9"
 uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-version = "1.8.0"
+version = "1.8.1"
 
 [[deps.FFTW_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "4d81ed14783ec49ce9f2e168208a12ce1815aa25"
 uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
-version = "3.3.10+1"
+version = "3.3.10+3"
 
 [[deps.FastGaussQuadrature]]
 deps = ["LinearAlgebra", "SpecialFunctions", "StaticArrays"]
@@ -472,9 +478,9 @@ weakdeps = ["PDMats", "SparseArrays", "Statistics"]
 
 [[deps.FiniteDiff]]
 deps = ["ArrayInterface", "LinearAlgebra", "Setfield"]
-git-tree-sha1 = "b10bdafd1647f57ace3885143936749d61638c3b"
+git-tree-sha1 = "84e3a47db33be7248daa6274b287507dd6ff84e8"
 uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
-version = "2.26.0"
+version = "2.26.2"
 
     [deps.FiniteDiff.extensions]
     FiniteDiffBandedMatricesExt = "BandedMatrices"
@@ -490,9 +496,9 @@ version = "2.26.0"
 
 [[deps.Fontconfig_jll]]
 deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Zlib_jll"]
-git-tree-sha1 = "db16beca600632c95fc8aca29890d83788dd8b23"
+git-tree-sha1 = "21fac3c77d7b5a9fc03b0ec503aa1a6392c34d2b"
 uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
-version = "2.13.96+0"
+version = "2.15.0+0"
 
 [[deps.Formatting]]
 deps = ["Logging", "Printf"]
@@ -512,15 +518,26 @@ weakdeps = ["StaticArrays"]
 
 [[deps.FreeType2_jll]]
 deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Zlib_jll"]
-git-tree-sha1 = "5c1d8ae0efc6c2e7b1fc502cbe25def8f661b7bc"
+git-tree-sha1 = "786e968a8d2fb167f2e4880baba62e0e26bd8e4e"
 uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
-version = "2.13.2+0"
+version = "2.13.3+1"
 
 [[deps.FriBidi_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "1ed150b39aebcc805c26b93a8d0122c940f64ce2"
+git-tree-sha1 = "846f7026a9decf3679419122b49f8a1fdb48d2d5"
 uuid = "559328eb-81f9-559d-9380-de523a88c83c"
-version = "1.0.14+0"
+version = "1.0.16+0"
+
+[[deps.FunctionWrappers]]
+git-tree-sha1 = "d62485945ce5ae9c0c48f124a84998d755bae00e"
+uuid = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
+version = "1.1.3"
+
+[[deps.Functors]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "64d8e93700c7a3f28f717d265382d52fac9fa1c1"
+uuid = "d9f16b24-f501-4c13-a1f2-28368ffc5196"
+version = "0.4.12"
 
 [[deps.Future]]
 deps = ["Random"]
@@ -553,27 +570,27 @@ version = "1.3.1"
 
 [[deps.Git_jll]]
 deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"]
-git-tree-sha1 = "ea372033d09e4552a04fd38361cd019f9003f4f4"
+git-tree-sha1 = "399f4a308c804b446ae4c91eeafadb2fe2c54ff9"
 uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb"
-version = "2.46.2+0"
+version = "2.47.1+0"
 
 [[deps.Glib_jll]]
 deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"]
-git-tree-sha1 = "674ff0db93fffcd11a3573986e550d66cd4fd71f"
+git-tree-sha1 = "b0036b392358c80d2d2124746c2bf3d48d457938"
 uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
-version = "2.80.5+0"
+version = "2.82.4+0"
 
 [[deps.Graphite2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011"
+git-tree-sha1 = "01979f9b37367603e2848ea225918a3b3861b606"
 uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472"
-version = "1.3.14+0"
+version = "1.3.14+1"
 
 [[deps.HarfBuzz_jll]]
 deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll"]
-git-tree-sha1 = "401e4f3f30f43af2c8478fc008da50096ea5240f"
+git-tree-sha1 = "55c53be97790242c29031e5cd45e8ac296dadda3"
 uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566"
-version = "8.3.1+0"
+version = "8.5.0+0"
 
 [[deps.HostCPUFeatures]]
 deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"]
@@ -612,9 +629,9 @@ version = "1.4.2"
 
 [[deps.IntelOpenMP_jll]]
 deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl"]
-git-tree-sha1 = "10bd689145d2c3b2a9844005d01087cc1194e79e"
+git-tree-sha1 = "0f14a5456bdc6b9731a5682f439a672750a09e48"
 uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
-version = "2024.2.1+0"
+version = "2025.0.4+0"
 
 [[deps.InteractiveUtils]]
 deps = ["Markdown"]
@@ -655,9 +672,9 @@ weakdeps = ["Dates", "Test"]
     InverseFunctionsTestExt = "Test"
 
 [[deps.InvertedIndices]]
-git-tree-sha1 = "0dc7b50b8d436461be01300fd8cd45aa0274b038"
+git-tree-sha1 = "6da3c4316095de0f5ee2ebd875df8721e7e0bdbe"
 uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
-version = "1.3.0"
+version = "1.3.1"
 
 [[deps.IrrationalConstants]]
 git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
@@ -676,9 +693,9 @@ version = "1.0.0"
 
 [[deps.JLLWrappers]]
 deps = ["Artifacts", "Preferences"]
-git-tree-sha1 = "be3dc50a92e5a386872a493a10050136d4703f9b"
+git-tree-sha1 = "a007feb38b422fbdab534406aeca1b86823cb4d6"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.6.1"
+version = "1.7.0"
 
 [[deps.JSON]]
 deps = ["Dates", "Mmap", "Parsers", "Unicode"]
@@ -686,12 +703,30 @@ git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a"
 uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 version = "0.21.4"
 
+[[deps.JSON3]]
+deps = ["Dates", "Mmap", "Parsers", "PrecompileTools", "StructTypes", "UUIDs"]
+git-tree-sha1 = "1d322381ef7b087548321d3f878cb4c9bd8f8f9b"
+uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
+version = "1.14.1"
+
+    [deps.JSON3.extensions]
+    JSON3ArrowExt = ["ArrowTypes"]
+
+    [deps.JSON3.weakdeps]
+    ArrowTypes = "31f734f8-188a-4ce0-8406-c8a06bd891cd"
+
 [[deps.KernelDensity]]
 deps = ["Distributions", "DocStringExtensions", "FFTW", "Interpolations", "StatsBase"]
 git-tree-sha1 = "7d703202e65efa1369de1279c162b915e245eed1"
 uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
 version = "0.6.9"
 
+[[deps.KernelFunctions]]
+deps = ["ChainRulesCore", "Compat", "CompositionsBase", "Distances", "FillArrays", "Functors", "IrrationalConstants", "LinearAlgebra", "LogExpFunctions", "Random", "Requires", "SpecialFunctions", "Statistics", "StatsBase", "TensorCore", "Test", "ZygoteRules"]
+git-tree-sha1 = "4a38fbd48503f2839b1d6033e0f7ffaec90dda33"
+uuid = "ec8451be-7e33-11e9-00cf-bbf324bd1392"
+version = "0.10.64"
+
 [[deps.LAME_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
 git-tree-sha1 = "170b660facf5df5de098d866564877e119141cbd"
@@ -712,9 +747,9 @@ version = "18.1.7+0"
 
 [[deps.LZO_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "854a9c268c43b77b0a27f22d7fab8d33cdb3a731"
+git-tree-sha1 = "1c602b1127f4751facb671441ca72715cc95938a"
 uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
-version = "2.10.2+1"
+version = "2.10.3+0"
 
 [[deps.LaTeXStrings]]
 git-tree-sha1 = "dda21b8cbd6a6c40d9d02a73230f9d70fed6918c"
@@ -774,9 +809,9 @@ version = "1.11.0"
 
 [[deps.Libffi_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "0b4a5d71f3e5200a7dff793393e09dfc2d874290"
+git-tree-sha1 = "27ecae93dd25ee0909666e6835051dd684cc035e"
 uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
-version = "3.2.2+1"
+version = "3.2.2+2"
 
 [[deps.Libgcrypt_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll"]
@@ -786,27 +821,27 @@ version = "1.11.0+0"
 
 [[deps.Libgpg_error_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "c6ce1e19f3aec9b59186bdf06cdf3c4fc5f5f3e6"
+git-tree-sha1 = "df37206100d39f79b3376afb6b9cee4970041c61"
 uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
-version = "1.50.0+0"
+version = "1.51.1+0"
 
 [[deps.Libiconv_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "61dfdba58e585066d8bce214c5a51eaa0539f269"
+git-tree-sha1 = "be484f5c92fad0bd8acfef35fe017900b0b73809"
 uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
-version = "1.17.0+1"
+version = "1.18.0+0"
 
 [[deps.Libmount_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "0c4f9c4f1a50d8f35048fa0532dabbadf702f81e"
+git-tree-sha1 = "89211ea35d9df5831fca5d33552c02bd33878419"
 uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
-version = "2.40.1+0"
+version = "2.40.3+0"
 
 [[deps.Libuuid_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "5ee6203157c120d79034c748a2acba45b82b8807"
+git-tree-sha1 = "e888ad02ce716b319e6bdb985d2ef300e7089889"
 uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
-version = "2.40.1+0"
+version = "2.40.3+0"
 
 [[deps.LineSearches]]
 deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"]
@@ -827,9 +862,9 @@ version = "2.1.2"
 
 [[deps.LogExpFunctions]]
 deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
-git-tree-sha1 = "a2d09619db4e765091ee5c6ffe8872849de0feea"
+git-tree-sha1 = "13ca9e2586b89836fd20cccf56e57e2b9ae7f38f"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-version = "0.3.28"
+version = "0.3.29"
 
     [deps.LogExpFunctions.extensions]
     LogExpFunctionsChainRulesCoreExt = "ChainRulesCore"
@@ -876,9 +911,9 @@ version = "0.2.1"
 
 [[deps.MKL_jll]]
 deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "oneTBB_jll"]
-git-tree-sha1 = "f046ccd0c6db2832a9f639e2c669c6fe867e5f4f"
+git-tree-sha1 = "ed4097130e3dd3721814b5f277da72f48905e80c"
 uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
-version = "2024.2.0+0"
+version = "2025.0.1+0"
 
 [[deps.MLJModelInterface]]
 deps = ["Random", "ScientificTypesBase", "StatisticalTraits"]
@@ -887,10 +922,9 @@ uuid = "e80e1ace-859a-464e-9ed9-23947d8ae3ea"
 version = "1.11.0"
 
 [[deps.MacroTools]]
-deps = ["Markdown", "Random"]
-git-tree-sha1 = "2fa9ee3e63fd3a4f7a9a4f4744a52f4856de82df"
+git-tree-sha1 = "72aebe0b5051e5143a079a4685a46da330a40472"
 uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
-version = "0.5.13"
+version = "0.5.15"
 
 [[deps.ManualMemory]]
 git-tree-sha1 = "bcaef4fc7a0cfe2cba636d84cda54b5e4e4ca3cd"
@@ -909,10 +943,10 @@ uuid = "d0879d2d-cac2-40c8-9cee-1863dc0c7391"
 version = "0.1.2"
 
 [[deps.MathOptInterface]]
-deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays", "SpecialFunctions", "Test", "Unicode"]
-git-tree-sha1 = "e065ca5234f53fd6f920efaee4940627ad991fb4"
+deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "DataStructures", "ForwardDiff", "JSON3", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays", "SpecialFunctions", "Test", "Unicode"]
+git-tree-sha1 = "f346e3b50c8bb62b7a5c4961ebdbaef65d1ce1e2"
 uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
-version = "1.34.0"
+version = "1.35.1"
 
 [[deps.MbedTLS_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -941,9 +975,9 @@ version = "2023.12.12"
 
 [[deps.MutableArithmetics]]
 deps = ["LinearAlgebra", "SparseArrays", "Test"]
-git-tree-sha1 = "90077f1e79de8c9c7c8a90644494411111f4e07b"
+git-tree-sha1 = "43122df26d27424b23577d59e2d8020f28386516"
 uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
-version = "1.5.2"
+version = "1.6.2"
 
 [[deps.NLSolversBase]]
 deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"]
@@ -953,9 +987,9 @@ version = "7.8.3"
 
 [[deps.NaNMath]]
 deps = ["OpenLibm_jll"]
-git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4"
+git-tree-sha1 = "fe891aea7ccd23897520db7f16931212454e277e"
 uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
-version = "1.0.2"
+version = "1.1.1"
 
 [[deps.NaturalSort]]
 git-tree-sha1 = "eda490d06b9f7c00752ee81cfa451efe55521e21"
@@ -967,9 +1001,9 @@ uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 version = "1.2.0"
 
 [[deps.OffsetArrays]]
-git-tree-sha1 = "1a27764e945a152f7ca7efa04de513d473e9542e"
+git-tree-sha1 = "5e1897147d1ff8d98883cda2be2187dcf57d8f0c"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.14.1"
+version = "1.15.0"
 weakdeps = ["Adapt"]
 
     [deps.OffsetArrays.extensions]
@@ -983,9 +1017,9 @@ version = "1.3.5+1"
 
 [[deps.OpenBLAS32_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "dd806c813429ff09878ea3eeb317818f3ca02871"
+git-tree-sha1 = "ece4587683695fe4c5f20e990da0ed7e83c351e7"
 uuid = "656ef2d0-ae68-5445-9ca0-591084a874a2"
-version = "0.3.28+3"
+version = "0.3.29+0"
 
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
@@ -1001,19 +1035,19 @@ version = "0.8.1+2"
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
 git-tree-sha1 = "7493f61f55a6cce7325f197443aa80d32554ba10"
 uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
-version = "3.0.15+1"
+version = "3.0.15+3"
 
 [[deps.OpenSpecFun_jll]]
-deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "1346c9208249809840c91b26703912dff463d335"
 uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
-version = "0.5.5+0"
+version = "0.5.6+0"
 
 [[deps.Optim]]
 deps = ["Compat", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"]
-git-tree-sha1 = "d9b79c4eed437421ac4285148fcadf42e0700e89"
+git-tree-sha1 = "ab7edad78cdef22099f43c54ef77ac63c2c9cc64"
 uuid = "429524aa-4258-5aef-a3af-852621145aeb"
-version = "1.9.4"
+version = "1.10.0"
 weakdeps = ["MathOptInterface"]
 
     [deps.Optim.extensions]
@@ -1026,9 +1060,9 @@ uuid = "91d4177d-7536-5919-b921-800302f37372"
 version = "1.3.3+0"
 
 [[deps.OrderedCollections]]
-git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5"
+git-tree-sha1 = "12f1439c4f986bb868acda6ea33ebc78e19b95ad"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.6.3"
+version = "1.7.0"
 
 [[deps.PCRE2_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -1037,9 +1071,9 @@ version = "10.42.0+1"
 
 [[deps.PDMats]]
 deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "949347156c25054de2db3b166c52ac4728cbad65"
+git-tree-sha1 = "966b85253e959ea89c53a9abebbf2e964fbf593b"
 uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
-version = "0.11.31"
+version = "0.11.32"
 
 [[deps.Parameters]]
 deps = ["OrderedCollections", "UnPack"]
@@ -1132,9 +1166,9 @@ uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
 version = "1.10.2"
 
 [[deps.PtrArrays]]
-git-tree-sha1 = "77a42d78b6a92df47ab37e177b2deac405e1c88f"
+git-tree-sha1 = "1d36ef11a9aaf1e8b74dacc6a731dd1de8fd493d"
 uuid = "43287f4e-b6f4-7ad1-bb20-aadabca52c3d"
-version = "1.2.1"
+version = "1.3.0"
 
 [[deps.PyCall]]
 deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"]
@@ -1210,6 +1244,12 @@ git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
 uuid = "ae029012-a4dd-5104-9daa-d747884805df"
 version = "1.3.0"
 
+[[deps.ReverseDiff]]
+deps = ["ChainRulesCore", "DiffResults", "DiffRules", "ForwardDiff", "FunctionWrappers", "LinearAlgebra", "LogExpFunctions", "MacroTools", "NaNMath", "Random", "SpecialFunctions", "StaticArrays", "Statistics"]
+git-tree-sha1 = "cc6cd622481ea366bb9067859446a8b01d92b468"
+uuid = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
+version = "1.15.3"
+
 [[deps.Rmath]]
 deps = ["Random", "Rmath_jll"]
 git-tree-sha1 = "f65dcb5fa46aee0cf9ed6274ccbd597adc49aa7b"
@@ -1224,9 +1264,9 @@ version = "0.4.3+0"
 
 [[deps.SCS]]
 deps = ["MathOptInterface", "Requires", "SCS_jll", "SparseArrays"]
-git-tree-sha1 = "0dfe49eaa058ce905a4199af379b8e411e6126e5"
+git-tree-sha1 = "aa3fcff53da363b4ba4b54d4ac4c9186ab00d703"
 uuid = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
-version = "2.0.1"
+version = "2.0.2"
 
     [deps.SCS.extensions]
     SCSSCS_GPU_jllExt = ["SCS_GPU_jll"]
@@ -1238,9 +1278,9 @@ version = "2.0.1"
 
 [[deps.SCS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl", "OpenBLAS32_jll"]
-git-tree-sha1 = "668bcf4b25cf992564321ccb70b205f9a7487cfa"
+git-tree-sha1 = "902cc4e42ecca21bbd74babf899b2a5b12add323"
 uuid = "f4f2fc5b-1d94-523c-97ea-2ab488bedf4b"
-version = "3.2.6+0"
+version = "3.2.7+0"
 
 [[deps.SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
@@ -1276,9 +1316,9 @@ version = "0.5.0"
 
 [[deps.SentinelArrays]]
 deps = ["Dates", "Random"]
-git-tree-sha1 = "d0553ce4031a081cc42387a9b9c8441b7d99f32d"
+git-tree-sha1 = "712fb0231ee6f9120e005ccd56297abbc053e7e0"
 uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
-version = "1.4.7"
+version = "1.4.8"
 
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -1312,9 +1352,9 @@ version = "1.11.0"
 
 [[deps.SpecialFunctions]]
 deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "2f5d4697f21388cbe1ff299430dd169ef97d7e14"
+git-tree-sha1 = "64cca0c26b4f31ba18f13f6c12af7c85f478cfde"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.4.0"
+version = "2.5.0"
 weakdeps = ["ChainRulesCore"]
 
     [deps.SpecialFunctions.extensions]
@@ -1351,9 +1391,9 @@ weakdeps = ["OffsetArrays", "StaticArrays"]
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "777657803913ffc7e8cc20f0fd04b634f871af8f"
+git-tree-sha1 = "47091a0340a675c738b1304b58161f3b0839d454"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.8"
+version = "1.9.10"
 weakdeps = ["ChainRulesCore", "Statistics"]
 
     [deps.StaticArrays.extensions]
@@ -1405,6 +1445,12 @@ git-tree-sha1 = "a6b1675a536c5ad1a60e5a5153e1fee12eb146e3"
 uuid = "892a3eda-7b42-436c-8928-eab12a02cf0e"
 version = "0.4.0"
 
+[[deps.StructTypes]]
+deps = ["Dates", "UUIDs"]
+git-tree-sha1 = "159331b30e94d7b11379037feeb9b690950cace8"
+uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
+version = "1.11.0"
+
 [[deps.StyledStrings]]
 uuid = "f489334b-da3d-4c2e-b8f0-e476e12c162b"
 version = "1.11.0"
@@ -1440,6 +1486,12 @@ deps = ["ArgTools", "SHA"]
 uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
 version = "1.10.0"
 
+[[deps.TensorCore]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6"
+uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50"
+version = "0.1.1"
+
 [[deps.TerminalLoggers]]
 deps = ["LeftChildRightSiblingTrees", "Logging", "Markdown", "Printf", "ProgressLogging", "UUIDs"]
 git-tree-sha1 = "f133fab380933d042f6796eda4e130272ba520ca"
@@ -1533,74 +1585,80 @@ version = "1.0.0"
 
 [[deps.XML2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"]
-git-tree-sha1 = "6a451c6f33a176150f315726eba8b92fbfdb9ae7"
+git-tree-sha1 = "a2fccc6559132927d4c5dc183e3e01048c6dcbd6"
 uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
-version = "2.13.4+0"
+version = "2.13.5+0"
 
 [[deps.XSLT_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "XML2_jll", "Zlib_jll"]
-git-tree-sha1 = "a54ee957f4c86b526460a720dbc882fa5edcbefc"
+git-tree-sha1 = "7d1671acbe47ac88e981868a078bd6b4e27c5191"
 uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
-version = "1.1.41+0"
+version = "1.1.42+0"
 
 [[deps.Xorg_libX11_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
-git-tree-sha1 = "afead5aba5aa507ad5a3bf01f58f82c8d1403495"
+git-tree-sha1 = "9dafcee1d24c4f024e7edc92603cedba72118283"
 uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
-version = "1.8.6+0"
+version = "1.8.6+3"
 
 [[deps.Xorg_libXau_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "6035850dcc70518ca32f012e46015b9beeda49d8"
+git-tree-sha1 = "e9216fdcd8514b7072b43653874fd688e4c6c003"
 uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
-version = "1.0.11+0"
+version = "1.0.12+0"
 
 [[deps.Xorg_libXdmcp_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "34d526d318358a859d7de23da945578e8e8727b7"
+git-tree-sha1 = "89799ae67c17caa5b3b5a19b8469eeee474377db"
 uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
-version = "1.1.4+0"
+version = "1.1.5+0"
 
 [[deps.Xorg_libXext_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"]
-git-tree-sha1 = "d2d1a5c49fae4ba39983f63de6afcbea47194e85"
+git-tree-sha1 = "d7155fea91a4123ef59f42c4afb5ab3b4ca95058"
 uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
-version = "1.3.6+0"
+version = "1.3.6+3"
 
 [[deps.Xorg_libXrender_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"]
-git-tree-sha1 = "47e45cd78224c53109495b3e324df0c37bb61fbe"
+git-tree-sha1 = "a490c6212a0e90d2d55111ac956f7c4fa9c277a6"
 uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
-version = "0.9.11+0"
+version = "0.9.11+1"
 
 [[deps.Xorg_libpthread_stubs_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "8fdda4c692503d44d04a0603d9ac0982054635f9"
+git-tree-sha1 = "c57201109a9e4c0585b208bb408bc41d205ac4e9"
 uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
-version = "0.1.1+0"
+version = "0.1.2+0"
 
 [[deps.Xorg_libxcb_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
-git-tree-sha1 = "bcd466676fef0878338c61e655629fa7bbc69d8e"
+git-tree-sha1 = "1a74296303b6524a0472a8cb12d3d87a78eb3612"
 uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
-version = "1.17.0+0"
+version = "1.17.0+3"
 
 [[deps.Xorg_xtrans_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "e92a1a012a10506618f10b7047e478403a046c77"
+git-tree-sha1 = "6dba04dbfb72ae3ebe5418ba33d087ba8aa8cb00"
 uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
-version = "1.5.0+0"
+version = "1.5.1+0"
 
 [[deps.Zlib_jll]]
 deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
 version = "1.2.13+1"
 
+[[deps.ZygoteRules]]
+deps = ["ChainRulesCore", "MacroTools"]
+git-tree-sha1 = "434b3de333c75fc446aa0d19fc394edafd07ab08"
+uuid = "700de1a5-db45-46bc-99cf-38207098b444"
+version = "0.2.7"
+
 [[deps.libaom_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "1827acba325fdcdf1d2647fc8d5301dd9ba43a9d"
+git-tree-sha1 = "522c1df09d05a71785765d19c9524661234738e9"
 uuid = "a4ae2306-e953-59d6-aa16-d00cac43593b"
-version = "3.9.0+0"
+version = "3.11.0+0"
 
 [[deps.libass_jll]]
 deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Zlib_jll"]
@@ -1621,9 +1679,9 @@ version = "2.0.3+0"
 
 [[deps.libpng_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"]
-git-tree-sha1 = "b70c870239dc3d7bc094eb2d6be9b73d27bef280"
+git-tree-sha1 = "d7b5bbf1efbafb5eca466700949625e07533aff2"
 uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
-version = "1.6.44+0"
+version = "1.6.45+1"
 
 [[deps.libvorbis_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
diff --git a/docs/src/API/GaussianProcess.md b/docs/src/API/GaussianProcess.md
index 8778c835..6290887c 100644
--- a/docs/src/API/GaussianProcess.md
+++ b/docs/src/API/GaussianProcess.md
@@ -10,11 +10,11 @@ PredictionType
 GaussianProcess
 GaussianProcess(
     ::GPPkg;
-    ::Union{K, KPy, Nothing},
+    ::Union{K, KPy, AK, Nothing},
     ::Any,
     ::FT,
     ::PredictionType,
-) where {GPPkg <: GaussianProcessesPackage, K <: Kernel, KPy <: PyObject, FT <: AbstractFloat}
+) where {GPPkg <: GaussianProcessesPackage, K <: GaussianProcesses.Kernel, KPy <: PyObject, AK <:AbstractGPs.Kernel, FT <: AbstractFloat}
 build_models!(::GaussianProcess{GPJL}, ::PairedDataContainer{FT}) where {FT <: AbstractFloat}
 optimize_hyperparameters!(::GaussianProcess{GPJL})
 predict(::GaussianProcess{GPJL},  ::AbstractMatrix{FT}) where {FT <: AbstractFloat}
diff --git a/docs/src/API/MarkovChainMonteCarlo.md b/docs/src/API/MarkovChainMonteCarlo.md
index 4580ebf4..f6729b2d 100644
--- a/docs/src/API/MarkovChainMonteCarlo.md
+++ b/docs/src/API/MarkovChainMonteCarlo.md
@@ -22,8 +22,12 @@ MCMC sampling.
 
 ```@docs
 MCMCProtocol
+AutodiffProtocol
+ForwardDiffProtocol
+ReverseDiffProtocol
 RWMHSampling
 pCNMHSampling
+BarkerSampling
 MetropolisHastingsSampler
 ```
 
diff --git a/src/Emulator.jl b/src/Emulator.jl
index 186b290f..1904fd9b 100644
--- a/src/Emulator.jl
+++ b/src/Emulator.jl
@@ -35,7 +35,7 @@ include("RandomFeature.jl")
 function throw_define_mlt()
     throw(ErrorException("Unknown MachineLearningTool defined, please use a known implementation"))
 end
-function build_models!(mlt, iopairs)
+function build_models!(mlt, iopairs, mlt_kwargs...)
     throw_define_mlt()
 end
 function optimize_hyperparameters!(mlt)
@@ -104,6 +104,7 @@ function Emulator(
     standardize_outputs_factors::Union{AbstractVector{FT}, Nothing} = nothing,
     decorrelate::Bool = true,
     retained_svd_frac::FT = 1.0,
+    mlt_kwargs...,
 ) where {FT <: AbstractFloat}
 
     # For Consistency checks
@@ -150,7 +151,7 @@ function Emulator(
         training_pairs = PairedDataContainer(training_inputs, decorrelated_training_outputs)
         # [4.] build an emulator
 
-        build_models!(machine_learning_tool, training_pairs)
+        build_models!(machine_learning_tool, training_pairs; mlt_kwargs...)
     else
         if decorrelate || !isa(machine_learning_tool, VectorRandomFeatureInterface)
             throw(ArgumentError("$machine_learning_tool is incompatible with option Emulator(...,decorrelate = false)"))
@@ -158,7 +159,7 @@ function Emulator(
         decomposition = nothing
         training_pairs = PairedDataContainer(training_inputs, training_outputs)
         # [4.] build an emulator - providing the noise covariance as a Tikhonov regularizer
-        build_models!(machine_learning_tool, training_pairs, regularization_matrix = obs_noise_cov)
+        build_models!(machine_learning_tool, training_pairs; regularization_matrix = obs_noise_cov, mlt_kwargs...)
     end
 
     return Emulator{FT}(
@@ -192,11 +193,11 @@ Default is to predict in the decorrelated space.
 """
 function predict(
     emulator::Emulator{FT},
-    new_inputs::AbstractMatrix{FT};
+    new_inputs::AM;
     transform_to_real = false,
     vector_rf_unstandardize = true,
     mlt_kwargs...,
-) where {FT <: AbstractFloat}
+) where {FT <: AbstractFloat, AM <: AbstractMatrix}
     # Check if the size of new_inputs is consistent with the GP model's input
     # dimension. 
     input_dim, output_dim = size(emulator.training_pairs, 1)
diff --git a/src/GaussianProcess.jl b/src/GaussianProcess.jl
index c8980170..6a59669b 100644
--- a/src/GaussianProcess.jl
+++ b/src/GaussianProcess.jl
@@ -1,10 +1,12 @@
 
-import GaussianProcesses: predict
-using GaussianProcesses
-
 using EnsembleKalmanProcesses.DataContainers
-
 using DocStringExtensions
+
+# [1] For GaussianProcesses
+import GaussianProcesses: predict, get_params, get_param_names
+using GaussianProcesses
+
+# [2] For SciKitLearn
 using PyCall
 using ScikitLearn
 const pykernels = PyNULL()
@@ -14,23 +16,33 @@ function __init__()
     copy!(pyGP, pyimport_conda("sklearn.gaussian_process", "scikit-learn=1.5.1"))
 end
 
+# [3] For AbstractGPs
+using AbstractGPs
+using KernelFunctions
+
 #exports (from Emulator)
 export GaussianProcess
 
-export GPJL, SKLJL
+export GPJL, SKLJL, AGPJL
 export YType, FType
 
+export get_params
+export get_param_names
+
+
 """
 $(DocStringExtensions.TYPEDEF)
 
 Type to dispatch which GP package to use:
 
- - `GPJL` for GaussianProcesses.jl,
- - `SKLJL` for the ScikitLearn GaussianProcessRegressor.
+ - `GPJL` for GaussianProcesses.jl, [julia - gradient-free only]
+ - `SKLJL` for the ScikitLearn GaussianProcessRegressor, [python - gradient-free]
+ - `AGPJL` for AbstractGPs.jl, [julia - ForwardDiff compatible]
 """
 abstract type GaussianProcessesPackage end
 struct GPJL <: GaussianProcessesPackage end
 struct SKLJL <: GaussianProcessesPackage end
+struct AGPJL <: GaussianProcessesPackage end
 
 """
 $(DocStringExtensions.TYPEDEF)
@@ -55,9 +67,9 @@ $(DocStringExtensions.TYPEDFIELDS)
 """
 struct GaussianProcess{GPPackage, FT} <: MachineLearningTool
     "The Gaussian Process (GP) Regression model(s) that are fitted to the given input-data pairs."
-    models::Vector{Union{<:GaussianProcesses.GPE, <:PyObject, Nothing}}
+    models::Vector{Union{<:GaussianProcesses.GPE, <:PyObject, <:AbstractGPs.PosteriorGP, Nothing}}
     "Kernel object."
-    kernel::Union{<:Kernel, <:PyObject, Nothing}
+    kernel::Union{<:GaussianProcesses.Kernel, <:PyObject, <:AbstractGPs.Kernel, Nothing}
     "Learn the noise with the White Noise kernel explicitly?"
     noise_learn::Bool
     "Additional observational or regularization noise in used in GP algorithms"
@@ -79,14 +91,20 @@ $(DocStringExtensions.TYPEDSIGNATURES)
 """
 function GaussianProcess(
     package::GPPkg;
-    kernel::Union{K, KPy, Nothing} = nothing,
+    kernel::Union{K, KPy, AGPK, Nothing} = nothing,
     noise_learn = true,
     alg_reg_noise::FT = 1e-3,
     prediction_type::PredictionType = YType(),
-) where {GPPkg <: GaussianProcessesPackage, K <: Kernel, KPy <: PyObject, FT <: AbstractFloat}
+) where {
+    GPPkg <: GaussianProcessesPackage,
+    K <: GaussianProcesses.Kernel,
+    KPy <: PyObject,
+    AGPK <: AbstractGPs.Kernel,
+    FT <: AbstractFloat,
+}
 
     # Initialize vector for GP models
-    models = Vector{Union{<:GaussianProcesses.GPE, <:PyObject, Nothing}}(undef, 0)
+    models = Vector{Union{<:GaussianProcesses.GPE, <:PyObject, <:AbstractGPs.PosteriorGP, <:Nothing}}(undef, 0)
 
     # the algorithm regularization noise is set to some small value if we are learning noise, else
     # it is fixed to the correct value (1.0)
@@ -98,6 +116,21 @@ function GaussianProcess(
 end
 
 # First we create  the GPJL implementation
+"""
+Gets flattened kernel hyperparameters from a (vector of) `GaussianProcess{GPJL}` model(s). Extends GaussianProcess.jl method.
+"""
+function GaussianProcesses.get_params(gp::GaussianProcess{GPJL})
+    return [get_params(model.kernel) for model in gp.models]
+end
+
+"""
+Gets the flattened names of kernel hyperparameters from a (vector of) `GaussianProcess{GPJL}` model(s). Extends GaussianProcess.jl method.
+"""
+function GaussianProcesses.get_param_names(gp::GaussianProcess{GPJL})
+    return [get_param_names(model.kernel) for model in gp.models]
+end
+
+
 """
 $(DocStringExtensions.TYPEDSIGNATURES)
 
@@ -105,7 +138,8 @@ Method to build Gaussian process models based on the package.
 """
 function build_models!(
     gp::GaussianProcess{GPJL},
-    input_output_pairs::PairedDataContainer{FT},
+    input_output_pairs::PairedDataContainer{FT};
+    kwargs...,
 ) where {FT <: AbstractFloat}
     # get inputs and outputs 
     input_values = get_inputs(input_output_pairs)
@@ -224,7 +258,8 @@ predict(gp::GaussianProcess{GPJL}, new_inputs::AbstractMatrix{FT}) where {FT <:
 #now we build the SKLJL implementation
 function build_models!(
     gp::GaussianProcess{SKLJL},
-    input_output_pairs::PairedDataContainer{FT},
+    input_output_pairs::PairedDataContainer{FT};
+    kwargs...,
 ) where {FT <: AbstractFloat}
     # get inputs and outputs 
     input_values = permutedims(get_inputs(input_output_pairs), (2, 1))
@@ -296,3 +331,108 @@ function predict(gp::GaussianProcess{SKLJL}, new_inputs::AbstractMatrix{FT}) whe
 
     return μ, σ2
 end
+
+
+#We build the AGPJL implementation
+function build_models!(
+    gp::GaussianProcess{AGPJL},
+    input_output_pairs::PairedDataContainer{FT};
+    kernel_params = nothing,
+    kwargs...,
+) where {FT <: AbstractFloat}
+    # get inputs and outputs
+    input_values = permutedims(get_inputs(input_output_pairs), (2, 1))
+    output_values = get_outputs(input_output_pairs)
+
+    # Number of models (We are fitting one model per output dimension, as data is decorrelated)
+    models = gp.models
+    if length(gp.models) > 0 # check to see if gp already contains models
+        @warn "GaussianProcess already built. skipping..."
+        return
+    end
+
+    ##############################################################################
+    # Notes on borrowing hyperparameters optimised within GPJL:
+    # optimisation of the GPJL with default kernels produces kernel parameters
+    # in the way of [a b c], where:
+    # c is the log_const_value
+    # [a b] is the rbf_len: lengthscale parameters for SEArd kernel
+    # const_value = exp.(2 .* log_const_value)
+    ##############################################################################
+    ## For example A 2D->2D sinusoid input example:
+    #=
+    log_const_value = [2.9031145778344696; 3.8325906110973795]
+    rbf_len = [1.9952706691900783 3.066374123568536; 5.783676639895112 2.195849064147456]
+    =#
+    if isnothing(kernel_params)
+        throw(ArgumentError("""
+AbstractGP currently does not (yet) learn hyperparameters internally. The following can be performed instead:
+1. Create and optimize a GPJL emulator and default kernel. (here called gp_jl)
+2. Create the Kernel parameters as a vect-of-dict with
+   kernel_params = [
+       Dict(
+           "log_rbf_len" => model_params[1:end-2] # input-dim Vector,
+           "log_std_sqexp" => model_params[end-2] # Float,
+           "log_std_noise" => # Float,
+       )
+    for model_params in get_params(gp_jl)]
+    Note: get_params(gp_jl) returns `output_dim`-vector where each entry is [a, b, c] with:
+    - a is the `rbf_len`: lengthscale parameters for SEArd kernel [input_dim] Vector
+    - b is the `log_std_sqexp` of the SQexp kernel Float
+    - c is the `log_std_noise` of the noise kernel Float
+3. Build a new Emulator with kwargs `kernel_params=kernel_params`
+        """))
+    end
+
+
+    N_models = size(output_values, 1) #size(transformed_data)[1]
+    regularization_noise = gp.alg_reg_noise
+
+    # now obtain the values of the hyperparameters
+    if N_models == 1 && !(isa(kernel_params, AbstractVector)) # i.e. just a Dict
+        kernel_params_vec = [kernel_params]
+    else
+        kernel_params_vec = kernel_params
+    end
+
+    for i in 1:N_models
+        var_sqexp = exp.(2 .* kernel_params_vec[i]["log_std_sqexp"]) # Float
+        var_noise = exp.(2 .* kernel_params_vec[i]["log_std_noise"]) # Float
+        rbf_invlen = 1 ./ exp.(kernel_params_vec[i]["log_rbf_len"])# Vec
+
+        opt_kern =
+            var_sqexp * (KernelFunctions.SqExponentialKernel() ∘ ARDTransform(rbf_invlen[:])) +
+            var_noise * KernelFunctions.WhiteKernel()
+        opt_f = AbstractGPs.GP(opt_kern)
+        opt_fx = opt_f(input_values', regularization_noise)
+
+        data_i = output_values[i, :]
+        opt_post_fx = posterior(opt_fx, data_i)
+        println("optimised GP: ", i)
+        push!(models, opt_post_fx)
+        println(opt_post_fx.prior.kernel)
+    end
+
+end
+
+function optimize_hyperparameters!(gp::GaussianProcess{AGPJL}, args...; kwargs...)
+    @info "AbstractGP already built. Continuing..."
+end
+
+function predict(gp::GaussianProcess{AGPJL}, new_inputs::AM) where {AM <: AbstractMatrix}
+
+    N_models = length(gp.models)
+    N_samples = size(new_inputs, 2)
+    FTorD = eltype(new_inputs) # e.g. Float or Dual
+    μ = zeros(FTorD, N_models, N_samples)
+    σ2 = zeros(FTorD, N_models, N_samples)
+    for i in 1:N_models
+        pred_gp = gp.models[i]
+        pred = pred_gp(new_inputs)
+        μ[i, :] = mean(pred)
+        σ2[i, :] = var(pred)
+    end
+
+    σ2[:, :] .= σ2[:, :] .+ gp.alg_reg_noise
+    return μ, σ2
+end
diff --git a/src/MarkovChainMonteCarlo.jl b/src/MarkovChainMonteCarlo.jl
index dd84f8be..64faa5d5 100644
--- a/src/MarkovChainMonteCarlo.jl
+++ b/src/MarkovChainMonteCarlo.jl
@@ -11,6 +11,8 @@ using LinearAlgebra
 using Printf
 using Random
 using Statistics
+using ForwardDiff
+using ReverseDiff
 
 using MCMCChains
 import AbstractMCMC: sample # Reexport sample()
@@ -20,13 +22,18 @@ import AdvancedMH
 export EmulatorPosteriorModel,
     MetropolisHastingsSampler,
     MCMCProtocol,
+    GradFreeProtocol,
+    ForwardDiffProtocol,
+    ReverseDiffProtocol,
     RWMHSampling,
     pCNMHSampling,
+    BarkerSampling,
     MCMCWrapper,
     accept_ratio,
     optimize_stepsize,
     get_posterior,
-    sample
+    sample,
+    esjd
 
 # ------------------------------------------------------------------------------------------
 # Output space transformations between original and SVD-decorrelated coordinates.
@@ -68,6 +75,47 @@ end
 # only change is made. We do the latter here because doing the former would require more
 # boilerplate code (repeating AdvancedMH/src/proposal.jl for the new Proposals)).
 
+# some possible types of autodiff used
+"""
+$(DocStringExtensions.TYPEDEF)
+
+Type used to dispatch different autodifferentiation methods where different emulators have a different compatability with autodiff packages 
+"""
+abstract type AutodiffProtocol end
+
+"""
+$(DocStringExtensions.TYPEDEF)
+
+Type to construct samplers for emulators not compatible with autodifferentiation
+"""
+abstract type GradFreeProtocol <: AutodiffProtocol end
+
+"""
+$(DocStringExtensions.TYPEDEF)
+
+Type to construct samplers for emulators compatible with `ForwardDiff.jl` autodifferentiation
+"""
+abstract type ForwardDiffProtocol <: AutodiffProtocol end
+
+"""
+$(DocStringExtensions.TYPEDEF)
+
+Type to construct samplers for emulators compatible with `ReverseDiff.jl` autodifferentiation
+"""
+abstract type ReverseDiffProtocol <: AutodiffProtocol end
+# ...to be implemented...
+#=
+abstract type ZygoteProtocol <: AutodiffProtocol end
+abstract type EnzymeProtocol <: AutodiffProtocol end
+=#
+
+function _get_proposal(prior::ParameterDistribution)
+    # *only* use covariance of prior, not full distribution
+    Σsqrt = sqrt(ParameterDistributions.cov(prior)) # rt_cov * MVN(0,I) avoids the posdef errors for MVN in Julia Distributions   
+    return AdvancedMH.RandomWalkProposal(Σsqrt * MvNormal(zeros(size(Σsqrt)[1]), I))
+end
+
+
 """
 $(DocStringExtensions.TYPEDEF)
 
@@ -82,10 +130,12 @@ $(DocStringExtensions.TYPEDEF)
 [`MCMCProtocol`](@ref) which uses Metropolis-Hastings sampling that generates proposals for 
 new parameters as as vanilla random walk, based on the covariance of `prior`.
 """
-struct RWMHSampling <: MCMCProtocol end
+struct RWMHSampling{T <: AutodiffProtocol} <: MCMCProtocol end
 
-struct RWMetropolisHastings{D} <: AdvancedMH.MHSampler
-    proposal::D
+RWMHSampling() = RWMHSampling{GradFreeProtocol}()
+
+struct RWMetropolisHastings{PT, ADT <: AutodiffProtocol} <: AdvancedMH.MHSampler
+    proposal::PT
 end
 # Define method needed by AdvancedMH for new Sampler
 AdvancedMH.logratio_proposal_density(
@@ -94,12 +144,6 @@ AdvancedMH.logratio_proposal_density(
     candidate,
 ) = AdvancedMH.logratio_proposal_density(sampler.proposal, transition_prev.params, candidate)
 
-function _get_proposal(prior::ParameterDistribution)
-    # *only* use covariance of prior, not full distribution
-    Σsqrt = sqrt(ParameterDistributions.cov(prior)) # rt_cov * MVN(0,I) avoids the posdef errors for MVN in Julia Distributions   
-    return AdvancedMH.RandomWalkProposal(Σsqrt * MvNormal(zeros(size(Σsqrt)[1]), I))
-end
-
 """
 $(DocStringExtensions.TYPEDSIGNATURES)
 
@@ -115,7 +159,10 @@ Constructor for all `Sampler` objects, with one method for each supported MCMC a
     (possibly in a transformed space) and we assume enough fidelity in the Emulator that 
     inference isn't prior-dominated.
 """
-MetropolisHastingsSampler(::RWMHSampling, prior::ParameterDistribution) = RWMetropolisHastings(_get_proposal(prior))
+function MetropolisHastingsSampler(::RWMHSampling{T}, prior::ParameterDistribution) where {T <: AutodiffProtocol}
+    proposal = _get_proposal(prior)
+    return RWMetropolisHastings{typeof(proposal), T}(proposal)
+end
 
 """
 $(DocStringExtensions.TYPEDEF)
@@ -125,9 +172,10 @@ new parameters according to the preconditioned Crank-Nicholson (pCN) algorithm,
 usable for MCMC in the *stepsize → 0* limit, unlike the vanilla random walk. Steps are based 
 on the covariance of `prior`.
 """
-struct pCNMHSampling <: MCMCProtocol end
+struct pCNMHSampling{T <: AutodiffProtocol} <: MCMCProtocol end
+pCNMHSampling() = pCNMHSampling{GradFreeProtocol}()
 
-struct pCNMetropolisHastings{D} <: AdvancedMH.MHSampler
+struct pCNMetropolisHastings{D, T <: AutodiffProtocol} <: AdvancedMH.MHSampler
     proposal::D
 end
 # Define method needed by AdvancedMH for new Sampler
@@ -136,8 +184,78 @@ AdvancedMH.logratio_proposal_density(
     transition_prev::AdvancedMH.AbstractTransition,
     candidate,
 ) = AdvancedMH.logratio_proposal_density(sampler.proposal, transition_prev.params, candidate)
+function MetropolisHastingsSampler(::pCNMHSampling{T}, prior::ParameterDistribution) where {T <: AutodiffProtocol}
+    proposal = _get_proposal(prior)
+    return pCNMetropolisHastings{typeof(proposal), T}(proposal)
+end
+
+#------ The following are gradient-based samplers
+
+"""
+$(DocStringExtensions.TYPEDEF)
+
+[`MCMCProtocol`](@ref) which uses Metropolis-Hastings sampling that generates proposals for
+new parameters according to the Barker proposal.
+"""
+struct BarkerSampling{T <: AutodiffProtocol} <: MCMCProtocol end
+BarkerSampling() = BarkerSampling{ForwardDiffProtocol}()
+
+struct BarkerMetropolisHastings{D, T <: AutodiffProtocol} <: AdvancedMH.MHSampler
+    proposal::D
+end
+# Define method needed by AdvancedMH for new Sampler
+AdvancedMH.logratio_proposal_density(
+    sampler::BarkerMetropolisHastings,
+    transition_prev::AdvancedMH.AbstractTransition,
+    candidate,
+) = AdvancedMH.logratio_proposal_density(sampler.proposal, transition_prev.params, candidate)
+
+function MetropolisHastingsSampler(::BarkerSampling{T}, prior::ParameterDistribution) where {T <: AutodiffProtocol}
+    proposal = _get_proposal(prior)
+    return BarkerMetropolisHastings{typeof(proposal), T}(proposal)
+end
+
+
+## -------- Autodifferentiation procedures ------ ##
+
+function autodiff_gradient(model::AdvancedMH.DensityModel, params, autodiff_protocol)
+    if autodiff_protocol == ForwardDiffProtocol
+        return ForwardDiff.gradient(x -> AdvancedMH.logdensity(model, x), params)
+    elseif autodiff_protocol == ReverseDiffProtocol
+        return ReverseDiff.gradient(x -> AdvancedMH.logdensity(model, x), params)
+    else
+        throw(
+            ArgumentError(
+                "Calling `autodiff_gradient(...)` on a sampler with protocol $(autodiff_protocol) that has *no* gradient implementation.\n Please select from a protocol with a gradient implementation (e.g., `ForwardDiffProtocol`).",
+            ),
+        )
+    end
+
+end
+
+autodiff_gradient(model::AdvancedMH.DensityModel, params, sampler::MH) where {MH <: AdvancedMH.MHSampler} =
+    autodiff_gradient(model::AdvancedMH.DensityModel, params, typeof(sampler).parameters[2]) # hacky way of getting the "AutodiffProtocol"
+
+function autodiff_hessian(model::AdvancedMH.DensityModel, params, autodiff_protocol)
+    if autodiff_protocol == ForwardDiffProtocol
+        return Symmetric(ForwardDiff.hessian(x -> AdvancedMH.logdensity(model, x), params))
+    elseif autodiff_protocol == ReverseDiffProtocol
+        return Symmetric(ReverseDiff.hessian(x -> AdvancedMH.logdensity(model, x), params))
+    else
+        throw(
+            ArgumentError(
+                "Calling `autodiff_hessian(...)` on a sampler with protocol $(autodiff_protocol) that has *no* hessian implementation.\n Please select from a protocol with a hessian implementation (e.g., `ForwardDiffProtocol`).",
+            ),
+        )
+
+    end
+end
+
+autodiff_hessian(model::AdvancedMH.DensityModel, params, sampler::MH) where {MH <: AdvancedMH.MHSampler} =
+    autodiff_hessian(model, params, typeof(sampler).parameters[2])
+
+
 
-MetropolisHastingsSampler(::pCNMHSampling, prior::ParameterDistribution) = pCNMetropolisHastings(_get_proposal(prior))
 
 # ------------------------------------------------------------------------------------------
 # Use emulated model in sampler
@@ -216,7 +334,7 @@ function AdvancedMH.transition(
     model::AdvancedMH.DensityModel,
     params,
     log_density::Real,
-) where {MHS <: Union{pCNMetropolisHastings, RWMetropolisHastings}}
+) where {MHS <: Union{pCNMetropolisHastings, RWMetropolisHastings, BarkerMetropolisHastings}}
     return MCMCState(params, log_density, true)
 end
 
@@ -245,6 +363,22 @@ function AdvancedMH.propose(
     return ρ * current_state.params .+ sqrt(1 - ρ^2) * rand(rng, sampler.proposal)
 end
 
+# method extending AdvancedMH.propose() for the Barker proposal
+function AdvancedMH.propose(
+    rng::Random.AbstractRNG,
+    sampler::BarkerMetropolisHastings,
+    model::AdvancedMH.DensityModel,
+    current_state::MCMCState;
+    stepsize::FT = 1.0,
+) where {FT <: AbstractFloat}
+    # Livingstone and Zanella (2022)
+    # Compute the gradient of the log-density at the current state
+    n = length(current_state.params)
+    log_gradient = autodiff_gradient(model, current_state.params, sampler)
+    xi = rand(rng, sampler.proposal)
+    return current_state.params .+ (stepsize .* ((rand(rng, n) .< 1 ./ (1 .+ exp.(-log_gradient .* xi))) .* xi))
+end
+
 # Copy a MCMCState and set accepted = false
 reject_transition(t::MCMCState) = MCMCState(t.params, t.log_density, false)
 
@@ -408,12 +542,14 @@ decorrelation) that was applied in the Emulator. It creates and wraps an instanc
     fixed stepsize.
   - [`pCNMHSampling`](@ref): Metropolis-Hastings sampling using the preconditioned 
     Crank-Nicholson algorithm, which has a well-behaved small-stepsize limit.
+  - [`BarkerSampling`](@ref): Metropolis-Hastings sampling using the Barker
+    proposal, which has a robustness to choosing step-size parameters.
 
 - `obs_sample`: A single sample from the observations. Can, e.g., be picked from an 
   Observation struct using `get_obs_sample`.
 - `prior`: [`ParameterDistribution`](https://clima.github.io/EnsembleKalmanProcesses.jl/dev/parameter_distributions/) 
   object containing the parameters' prior distributions.
-- `em`: [`Emulator`](@ref) to sample from. 
+- `emulator`: [`Emulator`](@ref) to sample from.
 - `stepsize`: MCMC step size, applied as a scaling to the prior covariance.
 - `init_params`: Starting parameter values for MCMC sampling.
 - `burnin`: Initial number of MCMC steps to discard from output (pre-convergence).
@@ -422,13 +558,13 @@ function MCMCWrapper(
     mcmc_alg::MCMCProtocol,
     obs_sample::AbstractVector{FT},
     prior::ParameterDistribution,
-    em::Emulator;
+    emulator::Emulator;
     init_params::AbstractVector{FT},
     burnin::IT = 0,
     kwargs...,
 ) where {FT <: AbstractFloat, IT <: Integer}
-    obs_sample = to_decorrelated(obs_sample, em)
-    log_posterior_map = EmulatorPosteriorModel(prior, em, obs_sample)
+    obs_sample = to_decorrelated(obs_sample, emulator)
+    log_posterior_map = EmulatorPosteriorModel(prior, emulator, obs_sample)
     mh_proposal_sampler = MetropolisHastingsSampler(mcmc_alg, prior)
 
     # parameter names are needed in every dimension in a MCMCChains object needed for diagnostics
@@ -541,6 +677,7 @@ function optimize_stepsize(
     init_stepsize = 1.0,
     N = 2000,
     max_iter = 20,
+    target_acc = 0.25,
     sample_kwargs...,
 )
     increase = false
@@ -555,9 +692,9 @@ function optimize_stepsize(
         _find_mcmc_step_log(it, stepsize, acc_ratio, trial_chain)
 
         change_step = true
-        if acc_ratio < 0.15
+        if acc_ratio < target_acc - 0.1
             decrease = true
-        elseif acc_ratio > 0.35
+        elseif acc_ratio > target_acc + 0.1
             increase = true
         else
             change_step = false
@@ -569,9 +706,9 @@ function optimize_stepsize(
             decrease = false
         end
 
-        if acc_ratio < 0.15
+        if acc_ratio < target_acc - 0.1
             stepsize *= 2^(-factor[1])
-        elseif acc_ratio > 0.35
+        elseif acc_ratio > target_acc + 0.1
             stepsize *= 2^(factor[1])
         end
 
@@ -623,4 +760,25 @@ function get_posterior(mcmc::MCMCWrapper, chain::MCMCChains.Chains)
     return posterior_distribution
 end
 
+# other diagnostic utilities
+"""
+$(DocStringExtensions.TYPEDSIGNATURES)
+
+Computes the expected squared jump distance of the chain.
+"""
+function esjd(chain::MCMCChains.Chains)
+    samples = chain.value[:, :, 1] # N_samples x N_params x n_chains
+    n_samples, n_params = size(samples)
+    esjd = zeros(Float64, n_params)
+    for i in 2:n_samples
+        esjd .+= (samples[i, :] .- samples[i - 1, :]) .^ 2 ./ n_samples
+    end
+    return esjd
+
+end
+
+
+
+
+
 end # module MarkovChainMonteCarlo
diff --git a/src/ScalarRandomFeature.jl b/src/ScalarRandomFeature.jl
index a946fd9c..2333a022 100644
--- a/src/ScalarRandomFeature.jl
+++ b/src/ScalarRandomFeature.jl
@@ -306,7 +306,8 @@ Builds the random feature method from hyperparameters. We use cosine activation
 """
 function build_models!(
     srfi::ScalarRandomFeatureInterface,
-    input_output_pairs::PairedDataContainer{FT},
+    input_output_pairs::PairedDataContainer{FT};
+    kwargs...,
 ) where {FT <: AbstractFloat}
     # get inputs and outputs 
     input_values = get_inputs(input_output_pairs)
diff --git a/src/VectorRandomFeature.jl b/src/VectorRandomFeature.jl
index 686cd1a2..cad00bd6 100644
--- a/src/VectorRandomFeature.jl
+++ b/src/VectorRandomFeature.jl
@@ -354,6 +354,7 @@ function build_models!(
     vrfi::VectorRandomFeatureInterface,
     input_output_pairs::PairedDataContainer{FT};
     regularization_matrix::Union{M, Nothing} = nothing,
+    kwargs...,
 ) where {FT <: AbstractFloat, M <: AbstractMatrix}
 
     # get inputs and outputs 
diff --git a/test/GaussianProcess/runtests.jl b/test/GaussianProcess/runtests.jl
index 755bf9c2..f92f2439 100644
--- a/test/GaussianProcess/runtests.jl
+++ b/test/GaussianProcess/runtests.jl
@@ -1,6 +1,6 @@
 # Import modules
-using Random
 using Test
+using Random
 using GaussianProcesses
 using Statistics
 using Distributions
@@ -80,6 +80,7 @@ using CalibrateEmulateSample.DataContainers
         standardize_outputs = false,
         retained_svd_frac = 1.0,
     ) # check that gp1 does not get more models added under second call
+
     Emulator(
         gp1,
         iopairs,
@@ -98,6 +99,57 @@ using CalibrateEmulateSample.DataContainers
     @test size(μ1) == (1, 5)
     @test vec(σ1²) ≈ [0.017, 0.003, 0.004, 0.004, 0.009] atol = 1e-2
 
+    # GaussianProcess 1b: use GPJL to create an abstractGP dist.
+    agp = GaussianProcess(AGPJL(); noise_learn = true, alg_reg_noise = 1e-4, prediction_type = pred_type)
+    @test_throws ArgumentError Emulator(
+        agp,
+        iopairs,
+        obs_noise_cov = nothing,
+        normalize_inputs = false,
+        standardize_outputs = false,
+        retained_svd_frac = 1.0,
+    )
+
+    gp1_opt_params = Emulators.get_params(gp1)[1] # one model only
+    gp1_opt_param_names = get_param_names(gp1)[1] # one model only
+
+    kernel_params = Dict(
+        "log_rbf_len" => gp1_opt_params[1:(end - 2)],
+        "log_std_sqexp" => gp1_opt_params[end - 1],
+        "log_std_noise" => gp1_opt_params[end],
+    )
+
+    em_agp_from_gp1 = Emulator(
+        agp,
+        iopairs,
+        obs_noise_cov = nothing,
+        normalize_inputs = false,
+        standardize_outputs = false,
+        retained_svd_frac = 1.0,
+        kernel_params = kernel_params,
+    )
+    optimize_hyperparameters!(em_agp_from_gp1)
+    # skip rebuild:
+    @test_logs (:warn,) (:warn,) Emulator(
+        agp,
+        iopairs,
+        obs_noise_cov = nothing,
+        normalize_inputs = false,
+        standardize_outputs = false,
+        retained_svd_frac = 1.0,
+        kernel_params = kernel_params,
+    )
+
+
+    μ1b, σ1b² = Emulators.predict(em_agp_from_gp1, new_inputs)
+
+    # gp1 and agp_from_gp2 should give similar predictions
+    tol_small = 1e-12
+    @test all(isapprox.(μ1, μ1b, atol = tol_small))
+    @test size(μ1) == (1, 5)
+    @test all(isapprox.(σ1², σ1b², atol = tol_small))
+
+
     # GaussianProcess 2: GPJL, predict_f
     pred_type = FType()
 
@@ -118,7 +170,6 @@ using CalibrateEmulateSample.DataContainers
     # predict_y and predict_f should give the same mean
     @test μ2 ≈ μ1 atol = 1e-6
 
-
     # GaussianProcess 3: SKLJL
 
     gppackage = SKLJL()
@@ -251,4 +302,34 @@ using CalibrateEmulateSample.DataContainers
     @test all(isapprox.(σ4²_noise_from_Σ, σ4²_noise_learnt, rtol = 2 * tol_mu))
 
 
+    # GaussianProcess 4b: use GPJL to create an abstractGP dist.
+    agp4 = GaussianProcess(AGPJL(); noise_learn = true, prediction_type = pred_type)
+
+    gp4_opt_params = Emulators.get_params(gp4_noise_learnt)
+    gp4_opt_param_names = get_param_names(gp4_noise_learnt)
+
+    kernel_params = [
+        Dict(
+            "log_rbf_len" => model_params[1:(end - 2)],
+            "log_std_sqexp" => model_params[end - 1],
+            "log_std_noise" => model_params[end],
+        ) for model_params in gp4_opt_params
+    ]
+
+    em_agp_from_gp4 = Emulator(
+        agp4,
+        iopairs2,
+        obs_noise_cov = Σ,
+        normalize_inputs = true,
+        retained_svd_frac = 1.0,
+        kernel_params = kernel_params,
+    )
+
+    μ4b, σ4b² = Emulators.predict(em_agp_from_gp4, new_inputs, transform_to_real = true)
+
+    # gp1 and agp_from_gp2 should give similar predictions
+    tol_small = 1e-12
+    @test all(isapprox.(μ4b, μ4_noise_learnt, atol = tol_small))
+    @test all(isapprox.(σ4b², σ4²_noise_learnt, atol = tol_small))
+
 end
diff --git a/test/MarkovChainMonteCarlo/runtests.jl b/test/MarkovChainMonteCarlo/runtests.jl
index e6534630..d24e7cb1 100644
--- a/test/MarkovChainMonteCarlo/runtests.jl
+++ b/test/MarkovChainMonteCarlo/runtests.jl
@@ -82,6 +82,54 @@ function test_gp_1(y, σ2_y, iopairs::PairedDataContainer; norm_factor = nothing
     return em
 end
 
+function test_gp_and_agp_1(y, σ2_y, iopairs::PairedDataContainer; norm_factor = nothing)
+    gppackage = GPJL()
+    pred_type = YType()
+    # Construct kernel:
+    # Squared exponential kernel (note that hyperparameters are on log scale)
+    # with observational noise
+    GPkernel = SE(log(1.0), log(1.0))
+    gp = GaussianProcess(gppackage; kernel = GPkernel, noise_learn = true, prediction_type = pred_type)
+    em = Emulator(
+        gp,
+        iopairs;
+        obs_noise_cov = σ2_y,
+        normalize_inputs = false,
+        standardize_outputs = false,
+        standardize_outputs_factors = norm_factor,
+        retained_svd_frac = 1.0,
+    )
+    Emulators.optimize_hyperparameters!(em)
+
+    # now make agp from gp
+    agp = GaussianProcess(AGPJL(); noise_learn = true, prediction_type = pred_type)
+
+    gp_opt_params = Emulators.get_params(gp)
+    gp_opt_param_names = get_param_names(gp)
+
+    kernel_params = [
+        Dict(
+            "log_rbf_len" => model_params[1:(end - 2)],
+            "log_std_sqexp" => model_params[end - 1],
+            "log_std_noise" => model_params[end],
+        ) for model_params in gp_opt_params
+    ]
+
+    em_agp = Emulator(
+        agp,
+        iopairs,
+        obs_noise_cov = σ2_y,
+        normalize_inputs = false,
+        standardize_outputs = false,
+        retained_svd_frac = 1.0,
+        standardize_outputs_factors = norm_factor,
+        kernel_params = kernel_params,
+    )
+
+    return em, em_agp
+end
+
+
 function test_gp_2(y, σ2_y, iopairs::PairedDataContainer; norm_factor = nothing)
     gppackage = GPJL()
     pred_type = YType()
@@ -100,6 +148,7 @@ function test_gp_2(y, σ2_y, iopairs::PairedDataContainer; norm_factor = nothing
         retained_svd_frac = 0.9,
     )
     Emulators.optimize_hyperparameters!(em)
+
     return em
 end
 
@@ -112,21 +161,21 @@ function mcmc_test_template(
     init_params = 3.0,
     step = 0.25,
     rng = Random.GLOBAL_RNG,
+    target_acc = 0.25,
 )
     obs_sample = reshape(collect(obs_sample), 1) # scalar or Vector -> Vector
     init_params = reshape(collect(init_params), 1) # scalar or Vector -> Vector
     mcmc = MCMCWrapper(mcmc_alg, obs_sample, prior, em; init_params = init_params)
-
     # First let's run a short chain to determine a good step size
-    new_step = optimize_stepsize(mcmc; init_stepsize = step, N = 5000)
+    new_step = optimize_stepsize(mcmc; init_stepsize = step, N = 5000, target_acc = target_acc)
 
     # Now begin the actual MCMC, sample is multiply exported so we qualify
-    chain = MCMC.sample(rng, mcmc, 100_000; stepsize = new_step, discard_initial = 1000)
+    chain = MCMC.sample(rng, mcmc, 100_000; stepsize = new_step, discard_initial = 50000)
     posterior_distribution = get_posterior(mcmc, chain)
     #post_mean = mean(posterior, dims=1)[1]
     posterior_mean = mean(posterior_distribution)
 
-    return new_step, posterior_mean[1]
+    return new_step, posterior_mean[1], chain
 end
 
 @testset "MarkovChainMonteCarlo" begin
@@ -136,6 +185,7 @@ end
     mcmc_params =
         Dict(:mcmc_alg => RWMHSampling(), :obs_sample => obs_sample, :init_params => [3.0], :step => 0.5, :rng => rng)
 
+
     @testset "Constructor: standardize" begin
         em = test_gp_1(y, σ2_y, iopairs)
         test_obs = MarkovChainMonteCarlo.to_decorrelated(obs_sample, em)
@@ -155,20 +205,39 @@ end
 
     end
 
+
     @testset "Sine GP & RW Metropolis" begin
-        em_1 = test_gp_1(y, σ2_y, iopairs)
-        new_step, posterior_mean_1 = mcmc_test_template(prior, σ2_y, em_1; mcmc_params...)
-        @test isapprox(new_step, 0.5; atol = 0.5)
+        em_1, em_1b = test_gp_and_agp_1(y, σ2_y, iopairs)
+        new_step_1, posterior_mean_1, chain_1 = mcmc_test_template(prior, σ2_y, em_1; mcmc_params...)
+        esjd1 = esjd(chain_1)
+        @info "ESJD = $esjd1"
+        @test isapprox(new_step_1, 0.5; atol = 0.5)
         # difference between mean_1 and ground truth comes from MCMC convergence and GP sampling
         @test isapprox(posterior_mean_1, π / 2; atol = 4e-1)
 
+        # test the agp setup on without derivatives
+        new_step_1b, posterior_mean_1b, chain_1b = mcmc_test_template(prior, σ2_y, em_1b; mcmc_params...)
+        @test isapprox(new_step_1b, new_step_1; atol = 1e-1)
+        tol_small = 5e-2
+        @test isapprox(posterior_mean_1b, posterior_mean_1; atol = tol_small)
+        esjd1b = esjd(chain_1b)
+        @info "ESJD = $esjd1b"
+        @test all(isapprox.(esjd1, esjd1b, rtol = 0.1))
+
         # now test SVD normalization
         norm_factor = 10.0
         norm_factor = fill(norm_factor, size(y[:, 1])) # must be size of output dim
         em_2 = test_gp_2(y, σ2_y, iopairs; norm_factor = norm_factor)
-        _, posterior_mean_2 = mcmc_test_template(prior, σ2_y, em_2; mcmc_params...)
+        _, posterior_mean_2, chain_2 = mcmc_test_template(prior, σ2_y, em_2; mcmc_params...)
         # difference between mean_1 and mean_2 only from MCMC convergence
         @test isapprox(posterior_mean_2, posterior_mean_1; atol = 0.1)
+        # test diagnostic functions on the chain
+        esjd2 = esjd(chain_2)
+        @info "ESJD = $esjd2"
+        # approx [0.04190683285347798, 0.1685296224916364, 0.4129400000002722]
+        @test all(isapprox.(esjd1, esjd2, rtol = 0.1))
+
+
     end
 
     @testset "Sine GP & pCN" begin
@@ -180,18 +249,70 @@ end
             :rng => rng,
         )
 
-        em_1 = test_gp_1(y, σ2_y, iopairs)
-        new_step, posterior_mean_1 = mcmc_test_template(prior, σ2_y, em_1; mcmc_params...)
-        @test isapprox(new_step, 0.75; atol = 0.6)
+        em_1, em_1b = test_gp_and_agp_1(y, σ2_y, iopairs)
+        new_step_1, posterior_mean_1, chain_1 = mcmc_test_template(prior, σ2_y, em_1; mcmc_params...)
+        esjd1 = esjd(chain_1)
+        @info "ESJD = $esjd1"
+        @test isapprox(new_step_1, 0.75; atol = 0.6)
         # difference between mean_1 and ground truth comes from MCMC convergence and GP sampling
         @test isapprox(posterior_mean_1, π / 2; atol = 4e-1)
 
+        # test the agp setup on without derivatives
+        new_step_1b, posterior_mean_1b, chain_1b = mcmc_test_template(prior, σ2_y, em_1b; mcmc_params...)
+        @test isapprox(new_step_1b, new_step_1; atol = 1e-1)
+        tol_small = 5e-2
+        @test isapprox(posterior_mean_1b, posterior_mean_1; atol = tol_small)
+        esjd1b = esjd(chain_1b)
+        @info "ESJD = $esjd1b"
+        @test all(isapprox.(esjd1, esjd1b, rtol = 0.1))
+
         # now test SVD normalization
         norm_factor = 10.0
         norm_factor = fill(norm_factor, size(y[:, 1])) # must be size of output dim
         em_2 = test_gp_2(y, σ2_y, iopairs; norm_factor = norm_factor)
-        _, posterior_mean_2 = mcmc_test_template(prior, σ2_y, em_2; mcmc_params...)
+        _, posterior_mean_2, chain_2 = mcmc_test_template(prior, σ2_y, em_2; mcmc_params...)
         # difference between mean_1 and mean_2 only from MCMC convergence
         @test isapprox(posterior_mean_2, posterior_mean_1; atol = 0.1)
+
+        esjd2 = esjd(chain_2)
+        @info "ESJD = $esjd2"
+        # approx [0.03470825350663073, 0.161606734823579, 0.38970000000024896]
+
+        @test all(isapprox.(esjd1, esjd2, rtol = 0.1))
+
     end
+
+
+
+
+    mcmc_params = Dict(:obs_sample => obs_sample, :init_params => [3.0], :step => 0.01, :rng => rng, :target_acc => 0.6) # the target is usually higher in grad-based MCMC
+
+    em_1, em_1b = test_gp_and_agp_1(y, σ2_y, iopairs)
+    # em_1 cannot be used here
+
+    mcmc_algs = [
+        RWMHSampling(), # sanity-check
+        BarkerSampling(), # ForwardDiffProtocol by default
+        BarkerSampling{ReverseDiffProtocol}(), # scales to high dim better, but slow.
+    ]
+
+    bad_mcmc_alg = BarkerSampling{GradFreeProtocol}()
+    @test_throws ArgumentError mcmc_test_template(prior, σ2_y, em_1; mcmc_alg = bad_mcmc_alg, mcmc_params...)
+
+
+    # GPJL doesnt support ForwardDiff
+    @test_throws ErrorException mcmc_test_template(prior, σ2_y, em_1; mcmc_alg = mcmc_algs[2], mcmc_params...)
+
+    for alg in mcmc_algs
+        @info "testing algorithm: $(typeof(alg))"
+        new_step, posterior_mean, chain = mcmc_test_template(prior, σ2_y, em_1b; mcmc_alg = alg, mcmc_params...)
+        esjd_tmp = esjd(chain)
+        @info "ESJD = $esjd_tmp"
+        @info posterior_mean
+        @testset "Sine GP & ForwardDiff variant:$(typeof(alg))" begin
+            @test isapprox(posterior_mean, π / 2; atol = 4e-1)
+        end
+
+    end
+
 end